search

OSGeo / grass

GRASS GIS - free and open-source geospatial processing engine
https://grass.osgeo.org
Other
849 stars 308 forks source link

[Bug] r.proj fails with global dataset exceeding 90N #2278

Open isaacullah opened 2 years ago

isaacullah commented 2 years ago

Describe the bug It seems I am experiencing a similar issue to #2757 when trying to reproject a LongLat map with extents greater than 90N/S or 180E/W. The particular dataset is from CHELSA and is based on GTED2010. I can import into a WGS84 LongLat location (EPSG 4326) using r.in.gdal with no issue, but I cannot then reproject using r.proj into another CRS (in my case WGS84 UTM Z36N). I have tried using flags -l and -a in r.in.gdal, but that does not affect the r.proj error.

To Reproduce Steps to reproduce the behavior:

1) Download some global dataset from CHELSA https://chelsa-climate.org/downloads/

2) Import into a WGS84 longlat location (EPSG 4326) using r.import or r.in.gdal. Optionally, try flags -l and/or -a in r.in.gdal. Try also to set an initial region to crop map on import. All import attempts work fine to get the map into this longlat location.

3) Change to projected CRS location (in my case WGS84 UTMz36n) and attempt to reproject using r.proj: r.proj --overwrite location=world_latlong_wgs84 mapset=PERMANENT input=CHELSA_TraCE21k_AP_2kya dbase=/home/iullah/grassdata memory=1000 resolution=1000 Get following error: proj_create: Error -1 (no arguments in initialization list): Pipeline: Bad step definition: init=epsg:4326 (no arguments in initialization list) WARNING: proj_create() failed for '+proj=pipeline +step +inv +proj=utm +zone=36 +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000 +step +init=epsg:4326' ERROR: Unable to initialize coordinate transformation

4) Put in first operation in proj pipeline using r.proj "pipeline" variable: r.proj --overwrite location=world_latlong_wgs84 mapset=PERMANENT input=CHELSA_TraCE21k_AP_2kya dbase=/home/iullah/grassdata memory=1000 resolution=1000 pipeline=+proj=longlat Now get new error: WARNING: proj_trans() failed: latitude or longitude exceeded limits ERROR: unable to transform coordinates 733465, 3.4442e+06

5) r.proj again, but checking flag -n r.proj -n --overwrite location=world_latlong_wgs84 mapset=PERMANENT input=CHELSA_TraCE21k_AP_2kya dbase=/home/iullah/grassdata memory=1000 resolution=1000 pipeline=+proj=longlat Get following error:

Input:
Cols: 81 (original: 81)
Rows: 45 (original: 45)
North: 31.116528 (original: 31.116528)
South: 30.741528 (original: 30.741528)
West: 35.433333 (original: 35.433333)
East: 36.108333 (original: 36.108333)
EW-res: 0.008333
NS-res: 0.008333

Output:
Cols: 63 (original: 2092)
Rows: 39 (original: 1291)
North: 3444220.000000 (original: 3444220.000000)
South: 3405490.000000 (original: 3405490.000000)
West: 733450.000000 (original: 733450.000000)
East: 796210.000000 (original: 796210.000000)
EW-res: 996.190476
NS-res: 993.076923

Allocating memory and reading input raster map...
Projecting...
WARNING: proj_trans() failed: latitude or longitude exceeded limits
ERROR: Error in GPJ_transform()

Note that I have tried this with the entire global map, and it is clear that the extent boundaries exceed 90N. However, even if I "pre-cropped" this map upon initial import by setting a smaller region window, this error continues (as above).

Expected behavior I should be able to reproject this data from longlat to a projected CRS with no issues.

Screenshots If applicable, add screenshots to help explain your problem.

System description (please complete the following information):

Additional context I have had a similar problem with other global datasets, including WorldClim2.0

metzm commented 2 years ago

What is the version of PROJ you are using?

For the projected CRS location (in your case WGS84 UTMz36n), can you provide the output of g.proj -p, g.proj -j, g.proj -w?

Generally, it is mathematically not possible to reproject a global map from EPSG:4326 to a UTM zone.

petrasovaa commented 2 years ago

There is a similar problem in #1658.

Generally, it is mathematically not possible to reproject a global map from EPSG:4326 to a UTM zone.

Perhaps, but it is not obvious to the user, they might expect distorted results outside of the zone, but not complete failure. Any idea how does gdal deal with that?

