ใช้ node-red แสดงผลที่เรียกค้นจากฐานข้อมูล SQLite

(ต่อจาก blog ก่อนหน้านี้)

การใช้ SQLite สำหรับเก็บและจัดการข้อมูล GeoJSON (Part 3)

กรณี 2 เก็บข้อมูลเป็นรูปแบบที่คล้าย GeoJSON มากที่สุด แล้วใช้ฟังก์ชันพิเศษช่วยแปลงเป็น GeoJSON

ขอตั้งเป้าหมายว่าต้องการเก็บข้อมูลของ point feature ตามนี้
- ตัวอย่าง geojson: {"type":"Feature","properties":{"id":1001,"title":"apple","value":1.2},
  "geometry":{"type":"Point","coordinates":[100.71,13.82]}}
- ไม่เก็บข้อมูลซ้ำซ้อน เช่น id, title, value มีอยู่ภายใน geojson ของpoint feature แล้ว ต้องไม่มีที่อื่นอีก

สร้าง table ชื่อ jsonpoint ดังนี้
- มี 2 column
  - feat_ID เหมือนแบบแรก (INTEGER UNIQUE)
  - geojson เป็นชนิด json (เก็บจริงเป็น text string)

sqlite> CREATE TABLE jsonpoint (feat_ID INTEGER UNIQUE, geojson json);

ตรวจดูผลลัพธ์ (มองหา jsonpoint พบว่าอยู่ที่บรรทัด 3)
sqlite> .schema
CREATE TABLE sqlite_sequence(name,seq);
CREATE TABLE jpoint (feat_ID INTEGER UNIQUE, title TEXT, value REAL, lon REAL, lat REAL);
CREATE TABLE jsonpoint (feat_ID INTEGER UNIQUE, geojson json);

ป้อนข้อมูล
sqlite> insert into jsonpoint values(1001,json('{"type":"Feature","properties":{"id":1001,"title":"apple","value":1.2},"geometry":{"type":"Point","coordinates":[100.71,13.82]}}'));
sqlite> insert into jsonpoint values(1002,json('{"type":"Feature","properties":{"id":1002,"title":"banana","value":2.3},"geometry":{"type":"Point","coordinates":[100.72,13.83]}}'));
sqlite> insert into jsonpoint values(1003,json('{"type":"Feature","properties":{"id":1003,"title":"carrot","value":3.4},"geometry":{"type":"Point","coordinates":[100.73,13.84]}}'));

เรียกค้นดูข้อมูลที่ป้อน
sqlite> select * from jsonpoint;
1001|{"type":"Feature","properties":{"id":1001,"title":"apple","value":1.2},"geometry":{"type":"Point","coordinates":[100.71,13.82]}}
1002|{"type":"Feature","properties":{"id":1002,"title":"banana","value":2.3},"geometry":{"type":"Point","coordinates":[100.72,13.83]}}
1003|{"type":"Feature","properties":{"id":1003,"title":"carrot","value":3.4},"geometry":{"type":"Point","coordinates":[100.73,13.84]}}

สถานการณ์
- โครงสร้างข้อมูลง่าย มีเพียง 2 คอลัมน์
- ค่า attributes ซ่อนอยู่ในคอลัมน์ 2 ใช้สืบค้นแบบปกติไม่ได้
  ใช้แบบนี้ไม่ได้: select * from jsonpoint where title='apple';
- คอลัมน์ 2 มีข้อมูลที่เป็น geojson แบบสมบูรณ์ที่ต้องการ

ล้วงลึกข้อมูลในคอลัมน์ geojson โดยใช้ตัวช่วย (ฟังก์ชันจัดการ JSON)
- เรียกดูเฉพาะ title (เส้นทางเข้าถึง: properties.title)

sqlite> SELECT json_extract(geojson, '$.properties.title') from jsonpoint;
apple
banana
carrot

- เรียกดูเฉพาะ value (เส้นทางเข้าถึง: properties.value)
sqlite> SELECT json_extract(geojson, '$.properties.value') from jsonpoint;
1.2
2.3
3.4

