0

I am projecting global rasters with the Homolosine. gdalwarp reprojects it nicely, but it does not compress the output, that must be done afterwards with gdal_translate. To reduce computation time, I tried applying gdalwarp into a VRT and then the gdal_translate into GeoTIFF.

The warping carries out fine:

$ gdalwarp -overwrite -s_srs EPSG:4326 -t_srs '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' -multi -co "BIGTIFF=YES" -te -20037500 -6729000 20037250 8600750 --config CHECK_WITH_INVERT_PROJ TRUE -tr 250 250 -of VRT input/TEST.tif output/TEST.tif.vrt
Creating output file that is 160299P x 61319L.
Processing input file input/TEST.tif.
Using internal nodata values (e.g. 65535) for image input/TEST.tif.
Copying nodata values from source input/TEST.tif to destination output/TEST.tif.vrt.

Checking with gdalinfo, everything looks in place:

$ gdalinfo output/TEST.tif.vrt
Driver: VRT/Virtual Raster
Files: output/TEST.tif.vrt
       input/TEST.tif
Size is 160299, 61319
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Interrupted_Goode_Homolosine"],
    UNIT["Meter",1]]
Origin = (-20037500.000000000000000,8600750.000000000000000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  SCALE=0.0001
Corner Coordinates:
Upper Left  (-20037500.000, 8600750.000) 
Lower Left  (-20037500.000,-6729000.000) 
Upper Right (20037250.000, 8600750.000) 
Lower Right (20037250.000,-6729000.000) 
Center      (    -125.000,  935875.000) (  0d19'37.24"W,  8d24'25.59"N)
Band 1 Block=512x128 Type=UInt16, ColorInterp=Gray
  NoData Value=65535

The final step with gdal_translate:

$ gdal_translate -co "COMPRESS=LZW" -co "PREDICTOR=2" -co "BIGTIFF=YES" -of GTiff output/TEST.tif.vrt output/TEST.tif
Input file size is 160299, 61319
0ERROR 1: Too many points (441 out of 441) failed to transform, unable to compute output bounds.
ERROR 1: output/TEST.tif.vrt, band 1: IReadBlock failed at X offset 0, Y offset 0
ERROR 1: GetBlockRef failed at X block offset 0, Y block offset 0

Is this behaviour expected? Could it be a bug?

  • 1
    I guess this is not releted with your problem at all, but if you do not have a special reason to create striped tiff with 160299 pixel wide stripes I would recommend to use creation option -co TILED=YES. – user30184 Feb 25 at 11:25
  • I haven't ever had problems compressing global rasters when projecting into Goode Homosline using gdalwarp. Is it giving you an error, or just not compressing the output? What is the command (with -co compress= option) you're feeding to gdalwarp? As for the second issue, I find that re-projecting to Goode Homosline often results in exactly this behaviour (points failing to transform), due to the fact that some portions of an image are just not transformable into Goode (at the global scale). One work around is to also include some [warp options](gdal.org/structGDALWarpOptions.html – scabecks Feb 26 at 4:25
  • @scabecks What leads you to think the warp is failing? It produces correct results into GeoTIFF or other hard format. Furthermore, the meta-data in the VRT created by the warp matches exactly that of a GeoTIFF output. Could you also elaborate on the claim that "some portions of the image are not transformable to the Homolosine"? Which portions are those? – Luís de Sousa Feb 26 at 11:05
  • @Luis de Sousa The warp isn't failing (when the re-projection is actually happening in the gdal_translate step), but there may be parts of the raster that fail to transform and are then not successfully written to the output raster - the result being that blocks of cells where there should be data have no data. The metadata will match exactly what it is supposed to be. The areas that this problem has occured for us are areas near the antimeridian (with standard WGS84) and at higher latitudes (where the non-linear orange-peel effect of Goode gets more pronounced). – scabecks Feb 27 at 8:00
  • @scabecks I have warped hundreds of rasters into the Homolosine and never got missing blocks. The error I report above concerns block input, not output, therefore I think your hypothesis can be ruled out. – Luís de Sousa Feb 28 at 14:33

Your Answer

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

Browse other questions tagged or ask your own question.