[GRASS-dev] custom CRS

GRASS GRASS developer
,
1

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

The corresponding proj.4 definition (as defined in the esri epsg-codes file in /usr/share/proj/)

+proj=longlat +a=6377301.243 +b=6356100.230165384 +no_defs

And when I create a location based on this shapefile, the projection information is:

-PROJ_INFO-------------------------------------------------
name : GCS_Everest_India_Nepal
a : 6377301.243
es : 0.006637846068423442
proj : ll
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : degree
units : degrees
meters : 1.0

Question 1: Now, I can import a shapefile or raster layer in latlon/WGS84 (EPGS: 4326) without reprojection. GRASS thus assumes the projection is identical, even if they are not? Why is that?

The proj.4 definition above is missing the shift w.r.t. the WGS84. The correct proj.4 definition for EPGS 37203 would be [2]

+proj=longlat +a=6377301.243 +b=6356100.2284 +towgs84=295,736,257,0,0,0,0 +no_defs

Question 2: In QGIS it is possible to define a custom CRS using the above. Can I do this in GRASS as well, and if so, how? Or can / should I add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Regards,

Paulo

[1] https://epsg.io/37203
[2] https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689

2

A link that provides further relevant information

···

On 18-03-17 10:06, Paulo van Breugel wrote:

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in the prj file as

GEOGCS[“GCS_Everest_India_Nepal”,DATUM[“D_Everest_India_Nepal”,SPHEROID[“Everest_Definition_1962”,6377301.243,300.8017255]],PRIMEM[“Greenwich”,0.0],UNIT[“Degree”,0.0174532925199433]]

The corresponding proj.4 definition (as defined in the esri epsg-codes file in /usr/share/proj/)

+proj=longlat +a=6377301.243 +b=6356100.230165384 +no_defs

And when I create a location based on this shapefile, the projection information is:

-PROJ_INFO-------------------------------------------------
name : GCS_Everest_India_Nepal
a : 6377301.243
es : 0.006637846068423442
proj : ll
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : degree
units : degrees
meters : 1.0

Question 1: Now, I can import a shapefile or raster layer in latlon/WGS84 (EPGS: 4326) without reprojection. GRASS thus assumes the projection is identical, even if they are not? Why is that?

The proj.4 definition above is missing the shift w.r.t. the WGS84. The correct proj.4 definition for EPGS 37203 would be [2]

+proj=longlat +a=6377301.243 +b=6356100.2284 +towgs84=295,736,257,0,0,0,0 +no_defs

Question 2: In QGIS it is possible to define a custom CRS using the above. Can I do this in GRASS as well, and if so, how? Or can / should I add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Regards,

Paulo

[1] https://epsg.io/37203
[2] https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689

https://github.com/klokantech/epsg.io/issues/49

3

Le 18 mars 2017 10:06:09 GMT+01:00, Paulo van Breugel <p.vanbreugel@gmail.com> a écrit :

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in
the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

The corresponding proj.4 definition (as defined in the esri epsg-codes
file in /usr/share/proj/)

+proj=longlat +a=6377301.243 +b=6356100.230165384 +no_defs

And when I create a location based on this shapefile, the projection
information is:

-PROJ_INFO-------------------------------------------------
name : GCS_Everest_India_Nepal
a : 6377301.243
es : 0.006637846068423442
proj : ll
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : degree
units : degrees
meters : 1.0

Question 1: Now, I can import a shapefile or raster layer in
latlon/WGS84 (EPGS: 4326) without reprojection. GRASS thus assumes the
projection is identical, even if they are not? Why is that?

The proj.4 definition above is missing the shift w.r.t. the WGS84.

That would be why: if no towgs84 info is given, the projection is assumed as using wgs84 as datum.

The
correct proj.4 definition for EPGS 37203 would be [2]

+proj=longlat +a=6377301.243 +b=6356100.2284
+towgs84=295,736,257,0,0,0,0 +no_defs

Question 2: In QGIS it is possible to define a custom CRS using the
above. Can I do this in GRASS as well, and if so, how?

