[PROJ] Always using tinshift when dealing with KKJ
Even Rouault
even.rouault at spatialys.com
Fri Aug 20 06:27:46 PDT 2021
Oskari,
I assume that you want to use the same exact transformation between KKJ
and WGS84 than between KKJ and ETRS89, that is you assume that WGS84 =
ETRS89 for your purpose. In EPSG, transformation EPSG:1149 does that
with an advertized accuracy of 1m (the issue here is that WGS 84 is a
time dependent datum and ETRS 89 is fixed to the plate, hence things are
drifting over time...)
Anyway if you're happy with that assumption,
you can see in
https://github.com/OSGeo/PROJ/blob/master/data/sql/other_transformation_custom.sql#L29
that the transformation between EPSG:4123 (KKJ) and EPSG:4258 (ETRS89)
is defined as a concatenated operation:
INSERT INTO "concatenated_operation" VALUES('PROJ','KKJ_TO_ETRS89','KKJ
to ETRS89 (using PROJ:YKJ_TO_ETRS35FIN)','Transformation based on a
triangulated irregular network','EPSG','4123','EPSG','4258',NULL,NULL,0);
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_ETRS89',1,'EPSG','18193');
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_ETRS89',2,'PROJ','YKJ_TO_ETRS35FIN');
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_ETRS89',3,'EPSG','16065');
INSERT INTO "usage" VALUES(
'PROJ',
'KKJ_TO_ETRS89_USAGE',
'concatenated_operation',
'PROJ',
'KKJ_TO_ETRS89',
'EPSG','3333', -- extent
'EPSG','1024' -- unknown
);
For KKJ to WGS84, we can restart from that and add a final 4th step for
the ETRS89 to WGS84 transformation:
So with the sqlite3 commandline utility:
sqlite3 proj.db
and then type
INSERT INTO "concatenated_operation" VALUES('PROJ','KKJ_TO_WGS84','KKJ
to WGS84 (using PROJ:YKJ_TO_ETRS35FIN and EPSG:1149)','Transformation
based on a triangulated irregular network for KKJ to ETRS89, and then
assuming ETRS89 ~= WGS84','EPSG','4123','EPSG','4326',NULL,NULL,0);
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_WGS84',1,'EPSG','18193');
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_WGS84',2,'PROJ','YKJ_TO_ETRS35FIN');
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_WGS84',3,'EPSG','16065');
INSERT INTO "concatenated_operation_step"
VALUES('PROJ','KKJ_TO_WGS84',4,'EPSG','1149');
INSERT INTO "usage" VALUES(
'PROJ',
'KKJ_TO_WGS84_USAGE',
'concatenated_operation',
'PROJ',
'KKJ_TO_WGS84',
'EPSG','3333', -- extent
'EPSG','1024' -- unknown
);
you'll aso need to degrade the accuracy of the registered transformation
between KKJ and WGS84 from 1m to 1.5m, as the above transformation will
have an accuracy of 1.1m (0.1 m for KKJ to ETRS89 and 1m for ETRS89 to
WGS84), so that our new transformation is picked up by default:
UPDATE helmert_transformation_table set accuracy=1.5 where
auth_name='EPSG' and code = 10099;
With that in place:
$ PROJ_LIB=.:../../PROJ-data/fi_nls src/projinfo -s EPSG:2393 -t
EPSG:4326 -o proj
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, KKJ / Finland Uniform Coordinate System to ETRS35FIN +
TM35FIN + ETRS89 to WGS 84 (1), 1.1 m, Finland - onshore.
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=tinshift +file=fi_nls_ykj_etrs35fin.json
+step +proj=utm +zone=35 +ellps=GRS80
-------------------------------------
Operation No. 2:
unknown id, Inverse of Finland Uniform Coordinate System + KKJ to WGS 84
(2), 1.5 m, Finland - onshore.
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +inv +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0
+ellps=intl
+step +proj=push +v_3
+step +proj=cart +ellps=intl
+step +proj=helmert +x=-96.062 +y=-82.428 +z=-121.753 +rx=-4.801
+ry=-0.345
+rz=1.376 +s=1.496 +convention=coordinate_frame
+step +inv +proj=cart +ellps=WGS84
+step +proj=pop +v_3
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
And ... dang ... this is wrong... The pipeline of operation no 1 is
incorrect. It does the UTM conversion in the right direction. It should
rather be (as output by projinfo -s EPSG:2393 -t EPSG:4258)
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=tinshift +file=fi_nls_ykj_etrs35fin.json
+step +inv +proj=utm +zone=35 +ellps=GRS80
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
I've filed this issue as https://github.com/OSGeo/PROJ/issues/2817. The
issue here is that in the concatenated_operation_step table we only
indicate the operation, but not its direction, and it is quite hard for
PROJ to figure out in some cases if it must be followed in the forward
or reverse direction...
So let's, start again from the fresh proj.db and add the 2 following
records using a PROJ string based transformation:
INSERT INTO other_transformation VALUES(
'PROJ','KKJ_TO_WGS84',
'KKJ to WGS84 (using PROJ:YKJ_TO_ETRS35FIN and EPSG:1149)',
'Transformation based on a triangulated irregular network for KKJ
to ETRS89, and then assuming ETRS89 ~= WGS84',
'PROJ','PROJString',
'+proj=pipeline ' ||
'+step +proj=axisswap +order=2,1 ' ||
' +step +proj=unitconvert +xy_in=deg +xy_out=rad ' ||
'+step +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0
+ellps=intl ' ||
'+step +proj=tinshift +file=fi_nls_ykj_etrs35fin.json ' ||
'+step +inv +proj=utm +zone=35 +ellps=GRS80 ' ||
'+step +proj=unitconvert +xy_in=rad +xy_out=deg ' ||
'+step +proj=axisswap +order=2,1',
'EPSG','4123',
'EPSG','4326',
0.99,
NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,0);
INSERT INTO "usage"
VALUES('PROJ','KKJ_TO_WGS84_USAGE','other_transformation','PROJ','KKJ_TO_WGS84','EPSG','3333','EPSG','1031');
(no need to tweak EPSG:10099 accuracy, as we can directly provide the
one for our transformation)
Now, we get the expected result:
$ PROJ_LIB=.:../../PROJ-data/fi_nls src/projinfo -s EPSG:2393 -t
EPSG:4326 -o proj
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, Inverse of Finland Uniform Coordinate System + KKJ to WGS84
(using PROJ:YKJ_TO_ETRS35FIN and EPSG:1149), 0.99 m, Finland - onshore.
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=tinshift +file=fi_nls_ykj_etrs35fin.json
+step +inv +proj=utm +zone=35 +ellps=GRS80
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
-------------------------------------
Operation No. 2:
unknown id, Inverse of Finland Uniform Coordinate System + KKJ to WGS 84
(2), 1 m, Finland - onshore.
~~~
Instead of hacking proj.db itself, you can indeed create an auxiliary
database
- Initiate it with : projinfo --dump-db-structure | sqlite3 aux.db
- then remove a trigger that does a check that would fail (the trigger
would check that EPSG:4123 exists, but it exists only in proj.db, not in
aux.db) : sqlite3 aux.db "DROP TRIGGER other_transformation_insert_trigger;"
- and then redo the 2 above INSERT INTO other_transformation and INSERT
INTO "usage" in aux.db itself
And then set the PROJ_AUX_DB environment variable to point to aux.db
$ PROJ_AUX_DB=aux.db PROJ_LIB=data:../../PROJ-data/fi_nls src/projinfo
-s EPSG:2393 -t EPSG:4326 -o proj
Candidate operations found: 2
-------------------------------------
Operation No. 1:
unknown id, Inverse of Finland Uniform Coordinate System + KKJ to WGS84
(using PROJ:YKJ_TO_ETRS35FIN and EPSG:1149), 0.99 m, Finland - onshore.
PROJ string:
+proj=pipeline
+step +proj=axisswap +order=2,1
+step +proj=tinshift +file=fi_nls_ykj_etrs35fin.json
+step +inv +proj=utm +zone=35 +ellps=GRS80
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1
-------------------------------------
Operation No. 2:
unknown id, Inverse of Finland Uniform Coordinate System + KKJ to WGS 84
(2), 1 m, Finland - onshore.
Best regards,
Even
Le 20/08/2021 à 13:00, Oskari Timperi a écrit :
> projinfo -s EPSG:2393 -t EPSG:4326 -o proj
--
http://www.spatialys.com
My software is free, but my time generally not.
More information about the PROJ
mailing list