[OSRS-PROJ] Re: Quick 'how-to' to on PROJ.4
Petri J. Riipinen
petri.riipinen at nic.fi
Thu Jul 26 13:44:10 PDT 2001
Precedence: bulk
Reply-To: osrs-proj at remotesensing.org
Frank & others,
Ok, I almost managed to do it, except that fl_transform returns an error
code (100)
First I want to mention that I found the Finnish KKJ-definitions in the
epsg-file and they are below:
# KKJ / Finland zone 1
<2391> +proj=tmerc +lat_0=0.000000000 +lon_0=21.000000000 +k=1.000000
+x_0=1500000.000 +y_0=0.000 +ellps=intl +units=m no_defs <>
# KKJ / Finland zone 2
<2392> +proj=tmerc +lat_0=0.000000000 +lon_0=24.000000000 +k=1.000000
+x_0=2500000.000 +y_0=0.000 +ellps=intl +units=m no_defs <>
# KKJ / Finland Uniform Coordinate System
<2393> +proj=tmerc +lat_0=0.000000000 +lon_0=27.000000000 +k=1.000000
+x_0=3500000.000 +y_0=0.000 +ellps=intl +units=m no_defs <>
# KKJ / Finland zone 4
<2394> +proj=tmerc +lat_0=0.000000000 +lon_0=30.000000000 +k=1.000000
+x_0=4500000.000 +y_0=0.000 +ellps=intl +units=m no_defs <>
I'm guessing that the Uniform Coordinate System is the one that I should
use with Shapefiles and then different zones with different "vertical
stripes" of Finland, when dealing with raster images.
Anyway, the translation doesn't work...:(
I define my projections like this:
const char KKJCoordinateSystem[] = "+proj=tmerc +lat_0=0.000000000
+lon_0=27.000000000 +k=1.000000 +x_0=3500000.000 +y_0=0.000
+towgs84=-87,-98,-121 +ellps=intl +units=m";
const char WGS84CoordinateSystem[] = "+proj=latlong +datum=WGS84";
and then do:
projPJ srcProjection = pj_init_plus(WGS84CoordinateSystem);
projPJ dstProjection = pj_init_plus(KKJCoordinateSystem);
I think that I should use the WGS84 as source and KKJ as destination, when
I want to transform from GPS to KKJ???
Then I transform e.g.
double dX = 6010.570;
double dY = 2455.930;
double dZ = 0;
status = pj_transform(pSourceProj, pDestProj, 1, 1, &dX, &dY, &dZ );
After this call status = 100 and dX, dY and dZ haven't changed.
I also found some accurate test data from Finnish National Land Survey:
GPS-measurements:
1 6010.570 2455.930 GPS: FINLANDIA-HOUSE
2 6012.467 2456.127 GPS: RAILWAY YARD
3 6013.045 2458.826 GPS: VIIKKI MUSEUM
Output in KKJ:
1 6674205. 2551915. GPS: FINLANDIA-HOUSE
2 6677730. 2552047. GPS: RAILWAY YARD
3 6678839. 2554525. GPS: VIIKKI MUSEUM
Do you any idea what might be the problem? Are the values '6010.570' etc in
correct form for input to pj_transform?
- Petri
At 23:18 25.7.2001 , Frank Warmerdam wrote:
>"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
<><><><><><><><>
Petri J. Riipinen
petri.riipinen at nic.fi
<><><><><><><><>
----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.
More information about the Proj
mailing list