The location creation wizard allows you to define a location using a proj.4 string.

Or can / should
I
add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Not sure about that file, but you could also add it directly in the /use/share/project/epsg or equivalent.

Moritz

Regards,

Paulo

[1] https://epsg.io/37203
[2]
https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

4

Or can / should I
add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

quite interesting as I don't have here $GISBASE/etc/proj/ogr_csv

from a fresh svn up and compilation

nada:~/dev/cpp/grass7_trunk/dist.x86_64-pc-linux-gnu/etc/proj$ ls -l
insgesamt 228
-rw-r--r-- 1 5725 Mär 12 22:33 datum.table
-rw-r--r-- 1 7556 Mär 12 22:33 datumtransform.table
-rw-r--r-- 1 1465 Mär 12 22:34 desc.table
-rw-r--r-- 1 5969 Mär 12 22:33 ellipse.table
-rw-r--r-- 1 8292 Mär 12 22:33 ellipse.table.solar.system
-rw-r--r-- 1 120559 Mär 12 22:33 FIPS.code
-rw-r--r-- 1 7939 Mär 12 22:34 parms.table
-rw-r--r-- 1 2950 Mär 12 22:33 projections
-rw-r--r-- 1 25112 Mär 12 22:33 state27
-rw-r--r-- 1 20617 Mär 12 22:33 state83
-rw-r--r-- 1 464 Mär 12 22:34 units.table

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/custom-CRS-tp5312992p5313012.html
Sent from the Grass - Dev mailing list archive at Nabble.com.

5

On 18-03-17 13:52, Helmut Kudrnovsky wrote:

Or can / should I
add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

quite interesting as I don't have here $GISBASE/etc/proj/ogr_csv