- เรียกดูเฉพาะ geojson (เส้นทางเข้าถึง: properties.title='banana')
sqlite> SELECT geojson from jsonpoint where json_extract(geojson, '$.properties.title') = 'banana';
{"type":"Feature","properties":{"id":1002,"title":"banana","value":2.3},"geometry":{"type":"Point","coordinates":[100.72,13.83]}}

สรุป
- ฟังก์ชันจัดการ JSON ที่ใช้
  - json() ใช้กลั่นกรองตอนป้อนข้อมูลไม่ให้ผิดพลาด
  - json_extract() ใช้ล้วงลึกข้อมูลภายใน JSON ออกมาใช้งาน
- ฟังก์ชันจัดการ JSON ทำให้เราสามารถใช้ sqlite จัดการกับข้อมูล geojson ได้ไม่ยากนัก

การใช้ SQLite สำหรับเก็บและจัดการข้อมูล GeoJSON (Part 2)

กรณี 1 เก็บข้อมูลแบบปกติทั่วไป (text, number, ...) แล้วใช้คำสั่งเรียกมาประกอบกันเป็น GeoJSON เมื่อต้องการ

เราจะใช้ sqlite3 ทำงานแบบโต้ตอบเพื่อสาธิตการทำงาน
หมายเหตุ
 jsDb    เป็น database file 
 ~#      คือ system prompt
 sqlite> คือ sqlite3 prompt

ใช้คำสั่ง sqlite3 jsDb เริ่มต้นใช้งาน sqlite3 กับ database file ชื่อ jsDb
(กรณีนี้เพิ่งเริ่มงาน ยังไม่มี jsDb อยู่ก่อน จึงถูกสร้างขึ้นใหม่)

~# sqlite3 jsDb
SQLite version 3.51.2 2026-01-09 17:27:48
Enter ".help" for usage hints.

ใช้คำสั่ง CREATE TABLE สร้างตาราง (table) ชื่อ jpoint สำหรับข้อมูลจำเป็นของ geoJSON ประเภท point
sqlite> CREATE TABLE jpoint (feat_ID INTEGER UNIQUE, title TEXT, value REAL, lon REAL, lat REAL);

ใช้คำสั่ง .schema ขอดูโครงสร้างข้อมูลของตารางที่สร้างไว้
sqlite> .schema
CREATE TABLE sqlite_sequence(name,seq);
CREATE TABLE jpoint (feat_ID INTEGER UNIQUE, title TEXT, value REAL, lon REAL, lat REAL);

ใช้คำสั่ง INSERT INTO ป้อนข้อมูล
sqlite> INSERT INTO jpoint values(1001, "apple", 1.2, 100.71 , 13.82);
sqlite> INSERT INTO jpoint values(1002, "banana", 2.3, 100.72 , 13.83);
sqlite> INSERT INTO jpoint values(1003, "carrot", 3.4, 100.73 , 13.84);

ใช้คำสั่ง SELECT สืบค้นข้อมูลทั้งหมด แสดงผลแบบ list
sqlite> SELECT * from jpoint;
1001|apple|1.2|100.71|13.82
1002|banana|2.3|100.72|13.83
1003|carrot|3.4|100.73|13.84

เปลี่ยนรูปแบบการแสดงผลเป็น "table" (ปกติเป็น list)
sqlite> .mode table
sqlite> select * from jpoint;
+---------+--------+-------+--------+-------+
| feat_ID | title  | value |  lon   |  lat  |
+---------+--------+-------+--------+-------+
| 1001    | apple  | 1.2   | 100.71 | 13.82 |
| 1002    | banana | 2.3   | 100.72 | 13.83 |
| 1003    | carrot | 3.4   | 100.73 | 13.84 |
+---------+--------+-------+--------+-------+

สั่งสืบค้นข้อมูลตามเงื่อนไข WHERE แบบพื้นฐาน เลือกระเบียนข้อมูลด้วยเงื่อนไข title='banana'

sqlite> .mode list
sqlite> select * from jpoint where title = 'banana';

