4

I am experiencing a strange problem by using the GRASS toolbox on QGIS 2.14.3 Essen.

Even if in settings>options>CRS I select to use the project CRS when creating new layers, the results of GRASS toolbox tools, like r.watershed and r.fill.dir, have an "user-defined CRS" with code 100000, even if the project CRS is EPSG 20137.

How can I solve this?

1

It is the result of the differences between the Proj4 parameters in QGIS and GRASS GIS:

1) If you look at the definition of the EPSG:20137 projection, the Proj4 string is:

+proj=utm +zone=37 +ellps=clrk80 +units=m +no_defs 

2) The GRASS7 commands of the Processing Toolbox use

a) the script /../python/plugins/processing/algs/grass7/Grass7Algorithm.py to obtain the projection of the layer. In the Python console

proj4 = unicode(layer.crs().toProj4())
print proj4 
+proj=utm +zone=37 +ellps=clrk80 +towgs84=-166,-15,204,0,0,0,0 +units=m +no_defs

b) Then the script /.../python/plugins/processing/algs/grass7/Grass7Utils.py creates the /.../.qgis2/processing/grass7_batch_job.sh (or bat) file which is processed by the Python module subprocess

Beginning of the grass7_batch_job.sh file

g.proj -c proj4="+proj=utm +zone=37 +ellps=clrk80 +towgs84=-166,-15,204,0,0,0,0 +units=m +no_defs"
v.in.ogr min_area=0.0001 snap=-1 input="/Users/Shared" layer=line 20137 output=tmp1468150743135 --overwrite -o
g.region -a n=6164723.90083 s=6163210.42107 e=-1860027.15695 w=-1862095.02489 res=100

3) If I use EPSG:20137 to create a LOCATION in GRASS GIS

GRASS 7.1 (test20137):~ > g.proj -g
name=Universal Transverse Mercator
proj=utm
ellps=clark80
zone=37
no_defs=defined
epsg=20137
unit=meter
units=meters
meters=1

And

GRASS 7.1 (test20137):~ > g.proj -j
+proj=utm
+no_defs
+zone=37
+a=6378249.145
+rf=293.465
+to_meter=1

4) If I use the QGIS Proj4 string to create a new GRASS GIS LOCATION

GRASS 7.1 (test_QGIS):~ > g.proj -g
name=Universal Transverse Mercator
proj=utm
a=6378249.145
es=0.006803481196021999
zone=37
towgs84=-166,-15,204,0,0,0,0
no_defs=defined
unit=Meter
units=Meters
meters=1

And

GRASS 7.1 (test_QGIS):~ > g.proj -j
+proj=utm
+no_defs
+zone=37
+a=6378249.145
+rf=293.4662999999974
+towgs84=-166,-15,204,0,0,0,0
+to_meter=1

5) Therefore for GRASS GIS the projections are not formally the same -> your "user-defined CRS" in QGIS.

.prj file of the original shapefile (EPSG:20137)

PROJCS["Adindan_UTM_zone_37N",GEOGCS["GCS_Adindan",DATUM["D_Adindan",SPHEROID["Clarke_1880_RGS",6378249.145,293.465]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

.prj file of the resulting shapefile (EPSG:????)

PROJCS["UTM_Zone_37_Northern_Hemisphere",GEOGCS["GCS_Unknown",DATUM["D_unknown",SPHEROID["Unknown",6378249.145,293.4662999999974]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",39],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

6) The solution ?

You can try to modify the Proj4 string in the file grass7_batch_job.sh or in the script Grass7Algorithm.py but is easier to simply modify the properties of the layer in QGIS.

  • Just one more question... But even if the CRS are not formally the same, are they practically the same? In this case I guess the best solution is to change the properties .If the CRS is not practically the same, the program will read the X coordinate of USER defined CRS as the 20137 one, leading to a mistake... only if the coordinates are the same, changing the CRS in the properties will not affect anything. Let me know if I have understood. Thanks – giulio Jul 10 '16 at 18:47
  • Yes, you have understood. – gene Jul 10 '16 at 19:22
  • The proj. definition at spatialreference.org is wrong, it misses the +towgs84 parameter. QGIS 2.14.3 and GRASS 7.04 use the correct values. I wonder how you get the GRASS GIS Location without +towgs84 parameters. – AndreJ Jul 11 '16 at 6:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.