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";
<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"
<mailto: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
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
Virus-free.www.avast.com
<https://www.avast.com/sig-email?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
