[Bug] r.proj from EPSG:3031 to EPSG:4326 returns one grid cell

[mankoff]

Describe the bug

I'm trying to reproject an Antarctic dataset on EPSG:3031 to EPSG:4326. I can successfully do Greenland EPSG:3413 to EPSG:4326, but the same commands for Antarctica return 1 grid cell and report

To Reproduce

1. Run the following commands

Set up 4326 mapset

grass --text -c EPSG:4326 ./G_4326
g.region res=1 -pa w=-180 e=180 s=-90 n=90 -pas # 1° x 1°
exit

3031 mapset

grass -c EPSG:3031 ./G_3031
# 500 m resolution, matching BedMachine boundaries
bound=3333250
g.region res=500 w=-${bound} e=${bound} s=-${bound} n=${bound} -ps
r.mapcalc "foo = 42"
exit

Reproject

grass ./G_4326/PERMANENT
r.proj location=./G_3031 mapset=PERMANENT input=foo # output: 1 column!

Expected behavior

The input dataset spans all longitudes and covers Antarctica. I'd expect the output dataset to do this, reprojected into EPSG:4326.

Screenshots

Screenshots come from earlier (see first version of this issue under the edit item above) using BedMachine. The current g.region matches that from BedMachine.

3031:

4326:

System description (please complete the following information):

