POINT (774135.822823993 1048589.68103062)However, correct result should be
POINT (774041.368 1048448.779)Source: Transformation service of Geoportal of Czech national mapping agency http://geoportal.cuzk.cz/Default.aspx?head_tab=sekce-01-gp&mode=TextMeta&text=wcts&menu=19
POINT (14.0 50.0000000647905)According my experience the coordinate conversion is normally done in two steps (reverse case):
1. projected coordinates in SRID 2065 are converted (by reverse Krovak projection algorithm, coord_op_id=9819) to geographic coordinates in S-JTSK/GEOGRAPHIC2D (SRID 4156), datum "Jednotne Trigonometricke Site Katastralni" (EPSG ID 6156) based on Bessel Ellipsoid (EPSG ID 7004)
2. then the geographic coordinates should be transformed using Position Vector 7-param. transformation (coord_op_id=9606) to ETRS-89 (SRID 4258)Assuming that above two steps are true in Oracle, and the projection algorithm (coord_op_id=9819) is correctly implemented (I wonder if it is, btw is it hardcoded?), then we should get the expected results.
select * from
MDSYS.SDO_COORD_OPS
where
coord_op_id=9819;
or
select * from
MDSYS.SDO_COORD_OPS
where
SOURCE_SRID = 2065 AND
TARGET_SRID = 4258;
I get zero results (I'm running Oracle Database 11g Enterprise Edition Release 11.2.0.1.0). When I run your sample transformations I'm getting exactly the same results, so the default Spatial Transformation is most likely not the one that the Czech national mapping agency is using.select * from MDSYS.SDO_COORD_OPS where source_srid in (2065,4156) and target_srid in (2065,4156) ;
select * from sdo_coord_ops where coord_op_name like '%JTSK%' or coord_op_name like '%Krovak%' ;However, I believe following two operations are those I would expect to be involved in the conversion chain (two steps mentioned above):
select * from sdo_coord_op_methods where coord_op_method_id in (9606,9819);
Now, when I try what I have indicated earlier:9606 Position Vector 7-param. transformation 1 EPSG guidance note #7.
9819 Krovak Oblique Conic Conformal 1 Research Institute for Geodesy Topography and Cartography (VUGTK); Prague.
select sdo_util.to_wktgeometry(sdo_cs.transform(sdo_cs.transform(sdo_geometry(2001,4258,sdo_point_type(14,50,null),null,null),4156),2065)) as wkt_geom from dual;I get:
..which is again, very close to my initial unexpected result.POINT (774135.825117089 1048589.67606019)
Vlad Pek wrote:Maybe try SDO_CS.DETERMINE_DEFAULT_CHAIN(2056, 4156); and see what that gives. Maybe that gives us a clue as to where this is going wrong? Since the values are 94 and 141 meter apart, I am thinking this has something to do with the parameters for the datum, but I can't find any differences. But we are not entirely sure which parameters VUGTK Prague is using, are we?
What I am not sure is the fact, that no sdo_coord_op seems to exist between SRID 2065 and 4156:
Vlad Pek wrote:Me too :-)
So, at this point I am a bit confused...
SQL*Plus: Release 11.2.0.1.0 Production on Tue Jan 29 20:34:08 2013
Copyright (c) 1982, 2010, Oracle. All rights reserved.
Connected to:
Oracle Database 11g Enterprise Edition Release 11.2.0.1.0 - Production
With the Partitioning, OLAP, Data Mining and Real Application Testing options
SQL> select '4156:' as txt, sdo_util.to_wktgeometry(sdo_cs.transform(sdo_geometry(2001,4258,sdo_point_type(14,50,null),null,null),4156)) as wkt_geom from dual
2 union all
3 select '4818:', sdo_util.to_wktgeometry(sdo_cs.transform(sdo_geometry(2001,4258,sdo_point_type(14,50,null),null,null),4818)) from dual;
TXT WKT_GEOM
-------------------------------------------------------------------------------------
4156: POINT (14.0010258818362 50.0007861342037)
4818: POINT (31.666667 49.9994109369388)
SQL>
Now this becomes very interesting: 31.666 - 14.001 = 17.665. That is a bit too close to the Prime Meridian 8909 to be coincidence. It's almost as if Oracle is not even looking at the datums in it's calculations....Vlad Pek wrote:Aha! That explains a lot. I'm not sure that those types of corrections could be implemented in Oracle.
Also, I have the answer available from VUGTK's expert, who is involved in coordinate transfomations, and who was my professor during my university studies. He said that they do not use Oracle for transfoming data, because they use own precise algorithms. For the purpose of the precise geodetic works, they need to reduce the local deformations introduced historically into the cadastral mapping data (long story), by implementing a special correction grid over the area of the Czech republic.
Vlad Pek wrote:Hmmm. To be honest, I'm not sure how to implement a projection algorithm. It's pretty well described how to define your own CS, but the projection algorithm ... I think I would consider developing my own transformation procedure. It depends a bit on how often you have to transform your geometries, but seeing as it's for INSPIRE I don't think transformations need to happen realtime, right? So that would enable you to do the transformation yourself, using the proper algorithm. Don't know how good your PL/SQL knowledge is, but I would consider it. especially since
But, if we omit these deformations, we can still use Oracle based transformation, assuming that Krovak's Projection algortihm would be implemented correctly in Oracle for our purpose
Krovak projection is specific and quite complex, it in fact involves four steps, I wonder whether it is possible to realize it simply by inserting appropriate parameters to Oracle system tables. So, in case it is not, we can provide a complete algorithm of conversion from JTSK projected coordinates to SRID:4156.It may not be possible at all, allthough the Dutch projection is in there and that is a three-step projection. But I'm not familiar enough with Oracle's way of doing things to be able to tell if the Krovak Projection is possible in oracle.
Anyway, I'm raising my hand and would be happy to help in case Oracle is willing to file a CR (what is a standard way to file a CR?, do we need to go through our local support here in the Czech republic?).I'm not sure. I know that Oracle employees visit these forums, and sometimes react and file problem reports based on discussions here. But formally I would think you would have to go through your local Czech Oracle office.