Hi all,
I'm being pulled out of retirement to fix something I created years ago.
I'm probably doing something stupid but I haven't been able to sort it
out, so ask for some help.
System: Ubuntu-22.04
GDAL 3.4.1, released 2021/12/27
OpenLayers 5.3.3
MapServer version 7.6.4 OUTPUT=PNG OUTPUT=JPEG OUTPUT=KML SUPPORTS=PROJ
SUPPORTS=AGG SUPPORTS=FREETYPE SUPPORTS=CAIRO SUPPORTS=SVG_SYMBOLS
SUPPORTS=RSVG SUPPORTS=ICONV SUPPORTS=FRIBIDI SUPPORTS=WMS_SERVER
SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT
SUPPORTS=WCS_SERVER SUPPORTS=SOS_SERVER SUPPORTS=FASTCGI
SUPPORTS=THREADS SUPPORTS=GEOS SUPPORTS=POINT_Z_M SUPPORTS=PBF
INPUT=JPEG INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE
Right some of this is pretty old, but I suspect that is not the current
issue.
1. I extract some data from HYCOM and generate a GTiff
2. Manually generate a VRT to add a color table
3. Have a mapserver mapfile to display it
Step 1. seems to work fine, but when I do a mapserver query I get no results
Step 2. seems to look good and if a gdal_translate -of PNG -outsize 1%
1% source target the image displays
Step 3. just gives me empty transparent tiles in OpenLayers
Here is the gdalinfo results for the GTif:
$ gdalinfo -stats -hist HYCOM_tomorrow_mld.tif
Driver: GTiff/GeoTIFF
Files: HYCOM_tomorrow_mld.tif
Size is 4500, 4251
Coordinate System is:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.080000000000000,-0.040000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left (-180.0000000, -80.0400000) (180d 0' 0.00"W, 80d 2'24.00"S)
Upper Right ( 180.0000000, 90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -80.0400000) (180d 0' 0.00"E, 80d 2'24.00"S)
Center ( 0.0000000, 4.9800000) ( 0d 0' 0.01"E, 4d58'48.00"N)
Band 1 Block=4500x1 Type=Float32, ColorInterp=Gray
Minimum=2.000, Maximum=5000.000, Mean=96.114, StdDev=346.786
0...10...20...30...40...50...60...70...80...90...100 - done.
256 buckets from -7.8 to 5009.8:
1489490 2306537 1907189 2028756 996825 306674 493493 0 196427 0 94515
0 0 33174 0 22423 0 0 11571 0 5892 0 0 0 0 8326 0 0 0 0 0 8045 0 0 0 0
7615 0 0 0 0 6418 0 0 0 0 5695 0 0 0 0 4944 0 0 0 0 0 0 0 0 0 0 0 0
13673 0 0 0 0 0 0 0 0 0 0 0 17543 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 31132 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 30109
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 40660 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 28594 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1702
NoData Value=-30000
Overviews: 2250x2126, 1125x1063, 563x532, 282x266, 141x133, 71x67
Metadata:
STATISTICS_MAXIMUM=5000
STATISTICS_MEAN=96.114400289494
STATISTICS_MINIMUM=2
STATISTICS_STDDEV=346.78571541402
STATISTICS_VALID_PERCENT=52.78
My assumption based on this is that I want to scale the 5000 to 250 in
the VRT which will be Type=Byte, So here is the VRT file.
$ cat HYCOM_tomorrow_mld.vrt <VRTDataset
rasterXSize="72000" rasterYSize="68016">
<SRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS
84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]</SRS>
<GeoTransform> -180.0 ,0.005 , 0 , 90.0 , 0 , -0.0025 </GeoTransform>
<Metadata>
</Metadata>
<VRTRasterBand dataType="Byte" band="1">
<ColorInterp>Palette</ColorInterp>
<ColorTable>
<Entry c1="0" c2="0" c3="0" c4="255"/>
[Snip lots of entries]
<Entry c1="255" c2="255" c3="255" c4="255"/>
</ColorTable>
<ComplexSource resampling="bilinear">
<SourceFilename
relativeToVRT="1">HYCOM_tomorrow_mld.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SrcRect xOff="0" yOff="0" xSize="4500" ySize="4251"/>
<DstRect xOff="0" yOff="0" xSize="72000" ySize="68016"/>
<ScaleRatio>1</ScaleRatio>
<NODATA>0</NODATA>
</ComplexSource>
</VRTRasterBand>
</VRTDataset>
So to scale the data I set ScaleRatio to 0.05, but the PNG image doesn't
look right and setting it to "1" does look correct. I'm also confused as
to what to set NODATA value to. One source said/implied that the NoData
of the source file would get set to NODATA value in the VRT, but seem to
be countered by the GDAL online docs.
Here is the gdalinfo for the VRT:
$ gdalinfo -stats -hist HYCOM_tomorrow_mld.vrt
Driver: VRT/Virtual Raster
Files: HYCOM_tomorrow_mld.vrt
HYCOM_tomorrow_mld.tif
Size is 72000, 68016
Coordinate System is:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.005000000000000,-0.002500000000000)
Corner Coordinates:
Upper Left (-180.0000000, 90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left (-180.0000000, -80.0400000) (180d 0' 0.00"W, 80d 2'24.00"S)
Upper Right ( 180.0000000, 90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -80.0400000) (180d 0' 0.00"E, 80d 2'24.00"S)
Center ( 0.0000000, 4.9800000) ( 0d 0' 0.01"E, 4d58'48.00"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Palette
Minimum=2.000, Maximum=5000.000, Mean=96.114, StdDev=346.786
0...10...20...30...40...50...60...70...80...90...100 - done.
256 buckets from -0.5 to 255.5:
0 0 534720 0 330663 0 149196 0 241797 0 233114 0 285113 0 0 385666 0
0 0 0 574491 0 0 0 0 517919 0 0 0 0 543348 0 0 0 0 455442 0 0 0 0 487763
0 0 0 0 530058 0 0 0 0 433926 0 0 0 0 0 0 0 0 0 1117668 0 0 0 0 0 0 0 0
0 911088 0 0 0 0 0 0 0 0 0 602087 0 0 0 0 0 0 0 0 0 394738 0 0 0 0 0 0 0
0 0 306674 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 493493 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 196427 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 94515 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 33174 0 0 0 0 244342
Overviews: 36000x34016, 18000x17008, 9008x8512, 4512x4256, 2256x2128,
1136x1072
Metadata:
STATISTICS_MAXIMUM=5000
STATISTICS_MEAN=96.114400289494
STATISTICS_MINIMUM=2
STATISTICS_STDDEV=346.78571541402
Color Table (RGB with 256 entries)
0: 0,0,0,255
[snip additional entries]
And finally the mapfile that is getting used via OpenLayers:
MAP
NAME "HYCOM_2d_2"
STATUS ON
SIZE 950 600
EXTENT -180 -90 180 90
UNITS DD
IMAGECOLOR "#000000"
IMAGETYPE agg_qn
MAXSIZE 4096
CONFIG ON_MISSING_DATA "IGNORE"
# CONFIG MS_ERRORFILE "stderr"
# DEBUG 10
OUTPUTFORMAT
NAME "agg_qn"
DRIVER "AGG/PNG"
EXTENSION "png"
MIMETYPE "image/png"
IMAGEMODE RGBA
FORMATOPTION "INTERLACE=false"
FORMATOPTION "QUANTIZE_NEW=ON"
FORMATOPTION "QUANTIZE_FORCE=ON"
FORMATOPTION "QUANTIZE_DITHER=OFF"
FORMATOPTION "QUANTIZE_COLORS=256"
TRANSPARENT ON
FORMATOPTION "TRANSPARENT=ON"
END
OUTPUTFORMAT
NAME aggpng24
DRIVER AGG/PNG
MIMETYPE "image/png"
IMAGEMODE RGB
EXTENSION "png"
END
OUTPUTFORMAT
NAME jpeg
DRIVER AGG/JPEG
MIMETYPE "image/jpeg"
IMAGEMODE RGB
EXTENSION "jpg"
FORMATOPTION "GAMMA=0.75"
FORMATOPTION "QUALITY=75"
END
WEB
METADATA
labelcache_map_edge_buffer "-20"
"ows_title" "iMaptools - HYCOM data"
"ows_onlineresource"
"http://map01.saltwatercentral.com/cgi-bin/mapserv?MAP=/maps/wms/hycom_2d_2.map";
"ows_srs" "EPSG:4326 EPSG:900913 EPSG:3857"
"ows_contactperson" "Stephen Woodbridge"
"ows_contactorganization" "iMaptools.com"
"ows_contactposition" "Owner"
"ows_contactelectronicmailaddress" "i...@imaptools.com"
"ows_enable_request" "GetMap"
"ows_http_max_age" "3200"
END
END
PROJECTION "init=epsg:4326" END
LAYER
NAME "ssh"
STATUS ON
TYPE RASTER
PROJECTION "init=epsg:4326" END
DATA "/maps/wms/data/HYCOM/HYCOM_tomorrow_ssh.vrt"
PROCESSING "NODATA=-30000"
PROCESSING "SCALE=-1.0,1.0"
END
LAYER
NAME "mlt"
STATUS ON
TYPE RASTER
PROJECTION "init=epsg:4326" END
DATA "/maps/wms/data/HYCOM/HYCOM_tomorrow_mlt.vrt"
#PROCESSING "NODATA=1.2676506002282294e+30"
PROCESSING "NODATA=-30000"
PROCESSING "SCALE=0,250"
END
LAYER
NAME "query"
STATUS ON
TYPE raster
PROJECTION "init=epsg:4326" END
VALIDATION 'type' '.' END
TOLERANCE 0
TOLERANCEUNITS pixels
DATA "/maps/wms/data/HYCOM/HYCOM_tomorrow_%type%.tif"
PROCESSING "NODATA=-30000"
TEMPLATE "/maps/wms/pixel.value.html"
END
END
I've tried various options of including PROCESSING "NODATA=' and
PROCESSING "SCALE=" to no success. And the "ssh" layer which is nearly
identical seems to would fine.
It get displayed via OpenLayers using a layer definition like:
hycom_mlt_2: new TileLayer({
title: 'Mixed Layer Depth (0.3degC chg)',
clickable: 'hycom_mlt_2',
type: 'overlay',
visible: false,
source: new TileWMS({
projection: 'EPSG:3857',
urls: getMapUrls('//', map_hosts,
'/cgi-bin/mapserv?map=/maps/wms/hycom_2d_2.map'),
params: {
'LAYERS': 'mlt',
'format': 'image/png',
'version': '1.3.0'
}
})
}),
I would appreciate any assistance in figuring out what I'm doing wrong.
Thanks,
Steve
--
This email has been checked for viruses by Avast antivirus software.
www.avast.com_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