Same here, sorry, copied that line from the g.proj help file. Should probably be corrected in the g.proj help file (in the description section in the part about the epsg. But the path may differ depending on the OS I guess, so not sure what to change it to.

from a fresh svn up and compilation

nada:~/dev/cpp/grass7_trunk/dist.x86_64-pc-linux-gnu/etc/proj$ ls -l
insgesamt 228
-rw-r--r-- 1 5725 Mär 12 22:33 datum.table
-rw-r--r-- 1 7556 Mär 12 22:33 datumtransform.table
-rw-r--r-- 1 1465 Mär 12 22:34 desc.table
-rw-r--r-- 1 5969 Mär 12 22:33 ellipse.table
-rw-r--r-- 1 8292 Mär 12 22:33 ellipse.table.solar.system
-rw-r--r-- 1 120559 Mär 12 22:33 FIPS.code
-rw-r--r-- 1 7939 Mär 12 22:34 parms.table
-rw-r--r-- 1 2950 Mär 12 22:33 projections
-rw-r--r-- 1 25112 Mär 12 22:33 state27
-rw-r--r-- 1 20617 Mär 12 22:33 state83
-rw-r--r-- 1 464 Mär 12 22:34 units.table

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/custom-CRS-tp5312992p5313012.html
Sent from the Grass - Dev mailing list archive at Nabble.com.
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

6

On 18-03-17 16:30, Paulo van Breugel wrote:

On 18-03-17 13:52, Helmut Kudrnovsky wrote:

Or can / should I
add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

quite interesting as I don't have here $GISBASE/etc/proj/ogr_csv

Same here, sorry, copied that line from the g.proj help file. Should probably be corrected in the g.proj help file (in the description section in the part about the epsg. But the path may differ depending on the OS I guess, so not sure what to change it to.

OK, I am confusing things: I have the datumtransform.table and the datum.table in two folders:

/usr/local/grass7/grass-7.3.svn/etc/proj/
/usr/local/share/gdal/grass/etc/ (installed with grass-gdal plugin?)

The other files listed below are only in the folder /usr/local/grass7/grass-7.3.svn/etc/proj/. So the help file is correct in my case.

from a fresh svn up and compilation

nada:~/dev/cpp/grass7_trunk/dist.x86_64-pc-linux-gnu/etc/proj$ ls -l
insgesamt 228
-rw-r--r-- 1 5725 Mär 12 22:33 datum.table
-rw-r--r-- 1 7556 Mär 12 22:33 datumtransform.table
-rw-r--r-- 1 1465 Mär 12 22:34 desc.table
-rw-r--r-- 1 5969 Mär 12 22:33 ellipse.table
-rw-r--r-- 1 8292 Mär 12 22:33 ellipse.table.solar.system
-rw-r--r-- 1 120559 Mär 12 22:33 FIPS.code
-rw-r--r-- 1 7939 Mär 12 22:34 parms.table
-rw-r--r-- 1 2950 Mär 12 22:33 projections
-rw-r--r-- 1 25112 Mär 12 22:33 state27
-rw-r--r-- 1 20617 Mär 12 22:33 state83
-rw-r--r-- 1 464 Mär 12 22:34 units.table

-----
best regards
Helmut
--
View this message in context: http://osgeo-org.1560.x6.nabble.com/custom-CRS-tp5312992p5313012.html
Sent from the Grass - Dev mailing list archive at Nabble.com.
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

7

On 18-03-17 12:26, Moritz Lennert wrote:

Le 18 mars 2017 10:06:09 GMT+01:00, Paulo van Breugel <p.vanbreugel@gmail.com> a écrit :

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in
the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

The corresponding proj.4 definition (as defined in the esri epsg-codes
file in /usr/share/proj/)

+proj=longlat +a=6377301.243 +b=6356100.230165384 +no_defs

And when I create a location based on this shapefile, the projection
information is:

-PROJ_INFO-------------------------------------------------
name : GCS_Everest_India_Nepal
a : 6377301.243
es : 0.006637846068423442
proj : ll
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : degree
units : degrees
meters : 1.0

Question 1: Now, I can import a shapefile or raster layer in
latlon/WGS84 (EPGS: 4326) without reprojection. GRASS thus assumes the
projection is identical, even if they are not? Why is that?

The proj.4 definition above is missing the shift w.r.t. the WGS84.

That would be why: if no towgs84 info is given, the projection is assumed as using wgs84 as datum.

Ah, yes. Clearly I am not fully understanding how this projecting is handled. I was looking at the ellipsoid parameters, a and es which are different from EPSG 4326.

The
correct proj.4 definition for EPGS 37203 would be [2]

+proj=longlat +a=6377301.243 +b=6356100.2284
+towgs84=295,736,257,0,0,0,0 +no_defs

Question 2: In QGIS it is possible to define a custom CRS using the
above. Can I do this in GRASS as well, and if so, how?

The location creation wizard allows you to define a location using a proj.4 string.

Of course, overlooked that in the rush of the moment. Thanks.

Or can / should
I
add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Not sure about that file, but you could also add it directly in the /use/share/project/epsg or equivalent.

Ah ok, thanks. Will this work properly when reprojecting layers from other locations with datum transformation? I am asking because in the help file of r.proj it is mentioned that it supports general datum transformations, making use of the PROJ.4 co-ordinate system translation library. If I understand correctly the datum.table and datumtransform.table files in the $GISBASE/etc/proj/ogr_csv folder? And those tables do not include the concerning datum / transformation information.

Moritz

Regards,

Paulo

[1] https://epsg.io/37203
[2]
https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689
_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

8

On 18/03/17 12:15, Paulo van Breugel wrote:

On 18-03-17 10:06, Paulo van Breugel wrote:

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in
the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

[SNIP]

+proj=longlat +a=6377301.243 +b=6356100.2284
+towgs84=295,736,257,0,0,0,0 +no_defs

Are you sure that your "+towgs84" parameters are correct?
It seems to me that you intend to apply a three parameter
transformation, but you are specifying a seven parameter
transform with the scaling (last parameter) set to "0",
which might have unintended consequences.

To specify a three parameter datum transform:

+towgs84=295,736,257

If you really want to specify all seven parameters, but
want only the first three to take effect:

+towgs84=295,736,257,0,0,0,1

Best,

Ben

Question 2: In QGIS it is possible to define a custom CRS using the
above. Can I do this in GRASS as well, and if so, how? Or can / should
I add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Regards,

Paulo

[1] https://epsg.io/37203
[2]
https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689

A link that provides further relevant information

https://github.com/klokantech/epsg.io/issues/49

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

--
Dr. Benjamin Ducke
{*} Geospatial Consultant
{*} GIS Developer

Spatial technology for the masses, not the classes:
experience free and open source GIS at http://gvsigce.org

9

On 19/03/2017 12:27, Benjamin Ducke wrote:

On 18/03/17 12:15, Paulo van Breugel wrote:

On 18-03-17 10:06, Paulo van Breugel wrote:

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in
the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

[SNIP]

+proj=longlat +a=6377301.243 +b=6356100.2284
+towgs84=295,736,257,0,0,0,0 +no_defs

Are you sure that your "+towgs84" parameters are correct?
It seems to me that you intend to apply a three parameter
transformation, but you are specifying a seven parameter
transform with the scaling (last parameter) set to "0",
which might have unintended consequences.

Good question. I am not sure, I have copied that from https://github.com/klokantech/epsg.io/issues/49 (the same is given in [2]). This is not my area of expertise, so any tips (for dummies) how I can test or find out more about this would be great. I have just send an email forwarding your question to the person that proposed the proj definition, so hopefully I'll get some useful feedback.

  To specify a three parameter datum transform:

+towgs84=295,736,257

If you really want to specify all seven parameters, but
want only the first three to take effect:

+towgs84=295,736,257,0,0,0,1

Best,

Ben

Question 2: In QGIS it is possible to define a custom CRS using the
above. Can I do this in GRASS as well, and if so, how? Or can / should
I add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Regards,

Paulo

[1] https://epsg.io/37203
[2]
https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689

A link that provides further relevant information

https://github.com/klokantech/epsg.io/issues/49

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

10

On 19/03/17 20:11, Paulo van Breugel wrote:

On 19/03/2017 12:27, Benjamin Ducke wrote:

On 18/03/17 12:15, Paulo van Breugel wrote:

On 18-03-17 10:06, Paulo van Breugel wrote:

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in
the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

[SNIP]

+proj=longlat +a=6377301.243 +b=6356100.2284
+towgs84=295,736,257,0,0,0,0 +no_defs

Are you sure that your "+towgs84" parameters are correct?
It seems to me that you intend to apply a three parameter
transformation, but you are specifying a seven parameter
transform with the scaling (last parameter) set to "0",
which might have unintended consequences.

Good question. I am not sure, I have copied that from
https://github.com/klokantech/epsg.io/issues/49 (the same is given in
[2]). This is not my area of expertise, so any tips (for dummies) how I
can test or find out more about this would be great.

Because your source system does not use the WGS84
datum, you should specify a datum transform. This is
assuming that you wish to reproject your data at some
point from the Indian Grid System into something more
modern, such as UTM with WGS84 datum. If you will
never need to reproject, then you can just ignore all
of the following.

Re. the datum transformation:
You'll know there is a problem when you get a result
with rather lousy accuracy after reprojection. I would
expect a simple three parameter datum transform to give
you ca. 5m of result accuracy.

Simple datum transformations have three parameters:

1-3: x/y/z translations

... or an additional four (more accurate):

4-6: x/y/z rotations
7: scaling

It's a common mistake to set the seventh parameter
to "0", believing that this means "skip this"
(and by "common" I mean: I've done this myself before).
This is true for parameters 1-6 (translating or rotating
by "0" does "nothing"). But it is not true for 7!
Scaling a number by "0" sets its value to "0". Scaling
by "1" is the operation that "does nothing".

To be on the safe side, specify only the three parameters
that your datum transform actually has:

+towgs84=295,736,257

I have just send an
email forwarding your question to the person that proposed the proj
definition, so hopefully I'll get some useful feedback.

This seems to be a rather nice page:

http://kanchu_deep.tripod.com/igintro.htm

It has the three parameters for your datum transform
plus a bunch of very useful literature at the bottom
of the page.

Best,

Ben

  To specify a three parameter datum transform:

+towgs84=295,736,257

If you really want to specify all seven parameters, but
want only the first three to take effect:

+towgs84=295,736,257,0,0,0,1

Best,

Ben

Question 2: In QGIS it is possible to define a custom CRS using the
above. Can I do this in GRASS as well, and if so, how? Or can / should
I add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Regards,

Paulo

[1] https://epsg.io/37203
[2]
https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689

A link that provides further relevant information

https://github.com/klokantech/epsg.io/issues/49

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev

--
Dr. Benjamin Ducke
{*} Geospatial Consultant
{*} GIS Developer

Spatial technology for the masses, not the classes:
experience free and open source GIS at http://gvsigce.org

11

On 19-03-17 21:24, Benjamin Ducke wrote:

On 19/03/17 20:11, Paulo van Breugel wrote:

On 19/03/2017 12:27, Benjamin Ducke wrote:

On 18/03/17 12:15, Paulo van Breugel wrote:

On 18-03-17 10:06, Paulo van Breugel wrote:

Dear all,

I have a shapefile which is in CRS ESRI:37203 [1], which is defined in
the prj file as

GEOGCS["GCS_Everest_India_Nepal",DATUM["D_Everest_India_Nepal",SPHEROID["Everest_Definition_1962",6377301.243,300.8017255]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

[SNIP]

+proj=longlat +a=6377301.243 +b=6356100.2284
+towgs84=295,736,257,0,0,0,0 +no_defs

Are you sure that your "+towgs84" parameters are correct?
It seems to me that you intend to apply a three parameter
transformation, but you are specifying a seven parameter
transform with the scaling (last parameter) set to "0",
which might have unintended consequences.

Good question. I am not sure, I have copied that from
https://github.com/klokantech/epsg.io/issues/49 (the same is given in
[2]). This is not my area of expertise, so any tips (for dummies) how I
can test or find out more about this would be great.

Because your source system does not use the WGS84
datum, you should specify a datum transform. This is
assuming that you wish to reproject your data at some
point from the Indian Grid System into something more
modern, such as UTM with WGS84 datum. If you will
never need to reproject, then you can just ignore all
of the following.

Re. the datum transformation:
You'll know there is a problem when you get a result
with rather lousy accuracy after reprojection. I would
expect a simple three parameter datum transform to give
you ca. 5m of result accuracy.

Simple datum transformations have three parameters:

1-3: x/y/z translations

... or an additional four (more accurate):

4-6: x/y/z rotations
7: scaling

It's a common mistake to set the seventh parameter
to "0", believing that this means "skip this"
(and by "common" I mean: I've done this myself before).
This is true for parameters 1-6 (translating or rotating
by "0" does "nothing"). But it is not true for 7!
Scaling a number by "0" sets its value to "0". Scaling
by "1" is the operation that "does nothing".

To be on the safe side, specify only the three parameters
that your datum transform actually has:

+towgs84=295,736,257

I have just send an
email forwarding your question to the person that proposed the proj
definition, so hopefully I'll get some useful feedback.

This seems to be a rather nice page:

http://kanchu_deep.tripod.com/igintro.htm

It has the three parameters for your datum transform
plus a bunch of very useful literature at the bottom
of the page.

This is very useful information, thank you!

Best,

Ben

   To specify a three parameter datum transform:

+towgs84=295,736,257

If you really want to specify all seven parameters, but
want only the first three to take effect:

+towgs84=295,736,257,0,0,0,1

Best,

Ben

Question 2: In QGIS it is possible to define a custom CRS using the
above. Can I do this in GRASS as well, and if so, how? Or can / should
I add this to the datum information files in $GISBASE/etc/proj/ogr_csv?

Regards,

Paulo

[1] https://epsg.io/37203
[2]
https://disqus.com/home/discussion/qgistutorials/georeferencing_topo_sheets_and_scanned_maps/#comment-1426454689

A link that provides further relevant information

https://github.com/klokantech/epsg.io/issues/49

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev