1

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>
  • Most likely the .vrt file is created wrongly. Can you add the build commandline? – AndreJ Jun 20 '16 at 13:59
  • See updated question – rakbar Jun 20 '16 at 14:11
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.

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.

  • 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? – rakbar Jun 20 '16 at 14:23
  • 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? – user30184 Jun 20 '16 at 14:45
  • I've updated the post to show how the tif is created and the output of gdalinfo. Thanks! – rakbar Jun 20 '16 at 14:50
  • Also, the *.h5 has no meta data on geolocation, just arrays of lats and longs. – rakbar Jun 20 '16 at 14:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.