grass
grass copied to clipboard
OSGeo
"ERROR: Unable to initialize coordinate transformation" when reprojecting from one CRS to another CRS in a different mapset
I am trying to reproject a global raster (this: https://sedac.ciesin.columbia.edu/data/set/wildareas-v3-2009-human-footprint/data-download, I named it HFI) from Mollweide (its projection; EPSG:54009) to LAEA3035 (EPSG:3035) and clip it to Europe. Usually, these are the steps I follow:
-
I create a location with the same coordinate system of the file (in this case, I selected "choose the coordinate system from a referenced file", and used this raster as reference),
-
I switch to the Location LAEA3035, where I set the region and resolution:
g.region -p res=1000 raster=myEuropeRaster
projection: 99 (ETRS89_LAEA_Europe)
zone: 0
datum: etrs89
ellipsoid: grs80
north: 6824829.90382694
south: 1385914.39679116
west: 944055.47446584
east: 7602904.43906583
nsres: 999.98446535
ewres: 999.97731861
rows: 5439
cols: 6659
cells: 36218301
-
create a mask
r.mask raster=myEuropeRaster -
run
r.proj location=MOLL mapset=PERMANENT input=HFI output=HFIto load and reproject HFI raster from Mollweide location to LAEA3035 location, clipping it to Europe (thanks to the existing mask).
It isn't working. First, it throws two warnings:
WARNING: proj_trans() failed: acos/asin: |arg| >1.+1e-14
WARNING: proj_create() failed for '(null)'
and then an error:
ERROR: Unable to initialize coordinate transformation
This is the selected PROJ pipeline that appears:
Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000
+y_0=3210000 +ellps=GRS80 +step +proj=moll +lon_0=0 +x_0=0 +y_0=0
+ellps=WGS84
I am having the same problem with other rasters when I try to reproject them in the same way, and the problem persists if I use a different raster (always in Mollweide projection) to set up the Mollweide location. I searched online but couldn't find much information on how to solve this issue... I am using GRASS 7.8.6 on a server running Ubuntu 18.04.6 from a Windows 10 machine.
AFAIU the target CRS is not defined for global areas (see https://epsg.io/3035). Have you tried only cutting the piece you are interested in?
Hi veroandreo. If you meant that the target CRS is suitable only for European-level maps, I was aware of it. I usually load global maps into their appropriate location (I have had only WGS84 maps by far, this is the first time I deal with a different CRS) and reproject them to EPSG 3035, using a shapefile of Europe as a mask (to avoid the import of, for instance, parts of Greenland or Africa). In Mollweide case it is not working. After your comment I also tried to:
- reproject
myEuropeRasterto Mollweide location, - use it as a mask,
- load the global raster HFI (and it was clipped to
myEuropeRaster), - save this new cropped raster,
- switch to EPSG 3035 location,
- load and reproject the new cropped raster from Mollweide location.
but it throws the same error.
Are locations created with different versions of GRASS and PROJ? Could you list the content of the PERMANENT mapsets of both locations and post it here?
Hi ninsbl, the version of GRASS is always the same as above (7.8.6), and I don't know how to check PROJ version. I just switch from one PERMANENT mapset of a Location to another inside the same GRASS session, so I guess it is the same version. In both the PERMANENT mapsets there are only the two rasters (HFI for Mollweide EPSG:54009 Location, myEuropeRaster for LAEA3035 EPSG:3035 Location). However, I resolved using gdalwarp: I reprojected from Mollweide to WGS84, and then cropped to Europe while reprojecting into LAEA3035. But still couldn't understand the problem outlined in my question.
I don't know how to check PROJ version
You can run g.version -rge in a GRASS GIS terminal.
You cannot check back in time of course, but we may guess from the files in PERMANENT in the locations. If you did not update your libraries in between I may be on a wrong path...
Locations created with older versions of PROJ contain files like "PROJ_INFO" and "PROJ_UNITS", while locations created with newer version of PROJ also contain "PROJ_WKT" and probably "PROJ_SIRD".
My hunch is, that one of your locations lacks the PROJ_WKT file and that you are using a newer version of PROJ...
If you check the content of the directory on the file system that may shed light on the issue (or exclude possible problems).
BTW: the dataset is projected nicely in GRASS GIS 8 and 7.8.7 (OSGeo4W)... And I canot reproduce your issue either with Ubuntu 18.04 and GRASS 7.8.6 either...
r.proj --overwrite --verbose location=mollweide mapset=PERMANENT input=wildareas-v3-2009-human-footprint resolution=1000
Input CRS definition converted from '+proj=moll +lon_0=0 +x_0=0 +y_0=0
+no_defs +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +type=crs
' to '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
+type=crs'
Output CRS definition converted from '+proj=laea +lat_0=52 +lon_0=10
+x_0=4321000 +y_0=3210000 +no_defs +a=6378137 +rf=298.257222101
+towgs84=0.000,0.000,0.000 +type=crs ' to '+proj=laea +lat_0=52 +lon_0=10
+x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs'
Input Projection Parameters: +proj=moll +lon_0=0 +x_0=0 +y_0=0 +no_defs +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +type=crs
Input Unit Factor: 1
Output Projection Parameters: +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +no_defs +a=6378137 +rf=298.257222101 +towgs84=0.000,0.000,0.000 +type=crs
Output Unit Factor: 1
Input:
Cols: 3150 (original: 36081)
Rows: 1550 (original: 16382)
North: 7674957.052787 (original: 9018957.052787)
South: 6124957.052787 (original: -7363042.947213)
West: -330094.097520 (original: -18040094.097520)
East: 2819905.902480 (original: 18040905.902480)
EW-res: 1000.000000
NS-res: 1000.000000
Output:
Cols: 2394 (original: 2394368)
Rows: 1500 (original: 1500471)
North: 5074164.000000 (original: 5074164.000000)
South: 3573693.000000 (original: 3573693.000000)
West: 3650850.000000 (original: 3650850.000000)
East: 6045218.000000 (original: 6045218.000000)
EW-res: 1000.153718
NS-res: 1000.314000
100.00 percent are kept in memory
Allocating memory and reading input raster map...
100%
Projecting...
100%
r.proj complete.
It seems specific to your setup / system...
Could you also provide the output of:
g.proj -w
from both locations?
Hi ninsbl, I checked on my file browser and both Mollweide and LAEA3035 have the PROJ_INFO, PROJ_UNITS, and PROJ_WKT files in their PERMANENT mapsets.
This is the output of g.proj -w from LAEA3035:
PROJCRS["ETRS89_LAEA_Europe",
BASEGEOGCRS["ETRS89",
DATUM["European Terrestrial Reference System 1989",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ID["EPSG",6258]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",52,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",10,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",4321000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",3210000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
And from Mollweide:
PROJCRS["World_Mollweide",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Mollweide"],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
PROJ version is 7.0.0.
Thanks. Here is a reproducible example run on Ubuntu 18.04 with packages from ubuntugis-unstable:
mkdir ./tmp
grass -c ./wildareas-v3-2009-human-footprint.tif -e ./tmp/mol
grass -c EPSG:3035 -e ./tmp/laea
grass ./tmp/mol/PERMANENT --exec r.external --verbose input=/data/R/Kladd/sbl/wildareas-v3-2009-human-footprint.tif output=HFI
grass ./tmp/laea/PERMANENT --exec g.region -p n=6000000 e=6050000 s=3000000 w=3000000 res=1000
grass ./tmp/laea/PERMANENT --exec r.proj location=mol mapset=PERMANENT input=HFI
r.proj runs fine:
Input CRS definition converted from '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +no_defs +a=6378137 +rf=298.257223563 +towgs84=0.000,0.000,0.000 +type=crs ' to '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs'
Output CRS definition converted from '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +no_defs +a=6378137 +rf=298.257222101 +towgs84=0.000,0.000,0.000 +type=crs ' to '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs +type=crs'
Input:
Cols: 4636 (original: 36081)
Rows: 2720 (original: 16382)
North: 8339957.052787 (original: 9018957.052787)
South: 5619957.052787 (original: -7363042.947213)
West: -1540094.097520 (original: -18040094.097520)
East: 3095905.902480 (original: 18040905.902480)
EW-res: 1000.000000
NS-res: 1000.000000
Output:
Cols: 3050 (original: 3050)
Rows: 3000 (original: 3000)
North: 6000000.000000 (original: 6000000.000000)
South: 3000000.000000 (original: 3000000.000000)
West: 3000000.000000 (original: 3000000.000000)
East: 6050000.000000 (original: 6050000.000000)
EW-res: 1000.000000
NS-res: 1000.000000
Allocating memory and reading input raster map...
Projecting...
r.proj complete.
However, note that my GRASS version is 7.8.5 (not 7.8.6).
grass ./tmp/laea/PERMANENT --exec g.version -reg
version=7.8.5
date=2020
revision=exported
build_date=2020-12-22
build_platform=x86_64-pc-linux-gnu
build_off_t_size=8
libgis_revision=2020-12-22T15:55:09+00:00
libgis_date=2020-12-22T15:55:09+00:00
proj=7.0.0
gdal=3.0.4
geos=3.8.0
sqlite=3.22.0
And my CRS definition for mollweide is different:
grass ./tmp/mol/PERMANENT --exec g.proj -w
PROJCS["unknown",
GEOGCS["wgs84",
DATUM["WGS_1984",
SPHEROID["WGS_1984",6378137,298.257223563]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]]],
PROJECTION["Mollweide"],
PARAMETER["central_meridian",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]
There are several PROJ related changes in 7.8.6: https://trac.osgeo.org/grass/wiki/Release/7.8.6-News that may cause the difference here.