ผลที่ได้
1002|banana|2.3|100.72|13.83

สั่งสืบค้นข้อมูลตามเงื่อนไข WHERE โดยให้ผลออกมาเป็น geoJSON ชนิด point
- มีการใช้ฟังก์ชัน json_object() และ json_array() ช่วย
- ผลที่ได้จะเป็น text string ที่สามารถนำไป parse เป็น JSON object ได้

sqlite> SELECT json_object('type','Feature','properties',json_object('id',feat_ID,'title',title,'value',value),
'geometry',json_object('type','Point','coordinates',json_array(lon,lat))) FROM jpoint WHERE title = 'banana';

ผลที่ได้
{"type":"Feature","properties":{"id":1002,"title":"banana","value":2.3},"geometry":{"type":"Point","coordinates":[100.72,13.83]}}

ข้อสังเกต
- มีการใช้ฟังก์ชัน json เฉพาะตอนเรียกค้นข้อมูลเมื่อต้องการผลลัพธ์รูปแบบ geoJSON เท่านั้น: json_object(), json_array()
- ตัวอย่างที่แสดงเหมาะกับ point feature เพราะส่วนของ geometry เป็นข้อมูลพื้นฐาน (lon, lat)
- สำหรับ linestring หรือ polygon เป็นกรณีท้าทาย เพราะข้อมูลพิกัดตำแหน่งเป็นแบบเชิงซ้อน และอาจมีขนาดใหญ่

การใช้ฟังก์ชัน json_group_array() เพื่อเก็บ point feature ทั้งหมดมารวมกันเป็น FeatureCollection สามารถทำได้โดยสั่งประมวลผลแบบ batch processing
คำสั่งที่ใช้ค่อนข้างยาว ควรรวมกันเป็นบรรทัดเดียวก่อนใช้งาน  ถ้าใช้บ่อยให้เก็บไว้ในไฟล์ เช่น "featCollOutput.sql" และสั่งใช้งานดังนี้: .read featCollOutput.sql 
หรือ sqlite3 jsDb ".read featCollOutput.sql"

~# sqlite3 jsDb "SELECT json_object('type', 'FeatureCollection','features', json_group_array(json_object('type', 
 'Feature','properties', json_object('id', feat_ID, 'title', title, 'value', value),'geometry', json_object('type', 'Point',
 'coordinates', json_array(lon, lat))))) FROM jpoint;"
 
ผลลัพธ์
{"type":"FeatureCollection","features":[
  {"type":"Feature","properties":{"id":1001,"title":"apple","value":1.2},"geometry":{"type":"Point","coordinates":[100.71,13.82]}},
  {"type":"Feature","properties":{"id":1002,"title":"banana","value":2.3},"geometry":{"type":"Point","coordinates":[100.72,13.83]}},
  {"type":"Feature","properties":{"id":1003,"title":"carrot","value":3.4},"geometry":{"type":"Point","coordinates":[100.73,13.84]}}]}

ถ้าต้องการเก็บผลลัพธ์เป็นแฟ้มข้อมูล สามารถเปลี่ยนการแสดงผลให้ส่งออกไปยังแฟ้มข้อมูลที่ต้องการ
ทำได้โดยเพิ่มรหัสข้างล่างที่ท้ายคำสั่งเดิม
  > path_to_file.json
  
ถ้าเป็นการ append ข้อมูลเข้าไปในแฟ้มข้อมูลเดิม ให้ใช้ >> แทน >

การใช้ SQLite สำหรับเก็บและจัดการข้อมูล GeoJSON (Part 1)

SQLite: คือ ระบบจัดการฐานข้อมูล (Database Management System) ขนาดเล็กที่เก็บข้อมูลเป็นไฟล์เดียว ไม่ต้องมีเซิร์ฟเวอร์ นิยมใช้ในอุปกรณ์พกพาและระบบฝังตัว (Embedded Systems) เช่น OpenWrt, ในแอ็พลิเคชันที่ทำงานบนสมาร์ทโฟน ฯลฯ
GeoJSON: คือ รูปแบบการจัดเก็บข้อมูลทางภูมิศาสตร์ (Geographic Data) โดยใช้โครงสร้างแบบ JSON เพื่อระบุตำแหน่งบนแผนที่ เช่น จุด (Point), เส้น (Line) หรือพื้นที่ (Polygon) พร้อมข้อมูลคุณลักษณะต่างๆ

