[GRASS-dev] S-JTSK strikes back

GRASS GIS GRASS GIS developer
1

Hi,

for a project I need to define a S-JTSK location. I thought
that we resolved this recently, but apparently not:

# S-JTSK (Ferro) / Krovak
<2065> +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs <>

So I do:

g.proj epsg=2065

g.proj -w
PROJCS["Krovak",
    GEOGCS["bessel",
        DATUM["unknown",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.83333333333333],
    PARAMETER["azimuth",0],
    PARAMETER["pseudo_standard_parallel_1",0],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

This is not quite right. It should be
http://lists.maptools.org/pipermail/proj/2003-February/000626.html
something like this:
PROJCS["S-JTSK (Ferro) / Krovak",
     GEOGCS["S-JTSK (Ferro)",
         DATUM["S_JTSK_Ferro",
             SPHEROID["Bessel 1841",6377397.155,299.1528128,
                 AUTHORITY["EPSG","7004"]],
             AUTHORITY["EPSG","6818"]],
         PRIMEM["Ferro",-17.66666666666667,
             AUTHORITY["EPSG","8909"]],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4818"]],
     PROJECTION["Krovak"],
     PARAMETER["latitude_of_center",49.5],
     PARAMETER["longitude_of_center",42.5],
     PARAMETER["azimuth",30.28813972222222],
     PARAMETER["pseudo_standard_parallel_1",78.5],
     PARAMETER["scale_factor",0.9999],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     AUTHORITY["EPSG","2065"]]

Is any name remapping still missing in
lib/proj/convert.c
?

Markus

2

On Thu, May 03, 2007 at 11:50:22AM +0200, Markus Neteler wrote:

Hi,

for a project I need to define a S-JTSK location. I thought
that we resolved this recently, but apparently not:

# S-JTSK (Ferro) / Krovak
<2065> +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs <>

So I do:

g.proj epsg=2065

g.proj -w
PROJCS["Krovak",
    GEOGCS["bessel",
        DATUM["unknown",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.83333333333333],
    PARAMETER["azimuth",0],
    PARAMETER["pseudo_standard_parallel_1",0],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

This is not quite right.

Mhh, not sure about my statement. Anyway, the map which I have to import
looks like this:

PROJCS["S-JTSK_Krovak",GEOGCS["GCS_S_JTSK",DATUM["D_S_JTSK",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Krovak"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Pseudo_Standard_Parallel_1",78.5],PARAMETER["Scale_Factor",0.9999],PARAMETER["Azimuth",30.28813975277778],PARAMETER["Longitude_Of_Center",24.83333333333333],PARAMETER["Latitude_Of_Center",49.5],PARAMETER["X_Scale",1.0],PARAMETER["Y_Scale",1.0],PARAMETER["XY_Plane_Rotation",0.0],UNIT["Meter",1.0]]

And:

g.proj wkt=Area_topolcianky.prj
g.proj -w
PROJCS["Krovak",
    GEOGCS["bessel",
        DATUM["unknown",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.83333333333333],
    PARAMETER["azimuth",0],
    PARAMETER["pseudo_standard_parallel_1",0],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

Doesn't seem to correspond (see pseudo_standard_parallel_1, scale_factor etc).

Markus

3

Hi,

I tried

2007/5/3, Markus Neteler <neteler@itc.it>:

# S-JTSK (Ferro) / Krovak
<2065> +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs <>

So I do:

g.proj epsg=2065

g.proj epsg=2065 -c loc=s-jtsk datum=3

g.gisenv | grep LOC
LOCATION_NAME='s-jtsk';

g.proj -w

I got different values...

g.proj -w
        DATUM["unknown",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],

        DATUM["Militar_Geographische_Institut",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["ferro",-17.666666666668],
        UNIT["degree",0.0174532925199433]],

    PARAMETER["longitude_of_center",24.833333333332], <- This should
be 42.5 (because of Ferro)

I am wonder why g.proj -w gives different values...

g.proj -jf
+proj=krovak +lat_0=49.5 +lon_0=42.5 +k=0.9999 +x_0=0 +y_0=0 +pm=ferro
+no_defs +a=6377397.155 +rf=299.1528128
+towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56 +to_meter=1

Martin

This is not quite right. It should be
http://lists.maptools.org/pipermail/proj/2003-February/000626.html
something like this:
PROJCS["S-JTSK (Ferro) / Krovak",
     GEOGCS["S-JTSK (Ferro)",
         DATUM["S_JTSK_Ferro",
             SPHEROID["Bessel 1841",6377397.155,299.1528128,
                 AUTHORITY["EPSG","7004"]],
             AUTHORITY["EPSG","6818"]],
         PRIMEM["Ferro",-17.66666666666667,
             AUTHORITY["EPSG","8909"]],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4818"]],
     PROJECTION["Krovak"],
     PARAMETER["latitude_of_center",49.5],
     PARAMETER["longitude_of_center",42.5],
     PARAMETER["azimuth",30.28813972222222],
     PARAMETER["pseudo_standard_parallel_1",78.5],
     PARAMETER["scale_factor",0.9999],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     AUTHORITY["EPSG","2065"]]

Is any name remapping still missing in
lib/proj/convert.c
?

Markus

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *

4

On Thu, May 03, 2007 at 12:09:19PM +0200, Markus Neteler wrote:

cat Area_topolcianky.prj
PROJCS["S-JTSK_Krovak",GEOGCS["GCS_S_JTSK",DATUM["D_S_JTSK",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Krovak"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Pseudo_Standard_Parallel_1",78.5],PARAMETER["Scale_Factor",0.9999],PARAMETER["Azimuth",30.28813975277778],PARAMETER["Longitude_Of_Center",24.83333333333333],PARAMETER["Latitude_Of_Center",49.5],PARAMETER["X_Scale",1.0],PARAMETER["Y_Scale",1.0],PARAMETER["XY_Plane_Rotation",0.0],UNIT["Meter",1.0]]

g.proj wkt=Area_topolcianky.prj
WARNUNG: Datum 'Jednotne_Trigonometricke_Site_Katastralni' not recognised
         by GRASS and no parameters found.
g.proj -w
PROJCS["Krovak",
    GEOGCS["bessel",
        DATUM["unknown",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",24.83333333333333],
    PARAMETER["azimuth",0],
    PARAMETER["pseudo_standard_parallel_1",0],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

It has picked up at least 0.9999 (lost pseudo_standard_parallel_1 etc).

If I do the same with OGR, I get
ogrinfo -so Area_topolcianky.shp Area_topolcianky
INFO: Open of `Area_topolcianky.shp'
      using driver `ESRI Shapefile' successful.

Layer name: Area_topolcianky
Geometry: Polygon
Feature Count: 1
Extent: (1255514.000000, 473866.000000) - (1256400.000000, 474752.000000)
Layer SRS WKT:
PROJCS["S-JTSK_Krovak",
    GEOGCS["GCS_S_JTSK",
        DATUM["Jednotne_Trigonometricke_Site_Katastralni",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Pseudo_Standard_Parallel_1",78.5],
    PARAMETER["Scale_Factor",0.9999],
    PARAMETER["Azimuth",30.28813975277778],
    PARAMETER["Longitude_Of_Center",24.83333333333333],
    PARAMETER["Latitude_Of_Center",49.5],
    PARAMETER["X_Scale",1.0],
    PARAMETER["Y_Scale",1.0],
    PARAMETER["XY_Plane_Rotation",0.0],
    UNIT["Meter",1.0]]
Id: Integer (6.0)
AREA: Real (19.11)
PERIM: Real (19.11)

Apparently GRASS is still lacking something.

Markus

5

On Thu, 3 May 2007, Markus Neteler wrote:

Hi,

for a project I need to define a S-JTSK location. I thought
that we resolved this recently, but apparently not:

# S-JTSK (Ferro) / Krovak
<2065> +proj=krovak +lat_0=49.5 +lon_0=42.5 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=ferro +units=m +no_defs <>

So I do:

g.proj epsg=2065

Hello Markus,
As far as I can remember it was only resolved in so far as we decided that you couldn't create a Krovak location using EPSG codes and it needed to be done manually using g.setproj? This was because of problems translating WKT into a PROJ.4 description inside OGR, partly related to the fact that "proj=krovak" is still kind of temporary and has various paramters hard-coded in it and a good flexible translation will require to use the new "proj=kocc" when it is eventually merged into PROJ.4. But if you think there is more to it than this then I can look into it more.

Paul

g.proj -w
PROJCS["Krovak",
   GEOGCS["bessel",
       DATUM["unknown",
           SPHEROID["Bessel_1841",6377397.155,299.1528128]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Krovak"],
   PARAMETER["latitude_of_center",49.5],
   PARAMETER["longitude_of_center",24.83333333333333],
   PARAMETER["azimuth",0],
   PARAMETER["pseudo_standard_parallel_1",0],
   PARAMETER["scale_factor",0.9999],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   UNIT["Meter",1]]

This is not quite right. It should be
http://lists.maptools.org/pipermail/proj/2003-February/000626.html
something like this:
PROJCS["S-JTSK (Ferro) / Krovak",
    GEOGCS["S-JTSK (Ferro)",
        DATUM["S_JTSK_Ferro",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            AUTHORITY["EPSG","6818"]],
        PRIMEM["Ferro",-17.66666666666667,
            AUTHORITY["EPSG","8909"]],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4818"]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",49.5],
    PARAMETER["longitude_of_center",42.5],
    PARAMETER["azimuth",30.28813972222222],
    PARAMETER["pseudo_standard_parallel_1",78.5],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2065"]]

Is any name remapping still missing in
lib/proj/convert.c
?

Markus

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

6

OK I think I see now what you are noticing. This seems like a slightly
different problem again from the one we came across last time. In general
I think the Krovak support in OGR/PROJ needs completely re-done using
proj=kocc; it is not enough to keep coming across little issues like this
and filing bug reports for OGR but anyway let's have a look...

On Thu, 3 May 2007, Markus Neteler wrote:

On Thu, May 03, 2007 at 12:09:19PM +0200, Markus Neteler wrote:

cat Area_topolcianky.prj
PROJCS["S-JTSK_Krovak",GEOGCS["GCS_S_JTSK",DATUM["D_S_JTSK",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Krovak"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Pseudo_Standard_Parallel_1",78.5],PARAMETER["Scale_Factor",0.9999],PARAMETER["Azimuth",30.28813975277778],PARAMETER["Longitude_Of_Center",24.83333333333333],PARAMETER["Latitude_Of_Center",49.5],PARAMETER["X_Scale",1.0],PARAMETER["Y_Scale",1.0],PARAMETER["XY_Plane_Rotation",0.0],UNIT["Meter",1.0]]

Saving this to a file tmp.prj and running g.proj -p wkt=tmp.prj gives

-PROJ_INFO-------------------------------------------------
name : Krovak
proj : krovak
ellps : bessel
lat_0 : 49.5
lon_0 : 24.83333333333333
alpha : 30.28813975277778
k : 0.9999
x_0 : 0
y_0 : 0
no_defs : defined
-PROJ_UNITS------------------------------------------------
unit : Meter
units : Meters
meters : 1

g.proj wkt=Area_topolcianky.prj
WARNUNG: Datum 'Jednotne_Trigonometricke_Site_Katastralni' not recognised
        by GRASS and no parameters found.
g.proj -w
PROJCS["Krovak",
   GEOGCS["bessel",
       DATUM["unknown",
           SPHEROID["Bessel_1841",6377397.155,299.1528128]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Krovak"],
   PARAMETER["latitude_of_center",49.5],
   PARAMETER["longitude_of_center",24.83333333333333],
   PARAMETER["azimuth",0],
   PARAMETER["pseudo_standard_parallel_1",0],
   PARAMETER["scale_factor",0.9999],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   UNIT["Meter",1]]

It has picked up at least 0.9999 (lost pseudo_standard_parallel_1 etc).

This goes through two translations here (WKT to GRASS format and back to WKT), so is more complicated and I will come back to it later. But note that scale_factor =0.9999 is present in the orginal WKT above. I see Pseudo Standard Parallel is gone, but again have to ask - is it important for the PROJ.4 version of proj=krovak? I really have no idea.

If I do the same with OGR, I get
ogrinfo -so Area_topolcianky.shp Area_topolcianky
INFO: Open of `Area_topolcianky.shp'
     using driver `ESRI Shapefile' successful.

Layer name: Area_topolcianky
Geometry: Polygon
Feature Count: 1
Extent: (1255514.000000, 473866.000000) - (1256400.000000, 474752.000000)
Layer SRS WKT:
PROJCS["S-JTSK_Krovak",
   GEOGCS["GCS_S_JTSK",
       DATUM["Jednotne_Trigonometricke_Site_Katastralni",
           SPHEROID["Bessel_1841",6377397.155,299.1528128]],
       PRIMEM["Greenwich",0.0],
       UNIT["Degree",0.0174532925199433]],
   PROJECTION["Krovak"],
   PARAMETER["False_Easting",0.0],
   PARAMETER["False_Northing",0.0],
   PARAMETER["Pseudo_Standard_Parallel_1",78.5],
   PARAMETER["Scale_Factor",0.9999],
   PARAMETER["Azimuth",30.28813975277778],
   PARAMETER["Longitude_Of_Center",24.83333333333333],
   PARAMETER["Latitude_Of_Center",49.5],
   PARAMETER["X_Scale",1.0],
   PARAMETER["Y_Scale",1.0],
   PARAMETER["XY_Plane_Rotation",0.0],
   UNIT["Meter",1.0]]
Id: Integer (6.0)
AREA: Real (19.11)
PERIM: Real (19.11)

Apparently GRASS is still lacking something.

The WKT description here is just an exact copy of the one you posted at the top of the e-mail there and hasn't been processed at all. So no surprises that it hasn't changed.

Therefore it seems some information (the pseudo standard parallel value) is lost when converting from WKT into GRASS format, and some more information (the alpha value) is lost when converting from GRASS format back to WKT. Is this important? I don't know anyway how important alpha and pseudo standard parallel are. Are they important for proj=krovak or does it just ignore it? It seems alpha corresonds to Azimuth in the WKT and OSRImportFromProj4() isn't converting it properly, but I really don't know how important it is for proj=krovak so wouldn't like to get too deeply involved there... Likewise for Pseudo_Standard_Parallel_1 and OSRExportToProj4().

Paul

7

Hi Paul,

2007/5/4, Paul Kelly <paul-grass@stjohnspoint.co.uk>:

OK I think I see now what you are noticing. This seems like a slightly
different problem again from the one we came across last time. In general
I think the Krovak support in OGR/PROJ needs completely re-done using
proj=kocc; it is not enough to keep coming across little issues like this
and filing bug reports for OGR but anyway let's have a look...

Agree...

[snip]

> It has picked up at least 0.9999 (lost pseudo_standard_parallel_1 etc).

This goes through two translations here (WKT to GRASS format and back to
WKT), so is more complicated and I will come back to it later. But note
that scale_factor =0.9999 is present in the orginal WKT above. I see
Pseudo Standard Parallel is gone, but again have to ask - is it important
for the PROJ.4 version of proj=krovak? I really have no idea.

The WKT description here is just an exact copy of the one you posted at
the top of the e-mail there and hasn't been processed at all. So no
surprises that it hasn't changed.

Therefore it seems some information (the pseudo standard parallel value)
is lost when converting from WKT into GRASS format, and some more
information (the alpha value) is lost when converting from GRASS format
back to WKT. Is this important? I don't know anyway how important alpha
and pseudo standard parallel are. Are they important for proj=krovak or
does it just ignore it? It seems alpha corresonds to Azimuth in the WKT and
OSRImportFromProj4() isn't converting it properly, but I really don't know
how important it is for proj=krovak so wouldn't like to get too deeply
involved there... Likewise for Pseudo_Standard_Parallel_1 and
OSRExportToProj4().

Latitude of pseudo standard parallel (78d30m) is hard-coded in
PJ_krovak.c and alpha I guess is ignored.

I see problem using g.proj epsg=2065 -p / -w :

PRIMEM["ferro",-17.666666666668]
PARAMETER["longitude_of_center",24.833333333332], <- not Ferro, but Greenwich
x
pm : ferro
lon_0 : 42.5 <- OK

Martin

Paul

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *

8

Hello Martin

On Fri, 4 May 2007, Martin Landa wrote:

information (the alpha value) is lost when converting from GRASS format
back to WKT. Is this important? I don't know anyway how important alpha
and pseudo standard parallel are. Are they important for proj=krovak or
does it just ignore it? It seems alpha corresonds to Azimuth in the WKT and
OSRImportFromProj4() isn't converting it properly, but I really don't know
how important it is for proj=krovak so wouldn't like to get too deeply
involved there... Likewise for Pseudo_Standard_Parallel_1 and
OSRExportToProj4().

Latitude of pseudo standard parallel (78d30m) is hard-coded in
PJ_krovak.c and alpha I guess is ignored.

OK thanks for clearing this up a bit. So it isn't a big problem then. Just a bit messy.

I see problem using g.proj epsg=2065 -p / -w :

PRIMEM["ferro",-17.666666666668]
PARAMETER["longitude_of_center",24.833333333332], <- not Ferro, but Greenwich

Yes this is the problem we came across before. It is currently the subject of a pending bug report in OGR: http://trac.osgeo.org/gdal/ticket/367

Paul

x
pm : ferro
lon_0 : 42.5 <- OK

Martin

Paul

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *