<HTML dir=ltr><HEAD><TITLE>[Gvsig_english] reproject swiss epsg 21781 to wgs84 with gvsig</TITLE>
<META http-equiv=Content-Type content="text/html; charset=unicode">
<META content="MSHTML 6.00.6000.16945" name=GENERATOR></HEAD>
<BODY>
<DIV id=idOWAReplyText20265 dir=ltr>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>Hello, Bernie. I have done this test:</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- According to [1], these coordinates correspond to the same point:</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>EPSG:4326</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>LON = 8</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>LAT = 47</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>H = 900</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>(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 </FONT><FONT face="Times New Roman" color=#000000 size=3>zero)</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>EPSG:21781</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>EAST = 642695</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>NORTH = 205591</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>H = 851</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- I have written a little CSV file like this:</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>EAST;NORTH;LON;LAT</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>642695;205591;8;47</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>and created two shapefiles 4326.shp and 21781.shp extracting the data from those columns.</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- According to [2] the WKT for EPSG:21781 is:</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>PROJCS["CH1903 / LV03",GEOGCS["CH1903",DATUM["CH1903",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>["EPSG","7004"]],TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6149"]],PRIMEM["Greenwich",0,AUTHORITY</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4149"]],UNIT</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Hotine_Oblique_Mercator"],PARAMETER</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>["latitude_of_center",46.95240555555556],PARAMETER["longitude_of_center",7.439583333333333],PARAMETER</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>["azimuth",89.9999],PARAMETER["rectified_grid_angle",89.9999],PARAMETER["scale_factor",1],PARAMETER</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>["false_easting",600000],PARAMETER["false_northing",200000],AUTHORITY["EPSG","21781"],AXIS["Y",EAST],AXIS["X",NORTH]]</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>(this corresponds to the Swiss national coordinate system, I think, and I have replaced 90 with 89.9999)</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- Then in gvSIG create a custom CRS with the previous string (using the WKT option to describe a CRS)</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- Create a view in EPSG:4326 and add 4326.shp</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- 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:</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>674.374</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>15.056</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>405.346</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>0</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>0</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>0</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>0</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>- Then measure the distance between the two resulting points. In my test, the distance is 0.62 m.</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>[1] <A href="http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html">http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html</A></FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>[2] <A href="http://spatialreference.org/ref/epsg/21781/ogcwkt/">http://spatialreference.org/ref/epsg/21781/ogcwkt/</A></FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>Regards,</FONT></DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3></FONT> </DIV></DIV>
<DIV id=idSignature85352 dir=ltr>
<DIV><FONT face="Courier New" size=2><FONT face="Times New Roman" size=3>Juan Lucas Domínguez Rubio<BR></FONT>---</FONT></DIV>
<DIV><FONT face="Courier New" size=2><FONT face="Courier New" size=2><FONT face="Courier New" size=2>Prodevelop SL, Valencia (España)</FONT></DIV>
<DIV>
<DIV><FONT face="Courier New" size=2>Tlf.: 96.351.06.12 -- Fax: 96.351.09.68<BR></FONT><A href="http://www.prodevelop.es/"><FONT face="Courier New" size=2>http://www.prodevelop.es</FONT></A><BR><FONT face="Courier New" size=2>---</FONT></DIV></FONT></DIV></FONT></DIV>
<DIV dir=ltr><BR>
<HR tabIndex=-1>
<FONT face=Tahoma size=2><B>De:</B> gvsig_internacional-bounces@listserv.gva.es en nombre de Juan Lucas Dominguez Rubio<BR><B>Enviado el:</B> jue 17/12/2009 0:21<BR><B>Para:</B> Users and Developers mailing list; gvsig_internacional@listserv.gva.es<BR><B>Asunto:</B> Re: [Gvsig_english] reproject swiss epsg 21781 to wgs84 with gvsig<BR></FONT><BR></DIV>
<DIV dir=ltr>
<DIV id=idOWAReplyText24678 dir=ltr>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3>Hello, Bernie:</FONT></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr>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:</DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr><A href="http://lists.maptools.org/pipermail/proj/2009-January/004296.html">http://lists.maptools.org/pipermail/proj/2009-January/004296.html</A></DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr>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.</DIV>
<DIV dir=ltr> </DIV>
<DIV dir=ltr>Regards,</DIV>
<DIV dir=ltr><FONT face="Times New Roman" color=#000000 size=3></FONT> </DIV></DIV>
<DIV id=idSignature10972 dir=ltr>
<DIV><FONT face="Courier New" size=2><FONT face="Times New Roman" size=3>Juan Lucas Domínguez Rubio<BR></FONT>---</FONT></DIV>
<DIV><FONT face="Courier New" size=2><FONT face="Courier New" size=2><FONT face="Courier New" size=2>Prodevelop SL, Valencia (España)</FONT></DIV>
<DIV>
<DIV><FONT face="Courier New" size=2>Tlf.: 96.351.06.12 -- Fax: 96.351.09.68<BR></FONT><A href="http://www.prodevelop.es/"><FONT face="Courier New" size=2>http://www.prodevelop.es</FONT></A><BR><FONT face="Courier New" size=2>---</FONT></DIV></FONT></DIV></FONT></DIV>
<DIV dir=ltr><BR>
<HR tabIndex=-1>
<FONT face=Tahoma size=2><B>De:</B> gvsig_internacional-bounces@listserv.gva.es en nombre de Bernhard Draeyer<BR><B>Enviado el:</B> mié 16/12/2009 11:24<BR><B>Para:</B> gvsig_internacional@listserv.gva.es<BR><B>Asunto:</B> [Gvsig_english] reproject swiss epsg 21781 to wgs84 with gvsig<BR></FONT><BR></DIV>
<DIV><BR>
<P><FONT size=2>hello<BR><BR>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:<BR><BR>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).<BR><BR>2. open a view and apply the user crs to it.<BR><BR>3. add the shp-file to the layer. it displays with correct swiss coordinates.<BR><BR>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.<BR><BR>5. export the "..wgs84" layer as shp.<BR><BR>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).<BR><BR>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.<BR><BR>so where is the error?<BR>A) i do not use correctly the geoprocess manager.<BR>B) there is some transformation problem with the swiss coordinate system.<BR><BR>what looks strange to me:<BR><BR>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:<BR><BR>4149,CH1903,6149,CH1903,6149,9122,7004,8901,1,0,9603,674.374,15.056,405.346,,,,<BR><BR>compared to the line with gvsig:<BR><BR>4149,CH1903,6149,9122,7004,8901,1,0,,,,,,,,<BR><BR>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.<BR>-> there is no difference when reprojecting<BR><BR>->>> 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.<BR><BR>in this case, is the database accessible for changing entries or can it be forced to (re-?) read the csv.<BR><BR>now i'm out of ideas ... somebody here to help cracking the mystery or who had similar problems?<BR><BR>many thanks<BR><BR>bernie<BR><BR> <BR>_________________________________________________________________<BR>Windows Live: Make it easier for your friends to see what you’re up to on Facebook.<BR><A href="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">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</A><BR>_______________________________________________<BR>Gvsig_internacional mailing list<BR>Gvsig_internacional@listserv.gva.es<BR><A href="http://listserv.gva.es/cgi-bin/mailman/listinfo/gvsig_internacional">http://listserv.gva.es/cgi-bin/mailman/listinfo/gvsig_internacional</A><BR></FONT></P></DIV></DIV></BODY></HTML>