metzm commented 2 years ago

There is a similar problem in #1658.

Generally, it is mathematically not possible to reproject a global map from EPSG:4326 to a UTM zone.

Perhaps, but it is not obvious to the user, they might expect distorted results outside of the zone, but not complete failure. Any idea how does gdal deal with that?

The error comes from proj, not from gdal. I do expect complete failure in this case. By definition, it is simply not possible to reproject a global map to a UTM zone, or many other projected CRS's.

isaacullah commented 2 years ago

To expand, the expectation was not to reproject the entire global map to one UTM zone, but to reproject a portion of it (i.e., that which falls within my computational region) so that I can do math with the portion of these maps that fall within my UTM zone only. I don't have access to the computer on which I was doing this the other day (it's Spring Break now, and I am at home and not on campus), BUT, I did find a kind of work around late yesterday. I was able to import these global maps into the "Natural Earth" sample location, and checking the "override projection information (use location's projection)" check box. From, THIS location, I could then reproject to the UTM, but doing the same thing from the default "world_latlong_wgs84" location did NOT reproject. I know I should probably have just made a new location from scratch to try this out the first time, but I am doing this in the context of a class, where I am trying to teach novice students to import and reproject data. I thought, since the "world_latlong_wgs84" was already ostensibly in right projection (4326), that it would make sense for students to just use that instead of setting up a new blank lat-long location (and my understanding is that this new default first location is SUPPOSED to make it easier for newbies to get going). What I don't understand is why it could work the the natural_earth location and NOT the world_latlong_wgs84 location when they are both EPSG 4326?

petrasovaa commented 2 years ago

g.proj -w shows some differences, not sure how much they matter:

natural_earth_dataset (downloaded using GUI):

g.proj -w                                                                       
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["wgs84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS_1984",6378137,298.257223563,
                    LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8901]],
            CS[ellipsoidal,2],
                AXIS["longitude",east,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]],
                AXIS["latitude",north,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from wgs84 to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]

world_latlong_wgs84:

GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
metzm commented 2 years ago

A BOUNDCRS should not happen, this restricts the CRS in this case to a specific datum transformation.

What is the version of PROJ you are using?

How did you download the natural_earth_dataset using GUI, apparently it created a new location?

petrasovaa commented 2 years ago

I downloaded the natural earth dataset available through data catalog GUI, I assume it's the basic one here: https://grass.osgeo.org/download/data/

But @isaacullah had problem with the other one if I understood correctly...

isaacullah commented 2 years ago

Yes, I too downloaded the natural earth dataset through the GUI. THAT is the latlong location that I was able to "force" fit the CHELSA data into with a projection override, and then reprojection worked. The "standard" world latlong sample location that is now loaded by default on first launch of GRASS 8 did NOT work with projection override.

On Fri, Mar 25, 2022 at 1:42 PM Anna Petrasova @.***> wrote:

I downloaded the natural earth dataset available through data catalog GUI, I assume it's the basic one here: https://grass.osgeo.org/download/data/

But @isaacullah https://github.com/isaacullah had problem with the other one if I understood correctly...

— Reply to this email directly, view it on GitHub https://github.com/OSGeo/grass/issues/2278#issuecomment-1079419543, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACDYS5DUNT2S4LJDJBITAIDVBYQL5ANCNFSM5RR3IX6Q . You are receiving this because you were mentioned.Message ID: @.***>

--

Isaac I. Ullah, PhD Recent papers: Laborscapes and Archaeologies of Sustainability: http://dx.doi.org/10.1558/jma.39327 Water, dust, and loess: co-evolution of landscapes, farming, and human society: https://doi.org/10.1016/j.jaa.2019.101067 Modeling feedbacks between human and natural processes in the land system: https://doi.org/10.5194/esd-9-895-2018 https://doi.org/10.5194/esd-9-895-2018 https://doi.org/10.5194/esd-9-895-2018

metzm commented 2 years ago

Back to the initial error: proj_create: Error -1 (no arguments in initialization list): Pipeline: Bad step definition: init=epsg:4326 (no arguments in initialization list) WARNING: proj_create() failed for '+proj=pipeline +step +inv +proj=utm +zone=36 +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000 +step +init=epsg:4326' ERROR: Unable to initialize coordinate transformation

What is the version of PROJ you are using?

What is the output of g.proj -p, g.proj -j, g.proj -w for both the source and target GRASS location?