เครื่องมือและฟังก์ชันจัดการ JSON ใน SQLite

ฟังก์ชันพื้นฐาน (Essential Functions)

json(string): ตรวจสอบความถูกต้องและย่อขนาด JSON string

json_extract(json, '$path'): ดึงข้อมูลออกมาตาม Path (เช่น $.title)

json_insert / json_replace: ใช้เพิ่มหรือแทนที่ข้อมูลในชุด JSON

json_set: ใช้แก้ไขข้อมูล (ถ้ามีอยู่แล้วจะแทนที่ ถ้าไม่มีจะเพิ่มให้ใหม่)

json_remove: ใช้ลบ Key หรือสมาชิกใน Array ที่ไม่ต้องการ

ไวยากรณ์ลัด (Fast Access Syntax - SQLite 3.38+)

->: ส่งค่ากลับมาเป็น JSON string

->>: ส่งค่ากลับมาเป็น ข้อมูลดิบ (SQL value) เช่น ข้อความหรือตัวเลข

ฟังก์ชันสำหรับการรวมกลุ่ม (Aggregate Tools)

json_group_array(value): รวมข้อมูลจากหลายๆ แถว (Rows) ให้กลายเป็น JSON Array [...] (ใช้ทำ FeatureCollection)

json_group_object(name, value): รวมข้อมูลให้กลายเป็น JSON Object {"key": "value"}

ฟังก์ชันแบบตาราง (Table-Valued Functions)

json_each(json): แตกข้อมูล JSON ออกมาเป็นแถวเหมือนตาราง (เฉพาะชั้นบนสุด)

json_tree(json): แตกข้อมูล JSON ออกมาทุกชั้น (Recursive) เหมาะสำหรับการค้นหาข้อมูลที่ฝังอยู่ลึกๆ

บทความชุดนี้จะแสดงการเก็บข้อมูลใน SQLite ที่เหมาะสมกับข้อมูล GeoJSON ในสถานการณ์ต่างๆ 
ที่วิศวกรสำรวจและนักภูมิสารสนเทศอาจเคยประสบโดยแยกเป็น 2 กรณี

กรณี 1 เก็บข้อมูลแบบปกติทั่วไป (text, number, ...) แล้วใช้คำสั่งเรียกมาประกอบกันเป็น GeoJSON เมื่อต้องการ
กรณี 2 เก็บข้อมูลเป็นรูปแบบที่คล้าย GeoJSON มากที่สุด แล้วใช้ฟังก์ชันพิเศษช่วยแปลงเป็น GeoJSON

ทั้งสองกรณีจะอยู่ใน blog ถัดไปตามลำดับ
การใช้ฟังก์ชันจัดการ JSON จะมีแทรกอยู่ในขั้นตอนต่างๆ ของการทำงาน

Web Application สำหรับการปรับเปลี่ยนหลักฐานและกรอบพิกัดอ้างอิง (ตอน ๒)

การใช้สคริปต์ในการคำนวณหรือตรวจสอบด้วยโปรแกรมอื่น

(เนื้อหาในบล็อกนี้ต่อจากตอนแรก)

ผลลัพธ์จากการคำนวณส่วนที่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()

Web Application สำหรับการปรับเปลี่ยนหลักฐานและกรอบพิกัดอ้างอิง

ในการประชุมคณะอนุกรรมการบริหารจัดการโครงสร้างพื้นฐานระบบดาวเทียมนำทาง (GNSS) ของประเทศครั้งที่ ๑/๒๕๖๔
เนื้อหาสำคัญที่ขอยกมากล่าว
	...
	นาย... (อนุกรรมการและเลขานุการ) แจ้งให้ที่ประชุมทราบว่า การใช้ ITRF ๒๐๑๔
