Re: [GRASS-dev] custom CRS
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 ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
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
___
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
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 ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
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
___
grass-dev mailing list
grass-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
On 18-03-17 12:26, Moritz Lennert wrote: Le 18 mars 2017 10:06:09 GMT+01:00, Paulo van Breugel 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 ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
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 ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
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 ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
>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. ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
Le 18 mars 2017 10:06:09 GMT+01:00, Paulo van Breugel 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 ___ grass-dev mailing list grass-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/grass-dev
Re: [GRASS-dev] custom CRS
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 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
