การใช้สคริปต์ในการคำนวณหรือตรวจสอบด้วยโปรแกรมอื่น
(เนื้อหาในบล็อกนี้ต่อจากตอนแรก)
ผลลัพธ์จากการคำนวณส่วนที่2
(ส่วนนี้ผู้เขียนทำขึ้นเพื่อความสะดวกในการเรียนรู้ แบ่งปัน และตรวจสอบผลการคำนวณด้วยวิธีการอื่น)
ผลลัพธ์เป็นข้อความประกอบด้วยข้อมูลที่เป็นประโยชน์และสคริปต์สำหรับการคำนวณที่พร้อมใช้งานดังนี้
Extra information:
Provided by Swatchai Kriengkraipet.
Source's srs_def= epsg:7909
Target's srs_def= epsg:7912
[(Command-line) proj v8 usage]
(1) Run this 'proj.projinfo' (assuming proj is installed globally),
you will get PROJ string that can be used to build a script to compute
the coordinates transformation.
$ projinfo -s epsg:7909 -t epsg:7912
(2) Run this 'proj.cs2cs' to get the solution:
$ echo 14.0 100.0 0.0 2004.5 | cs2cs epsg:7909 epsg:7912 -d 8
The scripts in (1) and (2) above can be used in various ways:-
- for Windows users, they can be pasted and run in OSGeo4W Shell
- on systems that installed OSGeoLive 14.0, they can be pasted and
run in standard terminal
[(Python) pyproj v3.2 usage]
##_BEGIN_script
# The minimum code for round-trip computation
# Copy the lines of code below and run in your environment
from pyproj import Transformer
# Direct
geoSrc_geoTgt = Transformer.from_crs(7909,7912)
# Inverse
geoTgt_geoSrc = Transformer.from_crs(7912,7909)
# Direct computation
# Input (lat/long/h/t): 14.0, 100.0, 0.0, 2004.5
ans = geoSrc_geoTgt.transform(xx=14.0, yy=100.0, zz=0.0, tt=2004.5)
print("The transformation solution:",ans)
# Reverse-check
bns = geoTgt_geoSrc.transform(xx=ans[0], yy=ans[1], zz=ans[2], tt=ans[3])
print("Expect values close to the input:",bns)
##_END_script
The Python script can be pasted and run with/in:-
- Python session in OSGeo4W Shell
- Kaggle's Jupyter notebook
- Python running in standard terminal of OSGeoLive 14.0's users
- any Python session that required dependencies libraries/modules
are already installed
คำอธิบายแยกเป็นส่วนๆ ตามลำดับ
ส่วนที่ 1
Source's srs_def= epsg:7909
Target's srs_def= epsg:7912
เกี่ยวกับกรอบอ้างอิงจีออเดติกที่บันทึกไว้แล้วใน Database of coordinates systems (https://epsg.io)
กรอบอ้างอิง ITRF2000 มีการบันทึกไว้เป็นรหัส epsg:7909
กรอบอ้างอิง ITRF2014 มีการบันทึกไว้เป็นรหัส epsg:7912
รหัส epsg:7909 และ epsg:7912 เป็นกรณีระบบพิกัดจีออเดติก 3 มิติ (ยังมีแบบอื่นอีก 2 แบบ คือ
จีออเดติก 2 มิติ และแบบคาร์ทีเซียน 3 มิติ ซึ่งมีรหัส epsg กำหนดไว้แล้วเช่นกัน)
รหัส epsg เหล่านี้มีประโยชน์ในการใช้เขียนประกอบในสคริปต์เพื่อคำนวณการแปลงค่าพิกัดระหว่างกรอบอ้างอิง
ด้วยโปรแกรม PROJ ซึ่งทำงานอยู่เบื้องหลังของหลายๆโปรแกรมทางด้าน Geo-processing
การใช้รหัส epsg แทนการเขียนข้อความเต็มในการเขียนสคริปต์ ช่วยให้ได้สคริปต์สั้น และลดโอกาสผิดพลาด
ตาราง 1 (รหัส EPSG ของแต่ละ ITRF)
| ITRF |
Geodetic-2D |
Geodetic-3D |
Cartesian-3D |
| ITRF88 |
8988 |
7900 |
4910 |
| ITRF89 |
8989 |
7901 |
4911 |
| ITRF90 |
8990 |
7902 |
4912 |
| ITRF91 |
8991 |
7903 |
4913 |
| ITRF92 |
8992 |
7904 |
4914 |
| ITRF93 |
8993 |
7905 |
4915 |
| ITRF94 |
8994 |
7906 |
4916 |
| ITRF96 |
8995 |
7907 |
4917 |
| ITRF97 |
8996 |
7908 |
4918 |
| ITRF2000 |
8997 |
7909 |
4919 |
| ITRF2005 |
8998 |
7910 |
4896 |
| ITRF2008 |
8999 |
7911 |
5332 |
| ITRF2014 |
9000 |
7912 |
7789 |
ส่วนที่ 2
(1) Run this 'proj.projinfo' (assuming PROJ is installed globally),
you will get PROJ string that can be used to build a script to compute
the coordinates transformation.
$ projinfo -s epsg:7909 -t epsg:7912
หมายความว่าถ้าใช้สคริปต์
projinfo -s epsg:7909 -t epsg:7912
ไปรันใน OSGeo4W Shell หรือ เทอร์มินัลของ OSGeoLive 14.0 จะได้ผลลัพธ์ที่เรียกว่า "PROJ string"
และ "WKT2:2019 string" ซึ่งอธิบายขั้นตอนการประมวลผลแบบ step-by-step ในการแปลงค่าพิกัด
ขั้นตอนเหล่านี้คือการทำงานจริงแบบที่เรียกว่า "under the hood" ในการประมวลผลด้วยโปรแกรม
PROJ
อนึ่ง "
projinfo" เป็นหนึ่งใน command line application ที่เป็นส่วนย่อยของ PROJ
ตัวอย่าง
นำไปใช้ใน OSGeo4W Shell
C:\>projinfo -s epsg:7909 -t epsg:7912
นำไปใช้ใน Kaggle Jupyter notebook
!projinfo -s epsg:7909 -t epsg:7912
ผลลัพธ์ที่ได้ (ตัดส่วนของ "WKT2:2019 string" ทิ้งไปเพราะใช้เนื้อที่มาก)
Candidate operations found: 1
Note: using '--spatial-test intersects' would bring more results (6)
-------------------------------------
Operation No. 1:
unknown id, Conversion from ITRF2000 (geog3D) to ITRF2000 (geocentric) + ITRF2000 to ITRF2014 (1) + Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D), 0.01 m, World
PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0 +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0 +ds=-0.00011 +t_epoch=2010 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1
WKT2:2019 string:
CONCATENATEDOPERATION["Conversion from ITRF2000 (geog3D) to ITRF2000 (geocentric) + ITRF2000 to ITRF2014 (1) + Conversion from ITRF2014 (geocentric) to ITRF2014 (geog3D)",
...
ผู้อ่านสามารถใช้โปรแกรม projinfo เรียกดูรายละเอียดของระบบพิกัดที่เกี่ยวข้องกับการแปลงพิกัดได้โดยรันสคริปต์ต่อไปนี้
OSGeo4W Shell
projinfo epsg:7909
projinfo epsg:7912
Kaggle Jupyter notebook
!projinfo epsg:7909
!projinfo epsg:7912
ข้อความที่เป็น PROJ string สามารถนำมาประกอบเป็นสคริปต์ที่ใช้คำนวณการแปลงด้วยโปรแกรม cct ได้ ทั้งนี้โดยการเพิ่มข้อความ
"echo 14.0 100.0 0.0 2004.5 | cct " ไว้ข้างหน้าส่วนของ PROJ string
โปรแกรม cct เป็นหนึ่งใน command line application ของ PROJ
เมื่อนำไปรันใน OSGeo4W Shell หรือ เทอร์มินัลของ OSGeoLive 14.0 จะได้ผลลัพธ์ของการแปลงค่าพิกัด
จากกรอบ ITRF2000 ไป ITRF2014
C:\>echo 14.0 100.0 0.0 2004.5 | cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +z_in=m +xy_out=rad +z_out=m +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=-0.0007 +y=-0.0012 +z=0.0261 +rx=0 +ry=0 +rz=0 +s=-0.00212 +dx=-0.0001 +dy=-0.0001 +dz=0.0019 +drx=0 +dry=0 +drz=0 +ds=-0.00011 +t_epoch=2010 +convention=position_vector +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1
14.0000 100.0000 -0.0065 2004.5000
ข้อความที่เป็น WKT2:2019 string จะเป็นรายละเอียดทุกขั้นตอนของการประมวลผลแบบ descriptive
สำหรับผู้อ่านที่ตั้งคำถามในใจว่า
- การแปลงนี้มีขั้นตอนคำนวณแก้ Plate Motion ด้วยหรือไม่
- ใช้ทรัพยากร (เช่น model file ประเภท hgrid/vgrid สำหรับค่าแก้ตำแหน่งทางราบ/ทางดิ่ง) หรือเปล่า
สามารถตรวจสอบดูใน PROJ/WKT2:2019 string อาจพบคำตอบ เช่น ถ้ามีข้อความทำนองนี้
+proj=vgridshift +grids=xxxxxx.gtx
+proj=hgridshift +grids=yyyyyy.gsb
แสดงว่ามีขั้นตอนคำนวณแก้ค่าพิกัดทางราบ/ทางดิ่งโดยใช้ข้อมูลจาก model file รวมอยู่ด้วย
ส่วนที่มาของค่าแก้มีความเกี่ยวข้องกับ Plate Motion หรือไม่นั้นเป็นอีกประเด็นหนึ่ง
อาจต้องดู meta data ของ model file (.gtx .gsb) ที่นำมาใช้เพื่อยืนยัน
อนึ่งในการแปลงจาก ITRF หนึ่ง ไปอีก ITRF หนึ่งจะไม่มีการนำค่าแก้ระดับภูมิภาค/ท้องถิ่นมาใช้ เพราะ
กรอบอ้างอิง ITRF ทุกกรอบเป็นประเภทที่เรียกว่า global ทำไว้สำหรับภารกิจแบบ whole_world
ตัวอย่างที่แสดงได้ชัดเจนว่าระบบ ITRF ไม่เหมาะกับการใช้งานระดับประเทศ/ภูมิภาคที่มีแผ่นดินไหว หรือ
Plate Motion คือ ที่ประเทศญี่ปุ่นภายหลังเกิดเหตุการณ์แผ่นดินไหว/ซึนามิที่ประเทศญี่ปุ่น ค่าพิกัด ITRF
ของตำแหน่งในประเทศทั้งหมดเสียหายอย่างสิ้นเชิง ค่าพิกัดอื่นที่เป็นระบบท้องถิ่นซึ่งมีกรอบอ้างอิงอยู่ในพื้นที่
เสียหายน้อยกว่าและแก้ไขได้ง่าย
ส่วนที่ 3
(2) Run this 'proj.cs2cs' to get the solution:
$ echo 14.0 100.0 0.0 2004.5 | cs2cs epsg:7909 epsg:7912 -d 8
หมายความว่าถ้าใช้สคริปต์
echo 14.0 100.0 0.0 2004.5 | cs2cs epsg:7909 epsg:7912 -d 8
ไปรันใน OSGeo4W Shell หรือ เทอร์มินัลของ OSGeoLive 14.0 จะได้ผลลัพธ์ของการแปลงค่าพิกัด
จากกรอบ ITRF200 ไป ITRF2014
อนึ่ง "
cs2cs" เป็นหนึ่งใน command line application ที่เป็นส่วนย่อยของ PROJ
ตัวอย่างการนำไปใช้ใน OSGeo4W Shell จะได้ผลลัพธ์แบบเดียวกับการคำนวณด้วย Web Application ของผู้เขียน
อนึ่ง ถ้าต้องการจำนวนทศนิยมในผลลัพธ์ 9 ตำแหน่ง ให้แก้เลข 8 ในสคริปต์เป็น 9
C:\>echo 14.0 100.0 0.0 2004.5 | cs2cs epsg:7909 epsg:7912 -d 8
14.00000014 100.00000000 -0.00647074 2004.5
การแปลงพิกัดที่น่าสนใจ
คำนวณย้อนกลับเพื่อตรวจสอบผลลัพธ์ได้ด้วยการรันสคริปต์
C:\>echo 14.00000014 100.0000 -0.00647074 2004.5 | cs2cs epsg:7912 epsg:7909 -d 8
คำนวณใหม่ แต่ขอผลลัพธ์เป็นค่าพิกัด ECEF Cartesian ด้วยการรันสคริปต์
C:\>echo 14.0 100.0 0.0 2004.5 | cs2cs epsg:7912 epsg:7789 -d 2
ซึ่งจะได้ผลลัพธ์
-1074863.43 6095853.42 1532981.83 2004.5
ส่วนที่ 4
เป็นโปรแกรมภาษาไพทอนที่ใช้โมดูล pyproj ช่วยในการประมวลผล
[(Python) pyproj v3.2 usage]
##_BEGIN_script
# The minimum code for round-trip computation
# Copy the lines of code below and run in your environment
from pyproj import Transformer
# Direct
geoSrc_geoTgt = Transformer.from_crs(7909,7912)
# Inverse
geoTgt_geoSrc = Transformer.from_crs(7912,7909)
# Direct computation
# Input (lat/long/h/t): 14.0, 100.0, 0.0, 2004.5
ans = geoSrc_geoTgt.transform(xx=14.0, yy=100.0, zz=0.0, tt=2004.5)
print("The transformation solution:",ans)
# Reverse-check
bns = geoTgt_geoSrc.transform(xx=ans[0], yy=ans[1], zz=ans[2], tt=ans[3])
print("Expect values close to the input:",bns)
##_END_script
หมายความว่าถ้าใช้สคริปต์ซึ่งอยู่จากบรรทัด ##_BEGIN_script ถึง ##_END_script ไปรันใน
Python session ที่ได้ติดตั้ง PyProj module แล้ว
จะได้ผลลัพธ์ของการแปลงค่าพิกัดจากกรอบ ITRF2000 ไป ITRF2014
ตัวอย่างการนำไปใช้จริง สามารถ copy/paste ได้เลย
Python 3.7.3 (default, Jan 22 2021, 20:04:44)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.24.1 -- An enhanced Interactive Python. Type '?' for help.
##_BEGIN_script
# The minimum code for round-trip computation
# Copy the lines of code below and run in your environment
from pyproj import Transformer
# Direct
geoSrc_geoTgt = Transformer.from_crs(7909,7912)
# Inverse
geoTgt_geoSrc = Transformer.from_crs(7912,7909)
# Direct computation
# Input (lat/long/h/t): 14.0, 100.0, 0.0, 2004.5
ans = geoSrc_geoTgt.transform(xx=14.0, yy=100.0, zz=0.0, tt=2004.5)
print("The transformation solution:",ans)
# Reverse-check
bns = geoTgt_geoSrc.transform(xx=ans[0], yy=ans[1], zz=ans[2], tt=ans[3])
print("Expect values close to the input:",bns)
##_END_script
The transformation solution:
(14.000000138728968, 100.00000000241214, -0.006470744498074055, 2004.5)
Expect values close to the input:
(14.000000000000005, 100.0, -9.313225746154785e-10, 2004.5)
โปรแกรมที่ให้ไว้นี้สามารถดัดแปลงให้คำนวณครั้งละหลายๆ จุดได้โดยง่าย ดังนี้
# (Continue execution from above)
# Prep a list of [lat,long,height,epoch] for 3 points
xyzt = [
[14.0, 100.0, 0.0, 2004.5],
[14.5, 99.5, 5, 2002.5],
[14.6, 99.6, 6, 2002.6],
]
# Transformation of the 3 point
for ea in xyzt:
print([eb for eb in ea])
ans = geoSrc_geoTgt.transform(xx=ea[0], yy=ea[1], zz=ea[2], tt=ea[3])
print("The transformation solution:",ans)
# Reverse-check
bns = geoTgt_geoSrc.transform(xx=ans[0], yy=ans[1], zz=ans[2], tt=ans[3])
print("Reverse-check: close to the input values?:",bns)
print()
