[Gvsig_english] reproject swiss epsg 21781 to wgs84 with gvsig

Thu Dec 17 10:52:07 CET 2009

Hello, Bernie. I have done this test:

- According to [1], these coordinates correspond to the same point:

EPSG:4326
LON = 8
LAT = 47
H = 900
(I took the height from Google Maps. The height above sea level is not the ellipsoidal height, but I think it's a better approximation than zero)

EPSG:21781
EAST = 642695
NORTH = 205591
H = 851

- I have written a little CSV file like this:

EAST;NORTH;LON;LAT
642695;205591;8;47

and created two shapefiles 4326.shp and 21781.shp extracting the data from those columns.

- According to [2] the WKT for EPSG:21781 is:

PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["CH1903",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY
["EPSG","7004"]],TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6149"]],PRIMEM["Greenwich",0,AUTHORITY
["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4149"]],UNIT
["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER
["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER
["azimuth",89.9999],PARAMETER["rectified_grid_angle",89.9999],PARAMETER["scale_factor",1],PARAMETER
["false_easting",600000],PARAMETER["false_northing",200000],AUTHORITY["EPSG","21781"],AXIS["Y",EAST],AXIS["X",NORTH]]

(this corresponds to the Swiss national coordinate system, I think, and I have replaced 90 with 89.9999)

- Then in gvSIG create a custom CRS with the previous string (using the WKT option to describe a CRS)

- Create a view in EPSG:4326 and add 4326.shp

- Add the file 21781.shp indicating the custom CRS, and in the reprojection option, choose "Manual reprojection". Then you have to enter the parameters to transform from CH1903 to WGS84. As we can see in the WKT string, the seven parameters are:

674.374
15.056
405.346
0
0
0
0

- Then measure the distance between the two resulting points. In my test, the distance is 0.62 m.

[1] http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html
[2] http://spatialreference.org/ref/epsg/21781/ogcwkt/

Regards,

Juan Lucas Domínguez Rubio
---
Prodevelop SL, Valencia (España)
Tlf.: 96.351.06.12 -- Fax: 96.351.09.68
http://www.prodevelop.es <http://www.prodevelop.es/> 
---

________________________________

De: gvsig_internacional-bounces at listserv.gva.es en nombre de Juan Lucas Dominguez Rubio
Enviado el: jue 17/12/2009 0:21
Para: Users and Developers mailing list; gvsig_internacional at listserv.gva.es
Asunto: Re: [Gvsig_english] reproject swiss epsg 21781 to wgs84 with gvsig

Hello, Bernie:

According to this message, the grid chenxy06.gsb does not take you from EPSG:21781 to EPSG:4326. It takes you from the CH1903 datum to the CH1903+ datum:

http://lists.maptools.org/pipermail/proj/2009-January/004296.html

Some weeks ago, I read something about this issue and somebody proposed to use 89.99999999 instead of 90, but you have used 89.9 which maybe is not close enough to 90. I don't have gvSIG 1.9 now, I'll do some other test tomorrow.

Regards,

Juan Lucas Domínguez Rubio
---
Prodevelop SL, Valencia (España)
Tlf.: 96.351.06.12 -- Fax: 96.351.09.68
http://www.prodevelop.es <http://www.prodevelop.es/> 
---

________________________________

De: gvsig_internacional-bounces at listserv.gva.es en nombre de Bernhard Draeyer
Enviado el: mié 16/12/2009 11:24
Para: gvsig_internacional at listserv.gva.es
Asunto: [Gvsig_english] reproject swiss epsg 21781 to wgs84 with gvsig

hello

i try to reproject a shp-layer with the swiss epsg 21781 coordinate system to the epsg 4326 (wgs84) system. unfortunately the geometry is displayed some 200 m northeast from where it should be. i couldn't find out so far if the transformation doesn't work properly or if i do some wrong manipulation. this is how i proceeded so far:

1. create a user crs from a copy of the 21781 projected swiss grid. set azimuth from 90 to 89.9 and rectified height to 89.9 (since the original 21781 gives an error, when these params are set to 90).

2. open a view and apply the user crs to it.

3. add the shp-file to the layer. it displays with correct swiss coordinates.

4. open the geoprocess manager and select reproject. source prj is my user crs. target prj is epsg 4326 wgs84. there exists no epsg transformation in the table so the transformation was done by grids->chenxy06.gsb (which can be downloaded from the webpage of the swiss federal office of topography). a layer "..wgs84" is added to the view.

5. export the "..wgs84" layer as shp.

6. import the "..wgs84" shp in a new view with epsg4326 set as view crs. (could't figure out if i could do this step without exporting and importing the shp).

7. the new shp is displayed about 200 m to the northeast, where it should be. how do i know? because transformations with other tools give me correct results. for example the same shp transformed with arc map displays correct. also the transformation with ogr2ogr from fwtools works fine.

so where is the error?
A) i do not use correctly the geoprocess manager.
B) there is some transformation problem with the swiss coordinate system.

what looks strange to me:

fwtools uses the same gdal library as gvsig, right? so i compared some files. it seems, that the crs definitions are stored in the gcs.csv and pcs.csv files. these files in the gvsig folder are older, than the ones in the fwtools folder. in addition there exists a gcs.overrdide.csv in the fwtools with a replacement line for the 4149 CH1903 definition, which is referenced for the epsg 21781. it reads:

4149,CH1903,6149,CH1903,6149,9122,7004,8901,1,0,9603,674.374,15.056,405.346,,,,

compared to the line with gvsig:

4149,CH1903,6149,9122,7004,8901,1,0,,,,,,,,

so i replaced the epsg 4149 line in the gcs.csv from gvsig and also replaced the azimuth value from 90 to 89.9 in the epsg21781 line.
-> there is no difference when reprojecting

->>> this let me assume, that probably gvsig doesn't use the gcs.csv and pcs.csv files on startup but uses perhaps a seperate database.

in this case, is the database accessible for changing entries or can it be forced to (re-?) read the csv.

now i'm out of ideas ... somebody here to help cracking the mystery or who had similar problems?

many thanks

bernie

_________________________________________________________________
Windows Live: Make it easier for your friends to see what you're up to on Facebook.
http://www.microsoft.com/middleeast/windows/windowslive/see-it-in-action/social-network-basics.aspx?ocid=PID23461::T:WLMTAGL:ON:WL:en-xm:SI_SB_2:092009
_______________________________________________
Gvsig_internacional mailing list
Gvsig_internacional at listserv.gva.es
http://listserv.gva.es/cgi-bin/mailman/listinfo/gvsig_internacional

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://listserv.gva.es/pipermail/gvsig_internacional/attachments/20091217/7ed26aeb/attachment.htm 