[gdal-dev] Using geolocation arrays

[Laurențiu Nicola]

Hi,

I have a NetCDF-style product with a separate geolocation dataset I'm trying to convert to something more easy to work with. Not sure if it matters, but the geolocation arrays are scaled. I tried to unscale them first:

$ cat lat.vrt
<VRTDataset rasterXSize="4865" rasterYSize="4090">
  <VRTRasterBand dataType="Float32" band="1" blockXSize="973" blockYSize="818">
    <NoDataValue>-2147483648</NoDataValue>
    <Scale>1e-06</Scale>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">NetCDF:geo_coordinates.nc:latitude</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4865" RasterYSize="4090" DataType="Int32" BlockXSize="973" BlockYSize="818" />
      <SrcRect xOff="0" yOff="0" xSize="4865" ySize="4090" />
      <DstRect xOff="0" yOff="0" xSize="4865" ySize="4090" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

$ gdal_translate -unscale -ot Float32 lat.vrt lat.tif

# same for lon.{vrt,tif}

Then I made a VRT for the main dataset:

<VRTDataset rasterXSize="4865" rasterYSize="4090">
  <Metadata domain="GEOLOCATION">
    <MDI key="X_DATASET">lon.tif</MDI>
    <MDI key="X_BAND">1</MDI>
    <MDI key="Y_DATASET">lat.tif</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="PIXEL_OFFSET">0</MDI>
    <MDI key="LINE_OFFSET">0</MDI>
    <MDI key="PIXEL_STEP">1</MDI>
    <MDI key="LINE_STEP">1</MDI>
  </Metadata>
  <VRTRasterBand dataType="UInt16" band="1" blockXSize="1217" blockYSize="1023">
    <NoDataValue>65535</NoDataValue>
    <Scale>0.01394645962864161</Scale>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">Oa01_radiance.nc</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4865" RasterYSize="4090" DataType="UInt16" BlockXSize="1217" BlockYSize="1023" />
      <SrcRect xOff="0" yOff="0" xSize="4865" ySize="4090" />
      <DstRect xOff="0" yOff="0" xSize="4865" ySize="4090" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

The coordinates look fine in the two TIFFs:

$ gdalinfo -stats lat.tif
Driver: GTiff/GeoTIFF
Files: lat.tif
Size is 4865, 4090
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4090.0)
Upper Right ( 4865.0,    0.0)
Lower Right ( 4865.0, 4090.0)
Center      ( 2432.5, 2045.0)
Band 1 Block=4865x1 Type=Float32, ColorInterp=Gray
  Minimum=39.552, Maximum=52.456, Mean=46.179, StdDev=3.112

$ gdalinfo -stats lon.tif
Driver: GTiff/GeoTIFF
Files: lon.tif
Size is 4865, 4090
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4090.0)
Upper Right ( 4865.0,    0.0)
Lower Right ( 4865.0, 4090.0)
Center      ( 2432.5, 2045.0)
Band 1 Block=4865x1 Type=Float32, ColorInterp=Gray
  Minimum=16.249, Maximum=36.250, Mean=25.532, StdDev=4.917
  NoData Value=-2.1474836e+09
  Metadata:
    STATISTICS_MINIMUM=16.248615264893
    STATISTICS_MAXIMUM=36.24991607666
    STATISTICS_MEAN=25.532068327481
    STATISTICS_STDDEV=4.9167101776381
    STATISTICS_VALID_PERCENT=100

However, this fails with "Warning 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds.":

$ gdalwarp -geoloc -t_srs EPSG:4326 Oa01.vrt Oa01.tif

The files are available at https://nas.dend.ro/browse/datasets/Oa01.zip, if anyone can to take a look.

Thanks,
Laurentiu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20250926/38c4da80/attachment-0001.htm>