search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.75k stars 2.5k forks source link

Reprojection to Goode's Homolosine duplicates Iceland and eastern Greenland #959

Closed ldesousa closed 2 years ago

ldesousa commented 5 years ago

Expected behavior and actual behavior.

I am reprojecting a global raster that is defined on constant longitude-latitude intervals into Goode's projection (no datum shift). Things go more less well, but Iceland and eastern Greeland appear twice, and ouside the counter-domain of Goode's projection:

goodegdal

Steps to reproduce the problem.

This is what gdal_info says about the input raster:

$ gdalinfo GAUL_ADMIN0_landmask_250m.tif
Driver: GTiff/GeoTIFF
Files: GAUL_ADMIN0_landmask_250m.tif
       GAUL_ADMIN0_landmask_250m.tif.aux.xml
Size is 172800, 71698
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"]]
Origin = (-180.000000000000000,87.370000000000005)
Pixel Size = (0.002083333000000,-0.002083333000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0000000,  87.3700000) (180d 0' 0.00"W, 87d22'12.00"N)
Lower Left  (-180.0000000, -62.0008094) (180d 0' 0.00"W, 62d 0' 2.91"S)
Upper Right ( 179.9999424,  87.3700000) (179d59'59.79"E, 87d22'12.00"N)
Lower Right ( 179.9999424, -62.0008094) (179d59'59.79"E, 62d 0' 2.91"S)
Center      (  -0.0000288,  12.6845953) (  0d 0' 0.10"W, 12d41' 4.54"N)
Band 1 Block=172800x1 Type=Int16, ColorInterp=Gray
  Min=2.000 Max=276.000 
  Minimum=2.000, Maximum=276.000, Mean=138.243, StdDev=85.746
  NoData Value=-32768
  Metadata:
    STATISTICS_MAXIMUM=276
    STATISTICS_MEAN=138.24278392441
    STATISTICS_MINIMUM=2
    STATISTICS_STDDEV=85.745958703755

This is the gdalwarp command:

gdalwarp -overwrite -s_srs EPSG:4326 -t_srs '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' -co "COMPRESS=DEFLATE" -te -20037500 -6729000 20037250 8600750 -tr 250 250 -wt Int16 -ot Int16 GAUL_ADMIN0_landmask_250m.tif GAUL_ADMIN0_landmask_Goode.tif

And gdalinfoon the output:

$ gdalinfo GAUL_ADMIN0_landmask_Goode.tif
Driver: GTiff/GeoTIFF
Files: GAUL_ADMIN0_landmask_Goode.tif
       GAUL_ADMIN0_landmask_Goode.tif.aux.xml
Size is 160299, 34403
Coordinate System is:
PROJCS["Interrupted_Goode_Homolosine",
    GEOGCS["GCS_WGS_1984",
        DATUM["D_WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Interrupted_Goode_Homolosine"],
    UNIT["Meter",1]]
Origin = (-20037500.000000000000000,8600750.000000000000000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-20037500.000, 8600750.000) 
Lower Left  (-20037500.000,       0.000) (179d59'59.73"W,  0d 0' 0.01"N)
Upper Right (20037250.000, 8600750.000) 
Lower Right (20037250.000,       0.000) (179d59'51.65"E,  0d 0' 0.01"N)
Center      (    -125.000, 4300375.000) (  8d24'16.92"W, 38d37'51.33"N)
Band 1 Block=160299x1 Type=Int16, ColorInterp=Gray
  Min=1.000 Max=274.000 
  Minimum=1.000, Maximum=274.000, Mean=147.335, StdDev=80.657
  NoData Value=-32768
  Metadata:
    STATISTICS_MAXIMUM=274
    STATISTICS_MEAN=147.33473144028
    STATISTICS_MINIMUM=1
    STATISTICS_STDDEV=80.656559922178

Operating system

Ubuntu 16.04 64 bit Ubuntu 18.04 64 bit

GDAL version and provenance

2.2.2 version from ubuntugis-unstable PPA (Ubuntu 16.04) 2.2.3 version from ubuntugis-unstable PPA (Ubuntu 18.04)

Extra

The end result is the same with either version 2.2.2 and 2.2.3. However, with 2.2.2 gdalwarp takes about 5 minutes to complete and produces a 50 MB tif, whereas 2.2.3 takes 75 min and produces a 700 MB tif.

If you wish to obtain the input file, please contact me by e-mail or through gitter.im.

ldesousa commented 5 years ago

I created a vector file with the Homolosine counter-domain to verify the resulting raster against it. There are more replications than I initially reported. Screen captures from GRASS: goodeglobe goodewest goodeeast

This raster does not have Antarctica, but it is possible similar issues happen there.

rouault commented 5 years ago

I can't replicate with the small_world.tif from https://github.com/OSGeo/gdal/blob/master/autotest/gdrivers/data/small_world.tif

wget https://github.com/OSGeo/gdal/blob/master/autotest/gdrivers/data/small_world.tif?raw=true -O in.tif
gdalwarp in.tif out.tif -overwrite -t_srs "+proj=igh"

out

ldesousa commented 5 years ago

Here is an example with a public dataset. Ran on Ubuntu 18.04 with GDAL 2.2.3:

wget https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_sovereignty.zip
unzip ne_10m_admin_0_sovereignty.zip
gdal_rasterize -a NE_ID -a_nodata 0 -tr 0.01 0.01 -of GTiff ne_10m_admin_0_sovereignty.shp ne_10m_admin_0_sovereignty.tif
gdalwarp -overwrite -s_srs EPSG:4326 -t_srs '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs' -co "COMPRESS=DEFLATE" -te -20037500 -6729000 20037250 8600750 -tr 250 250 ne_10m_admin_0_sovereignty.tif ne_10m_admin_0_sovereignty_Goode.tif
rouault commented 5 years ago

OK so the problem is that the invert IGH projection doesn't check that the input (x,y) is in the domain of output of the forward path. Consequently points like -6367850 7051386 and 629308.892412171 7051386 map to the same point in (long,lat): -17.139173815748 65.529497121681

The solution is to run gdalwarp with "--config CHECK_WITH_INVERT_PROJ TRUE", so that after the inverse projection, a forward one is re-done.

The fact that the issue appears at some particular combination of input/output size/resolution comes from the chunking done by gdalwarp. Some chunks that are mostly outside the validity domain are rejected, wheras others that have mostly inside are still processed, but with those potential out of domain points.

ldesousa commented 5 years ago

I confirm that the parameter --config CHECK_WITH_INVERT_PROJ TRUE produces correct results. However, transformations made with programmes like QGis or GRASS produce the errors reported initially, meaning that the software is not conceived for cases like this.

Should perhaps the CHECK_WITH_INVERT_PROJ parameter be true by default?

thengl commented 5 years ago

I have made a function to reproject any large global layer to Goode Homolosine in parallel (see latlon2gh function. This seems to produce OK result (see: https://doi.org/10.5281/zenodo.3355006) although there are still some smaller pixels on the upper edges of the map that get lost in reprojection. I've used GDAL 2.3.1 for processing. Thanks @ldesousa for preparing the domain for GH and for sharing!

rouault commented 2 years ago

it should probably the job of PROJ's inverse Interrupted Goode Homolosine implementation to reject points that are outside the validity domain. If the issue still matters, please raise the issue on PROJ side