ถือเป็นโอกาสอันดีที่จะเปลี่ยนแปลงกรอบอ้างอิงค่าพิกัดของประเทศไทย เพื่อให้หน่วยงานใช้ระบบเดียวกันทั้งประเทศ 
เพื่อแก้ไขปัญหาการเป็นเอกภาพของค่าพิกัดในอดีตที่ผ่านมา และเห็นควรให้จัดสัมมนาแบบออนไลน์ 
เพื่อส่งเสริม ประชาสัมพันธ์การนำกรอบอ้างอิง ITRF ๒๐๑๔ และเวลามาตรฐาน UTC (NIMT) มาใช้ในการ-
อ้างอิงค่าพิกัดของประเทศไทย  รวมทั้งต้องเร่งดำเนินการพัฒนา Application 
สำหรับการปรับเปลี่ยนหลักฐานและกรอบพิกัดอ้างอิงเป็น ITRF ๒๐๑๔ เพื่อเตรียมความพร้อมสำหรับการใช้งาน
อย่างเป็นทางการ โดยกิจกรรมที่กล่าวถึงต้องขอความร่วมมือจากหน่วยงานที่มีสถานี CORS 
และองค์ความรู้ของผู้เชี่ยวชาญด้าน GNSS เป็นหลักในการดำเนินงานในครั้งนี้
	...
มติที่ประชุม
มติที่ประชุม ๖๑/๘ :
๑. เห็นชอบให้ศูนย์ข้อมูลค่าอ้างอิงพิกัดแบบต่อเนื่องแห่งชาติ ใช้กรอบอ้างอิง ITRF ๒๐๑๔ เป็นกรอบอ้างอิง
ค่าพิกัดของประเทศไทย และมาตรเวลา UTC (NIMT) สำหรับการอ้างอิงค่าเวลาของประเทศไทย
๒. มอบหมายให้ ฝ่ายเลขานุการฯ จัดสัมมนา เพื่อส่งเสริม ประชาสัมพันธ์การนำกรอบอ้างอิง ITRF ๒๐๑๔
และเวลามาตรฐาน UTC (NIMT) มาใช้ในการอ้างอิงค่าพิกัดของประเทศไทย
๓. มอบหมายให้ ฝ่ายเลขานุการฯ ประสานหน่วยงานที่เกี่ยวข้อง (กรมที่ดิน กรมโยธาธิการและผังเมือง
สถาบันมาตรวิทยาแห่งชาติ สถาบันสารสนเทศทรัพยากรน้ำ กรมแผนที่ทหาร และ สทอภ.) เตรียมความพร้อม
ในการประกาศใช้กรอบอ้างอิง ITRF ๒๐๑๔ เป็นกรอบอ้างอิงค่าพิกัดที่เป็นเอกภาพของประเทศไทยต่อไป
๔. มอบหมายให้หน่วยงานต่างๆ ที่มี CORS ใน NCDC ประสานงานกับสถาบันมาตรวิทยาแห่งชาติ
เพื่อดำเนินการให้เกิดระบบและกระบวนการเทียบเวลามาตรฐาน UTC (NIMT) เพื่อการบ่งชี้คุณลักษณะ ด้านเวลา
ของนาฬิกาภายในเครื่องรับสัญญาณ GNSS ของ NCDC และการอ้างอิงไปสู่เวลามาตรฐานระหว่างประเทศ

Geo-4d Coordinates Transformation Application ที่ผู้เขียนทำและเผยแพร่แล้ว
มติของคณะอนุกรรมการดังกล่าวเป็นสิ่งที่ถูกต้องเหมาะสมแล้ว บางเรื่องสามารถดำเนินการได้ทันที 
ผู้เขียนสงสัยว่าทำไมไม่มอบหมายไปเลยว่าให้หน่วยงานใดเป็นหน่วยงานหลักในการทำ Application
สำหรับตัวผู้เขียนได้ตระหนักถึงความสำคัญนี้มานานแล้ว และได้ทำ Web Application 
สำหรับงานแปลงค่าพิกัดระหว่างกรอบอ้างอิงออกเผยแพร่ในกลุ่มผู้ใกล้ชิดในปี พ.ศ.2564 
ผู้เขียนให้ชื่อผลงานนี้ว่า Geo-4d Coordinates Transformation

