Warped VRT of NETCDF file with geos projection

vincentsarago commented 5 years ago

Expected behavior and actual behavior.

I'm trying to use VRT to warp some GOES-16 files stored in GEOS projection

gdalinfo NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad
Driver: netCDF/Network Common Data Format
Files: OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc
Size is 10848, 10848
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unknown",
        DATUM["unknown",
            SPHEROID["Spheroid",6378137,298.2572221]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Geostationary_Satellite"],
    PARAMETER["central_meridian",-75],
    PARAMETER["satellite_height",35786023],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    EXTENSION["PROJ4","+proj=geos +lon_0=-75 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"]]

Direct warping of the data (gdalwarp -t_srs EPSG:4326 NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad out.tif) works fine but when trying to use VRT as temporary step I get IReadBlock failed at X offset 0, Y offset 0: Too many points (529 out of 529) failed to transform, unable to compute output bounds

Note: My real use case is for using rasterio's WarpedVRT but I figured I was getting the same error with GDAL directly.

Steps to reproduce the problem.

$ aws s3 cp s3://noaa-goes16/ABI-L1b-RadF/2019/295/14/OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc .

$ gdalwarp  -of VRT -t_srs EPSG:4326 NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad out.vrt
Creating output file that is 11173P x 10513L.
Processing NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad [1/1] : 0Using internal nodata values (e.g. 1023) for image NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad.
Copying nodata values from source NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad to destination out.vrt.
...10...20...30...40...50...60...70...80...90...100 - done.

$ gdal_translate out.vrt out.tif
Input file size is 11173, 10513
0ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds.
ERROR 1: out.vrt, band 1: IReadBlock failed at X offset 0, Y offset 0: Too many points (529 out of 529) failed to transform, unable to compute output bounds.

For example (please modify !!!): "gdalinfo myfile"

Operating system

Mac/Linux

GDAL version and provenance

gdal3.0.1 (sources) gdal2.4 (homebrew)

rouault commented 5 years ago

You probably need to define manually the target extent with -te. Geostationnary projection is hard to deal with since generally the corners of the image point to space...

vincentsarago commented 5 years ago

adding -te to gdalwarp didn't work

$ gdalwarp -t_srs EPSG:4326 NETCDF:OR_ABI-L1b-RadF-M6C01_G16_s20192951400357_e20192951410065_c20192951410137.nc:Rad out.vrt -of VRT -te -156.2813002427007 -76.40208979748526 6.118050388262418 76.40417290575606

$ gdal_translate out.vrt out.tif
Input file size is 11173, 10513
0ERROR 1: Too many points (529 out of 529) failed to transform, unable to compute output bounds.
ERROR 1: out.vrt, band 1: IReadBlock failed at X offset 0, Y offset 0: Too many points (529 out of 529) failed to transform, unable to compute output bounds.

rouault commented 5 years ago

Backported to 3.0 and 2.4 branches

OSGeo / gdal