[OSRS-PROJ] Re: Quick 'how-to' to on PROJ.4

[Frank Warmerdam]

Precedence: bulk
Reply-To: osrs-proj at remotesensing.org

"Petri J. Riipinen" wrote:
> 
> Hi Frank,
> 
> Remember me a while ago, we exchanged some email about compiling GDAL with
> Cygwin. Well, I finally got Linux, no more Cygwin to fool around with.
> 
> I have now reached a point, where I should start integrating GPS with my
> map, that I draw from Shapefiles (reading them with shapelib).
> 
> I know that the answer must be somewhere in the PROJ.4 docs, but I guess
> I'm too lazy to go through them (they look a bit scary, IMHO).
> 
> So, could you give me a quick 'how-to' on how to use PROJ.4 to translate
> from GPS-coordinates to Finnish coordinates. I really don't need to be an
> expert in "complete world of transformations", the Finnish part is enough.
> 
> I guess that these links should contain all the information that you need
> to come up with an answer:
> http://www.nls.fi/kartta/julkaisu/kkj.html
> http://www.nls.fi/kartta/palvelut/wgs84/wgs84win_e.html
> 
> Also, I got these files for testing, they might mean more to you then they
> mean to me:
> ------ name: 20p.tab -------
> !table
> !version 300
> !charset WindowsLatin1
> !Copyright National Land Survey of Finland
> !          Helsinki, 11.05.2000, for MapInfo 4.x
> Definition Table
> File "20p.tif"
> Type "Raster"
> (3419987.5,6800012.5) (0,0) Label "Pt 1",
> (3500012.5,6800012.5) (3200,0) Label "Pt 2",
> (3500012.5,6719987.5) (3200,3200) Label "Pt 3",
> (3419987.5,6719987.5) (0,3200) Label "Pt 4"
> CoordSys Earth Projection 24,28,"m",27,0,1,3500000,0

As best as I can figure out this coordinate system is:

 Tranverse Mercator
  Central Meridian: 27dE
  Center Latitude: 0dN
  Scale Factor: 1.0
  False Easting: 3500000
  False Northing: 0

The ellipse is (28 in mapinfo) is apparently Hayford 1924, or in PROJ.4
speak "+ellps=intl". 
  
To convert from KKJ meters to lat/long values in the same ellipsoid:

warmerda at gdal[160]% proj -I +proj=tmerc +lon_0=27 +x_0=3500000 +ellps=intl
3419987.5 6800012.5
25d30'25.592"E  61d18'0.427"N

To convert back:

warmerda at gdal[161]% proj +proj=tmerc +lon_0=27 +x_0=3500000 +ellps=intl
25d30'25.592"E  61d18'0.427"N
3419987.51      6800012.50

To convert to WGS84 we need to approximate the datum shift to WGS84.  
>From my mapinfo knowledge, datum 28 is european 1950 for which mapinfo
uses the shift values -87, -98, -121.

To do datum shifts from the commandline we need to use "cs2cs" which
includes datum shifting. 

warmerda at gdal[167]% cs2cs +proj=tmerc +lon_0=27 +x_0=3500000 +towgs84=-87,-98,-121 +ellps=intl +to  +proj=latlong +datum=WGS84
3419987.5 6800012.5
25d30'22.167"E  61d17'59.479"N 16.655

Notice there is a difference of several seconds with the datum shift. 

The equivelent to "cs2cs" can be accomplished from a code by calling
pj_translate() with coordinate system definitions created with
pj_init().  The KKJ definition is:

  +proj=tmerc +lon_0=27 +x_0=3500000 +towgs84=-87,-98,-121 +ellps=intl 

To WGS84 definition is:
  +proj=latlong +datum=WGS84

To go from WGS84 to KKJ just reverse the order of the coordinate
system definitions passed in pj_transform(). 

You should check the results against an accurate map to ensure I haven't
messed something up (like reversing the sense of the towgs84 parameters
or something). 

I would also encourage you to email such requests to the PROJ.4 mailing
list.  I am going to be away for a week shortly, and it gives others the
change to help out too. 

Best regards,

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent
----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.