[GIS] Understanding GDAL ERROR 1: Too many points… failed to transform

coordinate systemgdalwarpvrt

How can one trace the source of the following error when using gdalwarp? (I'm very new to GIS, GDAL and python!)

Warning 1: Unable to compute source region for output window 0,0,5789,4748, skipping.

ERROR 1: Too many points (439 out of 441) failed to transform,

These errors are output after using gdalwarp on a virtual raster (.vrt) with the following info (Updated entire info)

Driver: VRT/Virtual Raster
Files: 20150509_01439_D_sigma0_hh_fore_geo.vrt
       /mnt/c/Users/20150509/D/sigma0_hh_fore/20150509_01439_D_sigma0_hh_fore.tif
Size is 1091, 12698
Coordinate System is:
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"]]
Metadata:
  AREA_OR_POINT=Area
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  X_BAND=1
  X_DATASET=lon_20150509_01439_D_sigma0_hh_fore.vrt
  Y_BAND=1
  Y_DATASET=lat_20150509_01439_D_sigma0_hh_fore.vrt
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,12698.0)
Upper Right ( 1091.0,    0.0)
Lower Right ( 1091.0,12698.0)
Center      (  545.5, 6349.0)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  Description = 20150509_01439_D_sigma0_hh_fore
  Min=0.000 Max=9.821
  Minimum=0.000, Maximum=9.821, Mean=0.086, StdDev=0.107
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=9.8212261199951
    STATISTICS_MEAN=0.085568480608667
    STATISTICS_MINIMUM=1.6214236620016e-09
    STATISTICS_STDDEV=0.10704839226272

using the following tags

-geoloc -overwrite -co "COMPRESS=LZW" -tr 3000 3000 -srcnodata -9999 -dstnodata -9999 -r average -s_srs EPSG:4326 -t_srs '+proj=cea +x_0=0 +y_0=0 +lon_0=0 +lat_ts=30 +ellps=WGS84 +units=m +nodefs

Update and steps so far:

  1. open h5 file using h5py and performed some data filtering (numpy)

  2. saved numpy arrays to Tiff using the following commands

    srs = osr.SpatialReference()
    srs.SetWellKnownGeogCS("WGS84")
    nameDateOrbitBand = nameDateOrbit+"_"+inBandNames[i]
    print(' Creating Tiff for '+ nameDateOrbitBand)

    outFileName = nameDateOrbitBand +'.tif'
gdalDriver = getGDALFormatFromExt(outFileName)
format="GTiff"
driver = gdal.GetDriverByName(format)
metadata = driver.GetMetadata()
newDataset = driver.Create(outFileName, inXSize, inYSize, 1, gdal.GDT_Float32)
newDataset.SetProjection(srs.ExportToWkt())
temp_name = inBandNames[i]
temp_data = inDataLayers[i]

print(' Saving ' + temp_name)

try:
newDataset.GetRasterBand(1).SetDescription(nameDateOrbitBand)
newDataset.GetRasterBand(1).WriteArray(temp_data)
newDataset.GetRasterBand(1).SetNoDataValue(nodat)
stat = newDataset.GetRasterBand(1).ComputeStatistics(0)  # get the band statistics (min, max, mean, standard deviation)
newDataset.GetRasterBand(1).SetStatistics(stat[0], stat[1], stat[2], stat[3]) # set the stats we just got to the band

except KeyError:
print('Something went wrong\n')
except Exception as err:
print(err)

newDataset.FlushCache()
newDataset = None

  3. Followed similar steps given here (https://gis.stackexchange.com/a/154360/74531) to generate the .vrt files

lat and lon .vrt:

<VRTDataset rasterXSize="1091" rasterYSize="12698">
  <VRTRasterBand dataType="Float32" band="1">
    <Metadata>
      <MDI key="Sigma0_Data_cell_lat_long_name">Representative geodetic latitude of grid cell</MDI>
      <MDI key="Sigma0_Data_cell_lat_units">degrees_north</MDI>
      <MDI key="Sigma0_Data_cell_lat_valid_max">90 </MDI>
      <MDI key="Sigma0_Data_cell_lat_valid_min">-90 </MDI>
    </Metadata>
    <SimpleSource>
      <SourceFilename relativeToVRT="0">HDF5:20150509T231624_T12373_001.h5://Sigma0_Data/cell_lat</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1091" RasterYSize="12698" DataType="Float32" BlockXSize="1091" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="1091" ySize="12698" />
      <DstRect xOff="0" yOff="0" xSize="1091" ySize="12698" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

.vrt for data file:

<VRTDataset rasterXSize="1091" rasterYSize="12698">
  <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>
  <Metadata>
    <MDI key="AREA_OR_POINT">Area</MDI>
  </Metadata>
<Metadata domain="GEOLOCATION">
<MDI key="X_DATASET">lon_20150509_01439_D_sigma0_hh_fore.vrt</MDI>
<MDI key="X_BAND">1</MDI>
<MDI key="Y_DATASET">lat_20150509_01439_D_sigma0_hh_fore.vrt</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="Float32" band="1">
    <Metadata>
      <MDI key="STATISTICS_MAXIMUM">9.8212261199951</MDI>
      <MDI key="STATISTICS_MEAN">0.085568480608667</MDI>
      <MDI key="STATISTICS_MINIMUM">1.6214236620016e-09</MDI>
      <MDI key="STATISTICS_STDDEV">0.10704839226272</MDI>
    </Metadata>
    <Description>20150509_01439_D_sigma0_hh_fore</Description>
    <NoDataValue>-9.99900000000000E+03</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">20150509_01439_D_sigma0_hh_fore.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="1091" RasterYSize="12698" DataType="Float32" BlockXSize="1091" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="1091" ySize="12698" />
      <DstRect xOff="0" yOff="0" xSize="1091" ySize="12698" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Best Answer

The error itself means that GDAL fails when it tries to transform coordinates from the source srs into target srs. The source srs is supposed to be EPSG:4326 but by looking at the corner coordinates it can't be right.

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,12698.0)
Upper Right ( 1091.0,    0.0)
Lower Right ( 1091.0,12698.0)
Center      (  545.5, 6349.0)

For example the lower right corner at 1091.0,12698.0 is so much outsize the valid coordinate range of EPSG:4326 [-180,-90,180,90] that it is not a wonder that computing the transformation fails.

Before you can warp from one projection into another you must get the source projection right.

Related Question