Understanding GDAL ERROR 1: Too many points... failed to transform?

Question

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>

• AndreJ: Most likely the .vrt file is created wrongly. Can you add the build commandline?
• rakbar: See updated question

Answer 1

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.

• rakbar: Thank you for pointing that out! The corner coordinate values are apparently the dimensions of the rows and columns of the arrays (block size); not actual coordinates! How can I generate/extract the proper corners?
• user30184: The main tiff file should contain the georeferencing info or you should have received the metadata from website, as an attachment etc. What does gdalinfo /mnt/c/Users/20150509/D/sigma0_hh_fore/20150509_01439_D_sigma0_hh_fore.tifshow?
• rakbar: I've updated the post to show how the tif is created and the output of gdalinfo. Thanks!
• rakbar: Also, the *.h5 has no meta data on geolocation, just arrays of lats and longs.

Answer 2

In my answer to Unable to warp HDF5 files, I used this command line:

 gdalwarp -geoloc -t_srs EPSG:4326 chlor.vrt chlor-out.tif

This forces the use of the X_BAND and Y_BAND (which contain lonlat degree for every cell as cell values) in the VRT, and assigns EPSG:4326 to the output file.

Your commandline has -s_srs EPSG:4326 -t_srs '+proj=cea ..., which might kill the geoloc using the bands.

Try to do the reprojection to cea in a second step.