Comments (9)
@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.
from grass.
@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?
from grass.
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.
from grass.
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
from grass.
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.
from grass.
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!
from grass.
I note this is sensitive. If I change g.region res=500
to res=250
or res=1000
the bug does not appear.
from grass.
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.
from grass.
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.
from grass.
Related Issues (20)
- [Bug] PATH setting discrepancy between the main exectable and Python API for scripts on Windows
- [Bug] wxGUI: Nonsense commands end in Success in gconsole OnCmdDone method
- [Bug] Windows CI builds fail to compile r.flow HOT 6
- Renovate naming convention issue HOT 2
- Raster layers cannot be displayed[Bug] HOT 1
- [Bug] wxGUI/history: update cross to green check after updating to current region
- Tests are failing in CI with an ISO format error HOT 6
- r.recode results in integer when floating numbers are expected [Bug] HOT 6
- [Bug] Windows GUI can't launch for 8.4.0dev (nightly) HOT 7
- [Build] Exclude specific module (r.geomorphon) at compilation HOT 17
- [Feat] temporal: revert #3723 (fix for SQLite version < 3.33)
- [Bug] [Windows] Crash of the GRASS-GIS 8.3 `v.surf.rst.exe` (OSGeo4W) and Windows module `ntdll.dll` HOT 13
- [Bug] Windows 8.4.0dev cannot launch GUI with error in `gui\wxpython\lmgr\statusbar.py` HOT 12
- [Bug] Windows OSGeo4W builds seem to report outdated build info HOT 2
- import nc file in grass gis HOT 2
- [Feat] library: Add Standard parser options for date and time HOT 1
- [Feat] i.group: print image group content with semantic_labels in JSON HOT 3
- [Bug] db.univar fails with Python error when input table not existing
- [Feat] Rename location to project in images in docs HOT 2
- [Bug] Travis CI duplicated runs on renovate branches
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from grass.