t480:~ $ lsb_release -a
No LSB modules are available.
Distributor ID: Ubuntu
Description:    Ubuntu 22.04.4 LTS
Release:    22.04
Codename:   jammy
t480:~ $ uname -a
Linux t480 5.15.0-97-generic #107-Ubuntu SMP Wed Feb 7 13:26:48 UTC 2024 x86_64 x86_64 x86_64 GNU/Linux
t480:~ $ ```

version=8.3.2 date=2024 revision=exported build_date=2024-03-08 build_platform=x86_64-pc-linux-gnu build_off_t_size=8 libgis_revision=8.3.2 libgis_date=2024-03-08T09:26:59+00:00 proj=9.3.1 gdal=3.8.4 geos=3.12.1 sqlite=3.37.2

[dhdeangelis]

I tried but could not reproduce this. See below.

I did not use the dataset you mention because I did not manage to create a user login at Earthdata that is necessary for downloading. So instead I used old data I already had, using the same projections. It is of course not an exact attempt.

Also, my locations were created many years ago and perhaps if we are dealing with problems in how projection systems are defined and managed that could be a problem.

Can you reproject the grid using just GDAL (gdalwarp)?

---

System Info                                                                     
GRASS version: 8.4.0dev                                                         
Code revision: exported                                                         
Build date: 2024-04-26                                                          
Build platform: x86_64-pc-linux-gnu                                             
GDAL: 3.8.5                                                                     
PROJ: 9.4.0                                                                     
GEOS: 3.12.1                                                                    
SQLite: 3.45.2                                                                  
Python: 3.11.9                                                                  
wxPython: 4.2.1                                                                 
Platform: Linux-6.8.8-1-default-x86_64-with-glibc2.39                           

Input data: Polar Stereographic, 1 km cell size

> g.region -p
projection: 99 (Stereographic)
zone:       0
datum:      wgs84
ellipsoid:  a=6378137 es=0.00669437999014138
north:      2384000
south:      -2397000
west:       -2836000
east:       2814000
nsres:      1000
ewres:      1000
rows:       4781
cols:       5650
cells:      27012650

> g.proj -j
+proj=stere
+lat_0=-90
+lat_ts=-71
+lon_0=0
+k=1
+x_0=0
+y_0=0
+no_defs
+a=6378137
+rf=298.257223563
+towgs84=0.000,0.000,0.000
+type=crs 
+to_meter=1

Output project: LatLong 4326

> g.region -p
projection: 3 (Latitude-Longitude)
zone:       0
datum:      wgs84
ellipsoid:  wgs84
north:      60S
south:      90S
west:       180W
east:       180E
nsres:      0:15
ewres:      0:15
rows:       120
cols:       1440
cells:      172800

 g.proj -j
+proj=longlat
+a=6378137
+rf=298.257223563
+no_defs
+towgs84=0.000,0.000,0.000
+type=crs

Reprojection:

> r.proj project=antartida_1km mapset=PERMANENT input=accumulation output=test
Selected PROJ pipeline:
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step
+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +no_defs
+over +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000
************************
Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1
+x_0=0 +y_0=0 +no_defs +over +a=6378137 +rf=298.257223563
+towgs84=0.000,0.000,0.000 +step +proj=unitconvert +xy_in=rad +xy_out=deg
************************

Input:
Cols: 5736 (original: 5736)
Rows: 4916 (original: 4916)
North: 2458000.000000 (original: 2458000.000000)
South: -2458000.000000 (original: -2458000.000000)
West: -2868000.000000 (original: -2868000.000000)
East: 2868000.000000 (original: 2868000.000000)
EW-res: 1000.000000
NS-res: 1000.000000

Output:
Cols: 1440 (original: 1440)
Rows: 120 (original: 120)
North: -60.000000 (original: -60.000000)
South: -90.000000 (original: -90.000000)
West: -180.000000 (original: -180.000000)
East: 180.000000 (original: 180.000000)
EW-res: 0.250000
NS-res: 0.250000

Allocating memory and reading input raster map...
 100%
Selected PROJ pipeline:
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step
+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +no_defs
+over +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000
************************
Projecting...
 100%
r.proj complete.

[veroandreo]

I downloaded the dataset, followed the exact same steps described, and I can reproduce the behavior you observe @mankoff.

System info:

GRASS version: 8.4.0dev                                                         
Code revision: a351e5eb69                                                       
Build date: 2024-05-02                                                          
Build platform: x86_64-pc-linux-gnu                                             
GDAL: 3.7.1                                                                     
PROJ: 9.2.1                                                                     
GEOS: 3.12.0                                                                    
SQLite: 3.42.0                                                                  
Python: 3.11.6                                                                  
wxPython: 4.2.1                                                                 
Platform: Linux-6.5.0-28-generic-x86_64-with-glibc2.38                          

[mankoff]

@dhdeangelis

Can you reproject the grid using just GDAL (gdalwarp)?

Yes. With gdalwarp -t_srs EPSG:4326 NETCDF:"BedMachineAntarctica-v3.nc":mask output.tiff everything appears OK to me. I'm also happy to share the file if you want to help debug this, but don't want to post a public link to DropBox and am not sure how else to share it.

[mankoff]

Same issue when using GRASS 7.8.8 from mundialis/grass-py3-pdal:7.8.8-alpine https://grass.osgeo.org/download/docker/#GRASS-GIS-old

I've reproduced this without BedMachine. I apologize for the non-minimal example before. I include an MWE below and I will edit the original post with this simplified code.

Set up 4326 mapset

grass --text -c EPSG:4326 ./G_4326
g.region res=1 -pa w=-180 e=180 s=-90 n=90 -pas # 1° x 1°
exit

3031 mapset

grass -c EPSG:3031 ./G_3031
# 500 m resolution, matching BedMachine boundaries
bound=3333250
g.region res=500 w=-${bound} e=${bound} s=-${bound} n=${bound} -ps
r.mapcalc "foo = 42"
exit

Reproject

grass ./G_4326/PERMANENT
r.proj location=./G_3031 mapset=PERMANENT input=foo # output: 1 column!

[mankoff]

I note this is sensitive. If I change g.region res=500 to res=250 or res=1000 the bug does not appear.

[dhdeangelis]

I have downloaded the BedMachine dataset and could reproduce the error.

EDIT: I can also confirm that projecting from a 1 km resolution region to a LonLat 1/4° resolution region works well.

[dhdeangelis]

I may have found a solution to this. It involves using the -n option:

From the r.proj manual (https://grass.osgeo.org/grass84/manuals/r.proj.html):

When reprojecting whole-world maps the user should disable map-trimming with the -n flag. Trimming is not useful here because the module has the whole map in memory anyway. Besides that, world "edges" are hard (or impossible) to find in CRSs other than latitude-longitude so results may be odd with trimming.

I tried and worked for me using the 500 m region as source and a LonLat region as target:

> r.proj -n project=antartida_500m mapset=PERMANENT input=bm

Input:
Cols: 13333 (original: 13333)
Rows: 13333 (original: 13333)
North: 3333250.000000 (original: 3333250.000000)
South: -3333250.000000 (original: -3333250.000000)
West: -3333250.000000 (original: -3333250.000000)
East: 3333250.000000 (original: 3333250.000000)
EW-res: 500.000000
NS-res: 500.000000

Output:
Cols: 1440 (original: 1440)
Rows: 133 (original: 133)
North: -56.750000 (original: -56.750000)
South: -90.000000 (original: -90.000000)
West: -180.000000 (original: -180.000000)
East: 180.000000 (original: 180.000000)
EW-res: 0.250000
NS-res: 0.250000

Allocating memory and reading input raster map...
 100%
Selected PROJ pipeline:
+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step
+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84
************************
Projecting...
 100%
r.proj complete.

[dhdeangelis]

@mankoff @veroandreo is this still of interest? IMO the problem was that the reprojection involving the pole and -n seems to take care of that. Perhaps some clarification can be added to the manual and then this could be closed?

[mankoff]

@dhdeangelis thanks for keeping this in mind and helping me solve the issue. I've added a small PR to recommend looking more carefully at the Notes section of the document if working in a global lat-lon setup. As I say in the PR, I'm not sure this is needed - maybe I just need to read the docs more carefully.