ขีดความสามารถของ Geo-4d Coordinates Transformation (กำหนดชื่อย่อเป็น G4DX)
G4DX มีขีดความสามารถดังนี้:-
แปลงค่าพิกัดระหว่างกรอบอ้างอิง
สร้างสคริปต์สำหรับใช้ในการคำนวณหรือตรวจสอบด้วยโปรแกรมอื่น
ตัวอย่างการใช้งาน: การแปลงค่าพิกัดจากกรอบอ้างอิง ITRF2000 ไปยัง ITRF2014
คำอธิบายการคำนวณในรูปข้างล่าง

ข้อมูลป้อนเข้า
แปลงค่าพิกัดจากกรอบอ้างอิง ITRF2000 ไปยัง ITRF2014
พิกัดตำแหน่ง ละติจูด = 14°N, ลองจิจูด = 100°E
ความสูงเหนือสัณฐานทรงรีอ้างอิง = 0 m.
Epoch-of-observation = 2004.5 (ขณะเวลากลางปีคศ.2004)
ผลลัพธ์จากการคำนวณ
พิกัดตำแหน่ง ละติจูด = 14.00000014°N, ลองจิจูด = 100.00000000°E
ความสูงเหนือสัณฐานทรงรีอ้างอิง = -0.0065 m.
Epoch-of-observation = 2004.5 (carried on value)
Reference frame = ITRF2014

ข้อสังเกต 
ค่าของขณะเวลาสังเกตเป็นค่าปีคศ.เติมทศนิยม 1 ตำแหน่ง และไม่มีประโยชน์ที่จะใช้เกินกว่า 1 คำแหน่ง
ค่าขณะเวลาสังเกตมีบทบาทสำคัญในการคำนวณหาค่าพารามิเตอร์ในการแปลง 7 ตัวที่ใช้ประกอบในกระบวนวิธี
7-parameter Helmert Transformation ค่าขณะเวลาสังเกตที่ละเอียดเพียง 1 ตำแหน่งเพียงพอแล้วสำหรับใช้งาน
อนึ่ง ถ้าค่ามีทศนิยมเป็นศูนย์นิยมเขียนให้มีค่าทศนิยมด้วยอย่างสมบูรณ์ เช่น 2022.0

เกี่ยวกับ G4DX
- ทำงานแบบ client-server ทำการประมวลผลที่ server
- การทำงานที่เซิร์ฟเวอร์โดยใช้ Python + pyproj module
- ประมวลผลได้ครั้งละ 1 จุด ตามข้อจำกัดของ user interface
- ผลลัพธ์ที่เป็นข้อความในกรอบ textarea มีเนื้อหาและสคริปต์ที่มีประโยชน์ต่อการศึกษาและใช้งาน
- สคริปต์ดังกล่าวสามารถใช้งานได้โดย copy/paste/run นำไปดัดแปลงให้ตรงความต้องการผู้ใช้ได้ง่าย

(ตอนต่อไป: การใช้สคริปต์ที่เป็นผลลัพธ์ในการคำนวณหรือตรวจสอบด้วยโปรแกรมอื่น)

ได้ทดลองใช้ Codepen เขียนโค้ดทำแผนที่ไว้มากมายหลายแบบ แบบที่แสดงไว้ข้างล่างมีจุดเด่นคือมี basemap ของ OpenStreetMap บวกกับ อาคารสามมิติ ของ OSM Buildings ทำให้ได้มิติที่สามเพิ่มจากแผนที่แบบเดิมๆ ที่พิเศษสุดคือได้ลองเขียน 3D model ง่ายๆ ซึ่งก็คือเสาธงจุฬาฯครับ อนุญาตให้ fork เอาไปใช้งานต่อได้ครับ