คอมไพล์ Python Script เป็นไฟล์ Executable ด้วย PyInstaller

PyInstaller

คือเครื่องมือที่ช่วยการแปลงโปรแกรมที่เขียนด้วยไพทอนเป็น execute binary file  ที่สามารถนำไปรันได้โดยที่เครื่องคอมพิวเตอร์ปลายทางไม่ต้องติดตั้งไพทอน สำหรับ PyInstaller เป็น cross-platform สามารถใช้งานได้บนวินโดส์ แมค และลีนุกซ์ สนับสนุนไพทอนรุ่น 2.7 และ ไพทอน รุ่น 3.3 ถึง 3.6 จุดมุ่งหมายของ PyInstaller คือต้องการช่วยผู้ใช้ในการแปลงโปรแกรมไพทอน ที่ใช้โมดูลไลบรารีภายนอกเช่น Matplotlib, DJango, wxPython, PyQt เป็นต้น ให้สามารถทำได้ง่ายสะดวก

ติดตั้ง PyInstaller

ติดตั้งง่ายๆด้วยคำสั่ง pip ใน command prompt

pip install pyinstaller

ใช้งาน PyInstaller

การใช้งานสามารถใช้งานผ่าน command line ได้ แต่สำหรับโปรแกรมที่เรียกใช้โมดูลไลบรารีข้างนอกและต้องขนข้อมูล (data) ที่โมดูลไลบรารีนั้นๆต้องการใช้  ผมแนะนำให้ใช้ไฟล์สคริปท์ (Spec file) มาช่วยจะดีกว่า ปรับแต่งได้มากกว่า ตัว spec file จริงๆก็คือไฟล์สคริปท์ของไพทอนนั่นเอง กรณีที่ต้องใช้ Spec file อีกกรณีหนึ่งคือต้องการขนรันไทม์ไลบรารีเช่น .dll หรือ .so ไปแบบแมนวล กรณีที่ผมเจอคือผมใช้ PySide2 ที่รุ่นทางการจริงๆยังไม่ออกมา แต่ hook file ก็มีมาให้แล้วพร้อมกับ PyInstaller รุ่นใหม่ 3.3 แต่ผมใช้งานแล้วยังไม่สำเร็จ ดังนั้นจึงต้องใช้ Spec file นี้เป็นตัวช่วยในการขนรันไทม์ไลบรารีไป ส่วนเรื่อง hook file คืออะไรค่อยว่าอีกที

กรณีศึกษาด้วย Surveyor Pocket Tools บนวินโดส์

โปรแกรม Surveyor Pocket Tools พัฒนาด้วยไพทอน ปัจจุบันใช้ไพทอน รุ่น 3.6 ใช้โมดูลไลบรารีข้างนอกคือ openpyxl, pyproj, geographiclib, gmplot, simplekml, pyshp และที่ขาดไม่ได้คือ PySide2 ซึ่งสำหรับ openpyxl และ pyproj จะมีการขนข้อมูลไปด้วย ส่วน PySide2 ผมจะขนไฟล์ dll  ที่ต้องการด้วยมือล้วนๆ

Spec file ของ Surveyor Pocket Tools

มาดูไฟล์สคริปนี้ ผมตั้งชื่อว่า “setup.spec”

# -*- mode: python -*-
# -*- mode: python -*-
import sys
import PySide2
import os
block_cipher = None
 
dirname = os.path.dirname(PySide2.__file__)
plugins_path = os.path.join(dirname, 'plugins', '')
 
pyside2_plugins = [(plugins_path + 'iconengines/*', 'plugins/iconengines/'),
                  (plugins_path + 'imageformats/*', 'plugins/imageformats/'),
		  (plugins_path + 'platforms/*', 'plugins/platforms/'),
		  (plugins_path + 'printsupport/*', 'plugins/printsupport/'),
		  (plugins_path + 'sqldrivers/*', 'plugins/sqldrivers/')]
 
added_files = [('markers/*', 'markers/'),
               ('geoids/*', 'geoids/'),
	       ('database/*', 'database/'),
	       ('example data/*', 'example data/'),
	       ('qt.conf', ''), 
	       ('*.xml', '')]
 
a = Analysis(['main.py'],
             pathex=['D:\\sourcecodes\\python\\surveyor pocket tools'],
             binaries=None,
             datas=added_files + pyside2_plugins,
             hiddenimports=[],
             hookspath=[],
             runtime_hooks=[],
             excludes=[],
             win_no_prefer_redirects=False,
             win_private_assemblies=False,
             cipher=block_cipher)
pyz = PYZ(a.pure, a.zipped_data,
             cipher=block_cipher)
exe = EXE(pyz,
          a.scripts,
          exclude_binaries=True,
          name='surveyor pocket tools',
		  icon='Land Surveying-64.ico',
          debug=False,
          strip=False,
          upx=False,
          console=False )
coll = COLLECT(exe,
               a.binaries,
               a.zipfiles,
               a.datas,
               strip=False,
               upx=True,
               name='setup')

ลองมาดูโค้ดกัน เริ่มจาก import PySide2 เข้ามาเพื่อจะตรวจสอบว่า PySide2 ที่เราใช้งานเป็น 32 บิตหรือ 64 บิต เพื่อจะได้ขน .dll ไปถูกรุ่น จากนั้นเก็บไดเรคทอรีของ PySide2 เข้าเก็บใน dirname ผ่านฟังก์ชัน os.path.dirname() ที่นี้เราทราบว่าในไดเรคทอรีของ PySide2 จะมีไดเรคทอรีย่อยชื่อ “plugins” อยู่ ทำการเก็บไดเรคทอรีนี้ด้วยฟังก์ชั่น os.path.join() ไปเก็บไว้ในตัวแปร plugins_path

# -*- mode: python -*-
import sys
import PySide2
import os
block_cipher = None
 
dirname = os.path.dirname(PySide2.__file__)
plugins_path = os.path.join(dirname, 'plugins', '')

ต่อไปคือตัวแปร pyside2_plugins จะเป็นลิสต์เก็บ tuple โดยสมาชิกตัวแรกจะเก็บชื่อไฟล์ไดเรคทอรีต้นทาง ใช้เครื่องหมาย * เพราะต้องการทุกๆไฟล์ในไดเรคทอรีนี้ สมาชิกตัวที่สอง จะเก็บชื่อไดเรคทอรีปลายทางที่ต้องการไฟล์เหล่านี้ไปอยู่

 
pyside2_plugins = [(plugins_path + 'iconengines/*', 'plugins/iconengines/'),
                   (plugins_path + 'imageformats/*', 'plugins/imageformats/'),
		   (plugins_path + 'platforms/*', 'plugins/platforms/'),
		   (plugins_path + 'printsupport/*', 'plugins/printsupport/'),
		   (plugins_path + 'sqldrivers/*', 'plugins/sqldrivers/')]

ลองมาดูว่าไดเรคทอรี “plugins” ผมไฮไลท์ไว้เฉพาะไดเรคทอรีที่โปรแกรม Surveyor Pocket Tools ต้องการ

ต่อไปจะขนไฟล์ที่โปรแกรม Surveyor Pocket Tools ต้องการใช้ ให้ใส่ไว้ที่ตัวแปร added_files โครงสร้างเป็น tuple เหมือนกัน และขนไฟล์ชื่อ qt.conf ที่ PySide2 ต้องการไปด้วย

added_files = [('markers/*', 'markers/'),
               ('geoids/*', 'geoids/'),
	       ('database/*', 'database/'),
	       ('example data/*', 'example data/'),
	       ('qt.conf', ''), 
	       ('*.xml', '')]

มาดูไดเรคทอรีที่โปรแกรมต้องการดังนี้

ต่อไปมาดูโค้ดส่วนที่สำคัญมาก ‘main.py’ คือไฟล์สคริปท์หลักของโปรแกรม Surveyor Pocket Tools ต่อไปคือ pathex เป็นไดเรคทอรีของไฟล์ไพทอนสคริปท์ และ datas ที่ผมจัดการรวม added_files และ pyside2_plugins เข้าด้วยกัน สุดท้าย hookspath คือไดเรคทอรีที่เก็บไฟล์ hook ไว้ สำหรับไฟล์ hook นี้ PyInstaller จะอ่านสคริปท์นี้ทีละไฟล์มาตัดสินใจว่าจะขนข้อมูลไดเรคทอรีไหนไป ผมเลือกใช้ดีฟอลท์ครับคือปล่อยว่าง

a = Analysis(['main.py'],
             pathex=['D:\\sourcecodes\\python\\surveyor pocket tools'],
             binaries=None,
             datas=added_files + pyside2_plugins,
             hiddenimports=[],
             hookspath=[],
             runtime_hooks=[],
             excludes=[],
             win_no_prefer_redirects=False,
             win_private_assemblies=False,
             cipher=block_cipher)

สำหรับไดเรอทอรี hooks ที่เป็นดีฟอลท์มากับ PyInstaller ผมใช้ไฟล์เพียงสองไฟล์เท่านั้น ตามที่ไฮไลท์ไว้

ใช้ PyInstaller คอมไพล์ไฟล์ setup.spec

ผมใช้ Minoconda เมื่อจะคอมไพล์ก็เรียก command prompt มาดังนี้ ใช้คำสั่ง cd เข้ามาที่พาทของสคริปท์ของไพทอน ใช้คำสั่ง dir ดูไฟล์ setup.spec

ต่อไปทำการคอมไพล์ ด้วยคำสั่ง

pyinstaller setup.spec

ผลลัพธ์ของ PyInstaller

เมื่อคอมไพล์เสร็จแล้ว ไม่มี error จะได้ไดเรคทอรีมาสองคือ “build” และ “dist” เมื่อเข้าไปดูใน “dist” จะเห็นไดเรคทอรีย่อยช “setup” ชื่อไดเรคทอรีนี้ PyInstaller จะสร้างตามชื่อหน้าของไฟล์ setup.spec เมื่อเข้าดูที่ไดเรคทอรี “setup” จะเห็นไฟล์ต่างๆที่โปรแกรมต้องการ

ผมลองดับเบิ้ลคลิกไฟล์ “surveyor pocket tools.exe” ก็สามารถเปิดมาและทำงานได้ตามปกติ ลองดูชื่อไดเรคทอรีจะเห็นสองไดเรคทอรี ที่ได้จากไฟล์ hooks คือ openpyxl และ pyproj ลองเข้าไปดูในไดเรคทอรี จะเห็นข้อมูลที่ pyproj ขนไปใช้ หมายเหตุว่าข้อมูลนี้ pyproj จะนำไปเป็นฐานข้อมูลในการแปลงพิกัดตาม datum และ projection

ทำไฟล์ Setup ด้วย Inno Setup

จากนั้นผมจะ copy ไดเรคทอรีที่อยู่ใน “setup” ไปไว้อีกที่หนึ่ง พื้นที่นี้สำหรับใช้ Inno Setup มาทำไฟล์ติดตั้ง ลองดูไดเรคทอรี

ในไดเรคทอรีนี้ผมจะมีไฟล์ “surveyorpockettools64.iss” เป็นไฟล์สคริปท์ของ Inno Setup เพื่อสร้างไฟล์ติดตั้ง setup สำหรับวินโดส์ 64 บิต

#define MyAppName "Surveyor Pocket Tools"
#define MyAppEXE "Surveyor Pocket Tools.exe"
#define MyShortAppName "SurveyorPocketTools"
#define MyMainRoot "Survey Suite"
#define Developer "Prajuab Riabroy"
#define Version "0.98"
#define Build "573"
 
[Setup]
AppName={#MyAppName}
AppVerName={#MyAppName} V{#Version}
DefaultDirName={pf}\{#MyMainRoot}\{#MyAppName}
DefaultGroupName={#MyMainRoot}\{#MyAppName}
UseSetupLdr=yes
UninstallDisplayIcon={app}\{#MyAppEXE}
VersionInfoProductName={#MyAppName}
VersionInfoCompany=priabroy
VersionInfoCopyright=Copyright 2000-2017 by {#Developer}
VersionInfoDescription={#MyAppName}
VersionInfoProductVersion={#Version}
VersionInfoVersion={#Version}
OutputDir=Setup
OutputBaseFilename={#MyShortAppName}V{#Version}Build{#Build}Setup64
;OutputDir=TraverseProV250Setup64
; "ArchitecturesAllowed=x64" specifies that Setup cannot run on
; anything but x64.
ArchitecturesAllowed=x64
; "ArchitecturesInstallIn64BitMode=x64" requests that the install be
; done in "64-bit mode" on x64, meaning it should use the native
; 64-bit Program Files directory and the 64-bit view of the registry.
;ArchitecturesInstallIn64BitMode=x64
ArchitecturesInstallIn64BitMode=x64
AppPublisher={#Developer}
AppPublisherURL=https://www.surveyorpockettools.org
AppVersion={#Version}.{#Build}
LicenseFile = eula.txt
ChangesEnvironment=yes
SolidCompression=yes
Compression=lzma2/ultra64
LZMAUseSeparateProcess=yes
LZMADictionarySize=1048576
LZMANumFastBytes=273
 
[Files]
Source: "{#MyAppName}.exe"; DestDir: "{app}"
Source: "base_library.zip"; DestDir: "{app}" ;
Source: "plugins\*"; DestDir: "{app}\plugins\"; Flags: ignoreversion recursesubdirs
Source: "database\*"; DestDir: "{userappdata}\{#MyAppName}\database\";
Source: "geoids\*"; DestDir: "{userappdata}\{#MyAppName}\geoids\";
Source: "markers\*"; DestDir: "{userappdata}\{#MyAppName}\markers\";
Source: "example data\*"; DestDir:"{userappdata}\{#MyAppName}\example data\";
Source: "pyproj\data\*"; DestDir: "{app}\pyproj\data\";
Source: "openpyxl\*"; DestDir: "{app}\openpyxl\";
;Source: "requests\*"; DestDir: "{app}\requests\";
Source: "*.html"; DestDir: "{userappdata}\{#MyAppName}\";
Source: "*.dll"; DestDir: "{app}"
Source: "*.pyd"; DestDir: "{app}"
Source: "*.xml"; DestDir: "{userappdata}\{#MyAppName}\";
Source: "qt.conf"; DestDir: "{app}"
 
[Icons]
;create icon at start menu group
Name: "{group}\{#MyAppName}"; Filename: "{app}\{#MyAppExe}"
;create icon at desktop
Name: "{commondesktop}\{#MyAppName}"; FileName:"{app}\{#MyAppExe}"
 
[Registry]
; Start "Software\My Company\My Program" keys under HKEY_CURRENT_USER
; and HKEY_LOCAL_MACHINE. The flags tell it to always delete the
; "My Program" keys upon uninstall, and delete the "My Company" keys
; if there is nothing left in them.
Root: HKCU; Subkey: "Software\{#MyMainRoot}"; Flags: uninsdeletekeyifempty
Root: HKCU; Subkey: "Software\{#MyMainRoot}\{#MyAppName}"; Flags: uninsdeletekey
Root: HKLM; Subkey: "Software\{#MyMainRoot}"; Flags: uninsdeletekeyifempty
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}"; Flags: uninsdeletekey
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}\Settings"; ValueType: string; ValueName: "InstalledPath"; ValueData: "{app}"
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}\Settings"; ValueType: string; ValueName: "DevelopedBy"; ValueData: "{#Developer}"
Root: HKLM; Subkey: "Software\{#MyMainRoot}\{#MyAppName}\Settings"; ValueType: string; ValueName: "ApplicationName"; ValueData: "{#MyAppName}"
;Root: HKCU; Subkey: "Environment"; ValueType:string; ValueName:"PROJ_LIB"; ValueData:"{userappdata}\{#MyAppName}\geoidgrids\" ; Flags: preservestringtype ;

ถ้าเป็นไฟล์สำหรับ Surveyor Pocket Tools รุ่น 32 บิตเพียงใส่คอมเมนต์หน้า ;ArchitecturesInstallIn64BitMode=x64 ก็พอ สำหรับรายละเอียดสคริปท์ของ Inno Setup ผมจะไม่กล่าวถึงรายละเอียดในที่นี้ผู้อ่านที่สนใจสามารถศึกษาได้ครับ จากนั้นก็ใช้ Inno Setup ทำการ build ก็จะได้ไฟล์ Exe เดี่ยวๆ ที่สามารถ zip ไปให้ผู้ใช้ได้ download ต่อไป พบกันตอนหน้าครับ

สนุกกับโปรแกรมเครื่องคิดเลขสำหรับงานสำรวจ ตอนที่ 1 โปรแกรมแปลงพิกัด “Geo2UTM” บนเครื่องคิดเลข Casio FX 5800P

  • สวัสดีครับท่านผู้อ่านทุกท่าน ในโอกาสที่บ.เคเอ็นเอส เอ็นจิเนียริ่ง คอร์ปอเรชั่น จำกัด มีอายุใกล้จะขวบปีแล้ว ผมในฐานะวัยแล้วเกือบจะรุ่นพ่อของน้องๆชุดนี้แล้ว เห็นความตั้งใจของน้องๆ และก็ยินดีเป็นอย่างยิ่ง ในโอกาสนี้ด้วยความที่สนิทสนมกันก็ถูกลากมาให้เขียนบทความให้เพื่อฉลองครบรอบหนึ่งปี และลงที่นี่ ไม่ใช่ที่บล็อก priabroy.com ที่ประจำ ตอนแรกคิดๆอยู่ว่าจะเขียนเรื่องอะไรดี แต่สุดท้ายก็จะเขียนเรื่องโปรแกรมมิ่งเครื่องคิดเลข เพราะว่าเครื่องคิดเลขเป็นเครื่องมือพื้นฐานที่สุดสำหรับช่างสำรวจ ช่างโยธาของเรา
  • บทความนี้คงจะมีหลายตอนก็มาติดตามกัน

ย้อนอดึตแห่งความทรงจำ

  • เครื่องคิดเลข Casio นับว่าเป็นขวัญใจของเราช่างสำรวจ ช่างโยธา ตอนเรียนอยู่วิทยาลัยหรือมหาวิทยาลัย บางคนอาจจะมีประสบการณ์ในการพาเครื่องไปฝากไว้ที่โรงจำนำเพราะกลัวทำหาย ^-^ ตอนสมัยเรียนช่วงปี 30-33 ผมใช้ FX3800 ปัจจุบันเลิกผลิตไปนานแล้ว เหลือไว้แต่ความทรงจำ สมัยแต่ก่อนเขียนโปรแกรมลงเครื่องพวกนี้คงได้แค่โปรแกรมเล็กๆ เพราะเครื่องมื memory ที่จำกัดจำเขี่ยมากๆ
Casio Fx3800
Casio Fx3800
  • ต่อมาทำงานได้สักสองปีก็ไปถอยเอาเครื่องในตำนาน  FX-880P ได้มาเครื่องหนึ่ง รุ่นนี้ปัจจุบันขึ้นหิ้งเป็นตำนานไปแล้ว ทั้งที่ผ่านระยะเวลามายี่สิบกว่าปี ยังมีคนตามล่าหากัน รุ่นนี้เวลาพกพาก็เสียบไว้ที่กระเป๋าหลังของกางเกงยีนส์ ดูมันเท่ห์ แต่พกแบบนี้ บางคนเผลอนั่งทับจนเครื่องคิดเลขหักเป็นท่อน น้ำตาตกกันมาแล้วก็มี โปรแกรมมิ่งรุ่นนี้ใช้โปรแกรมภาษาเบสิคแบบมีหมายเลขบรรทัดกำกับ เปิดโลกโปรแกรมมิ่งไปอีกหลายระดับ เขียนโปรแกรมยากๆได้พอสมควร การตั้งตัวแปรใช้ตัวอักษรหลายตัวได้ เครื่องเดิมๆ มี memory 32 กิโลไบต์ มันเยอะพอสมควร
Casio FX - 880P
Casio FX – 880P
  • ด้านหลังเครื่องรุ่นนี้ยังมีช่องให้ใส่แรมเพิ่มได้อีก 32 กิโลไบต์ ผมอุตส่าห์ไปเดินแถวสะพานเหล็ก คลองถมจนได้มาหนึ่งอัน เมื่อใส่แล้วก็รวมกันได้ 64 กิโลไบต์ คิดเป็น 65536 ไบต์ เหลือเฟือ ขนาดเขียนวงรอบเล่นๆ ยังไม่เต็มเลย
s-l300
Memory Pack 32KB for FX-880P
  • มาดูโปรแกรมสำหรับเครื่องคิดเลข FX-880P เพื่อรำลึกความหลัง ผมเขียนโปรแกรมนี้หาจุดตัดระหว่างเส้นตรงกับเส้นตรง เส้นตรงกับวงกลม วงกลมกับวงกลม ผมลงมาให้เต็มๆแบบไม่มีตัดทอน ตอนนี้ผมไม่มีเครื่องคิดเลขรุ่นนี้ให้ลองแล้ว

5 ‘Intersection
10 CLS:BEEP:BEEP1:ANGLE 0:SET F5:PRINT CHR$(9);:CLS
20 PRINT ” ***Find Intersection point***”
30 CLS:PRINT “<<1:AZI#AZI>> <<2:AZI#DIST>>”;:PRINT
40 PRINT”<<3:DIST#DIST>> <>”;
50 T$=INKEY$
60 IF (T$=”1″) THEN 150
70 IF (T$=”2″) THEN 600
80 IF (T$=”3″) THEN 800
90 IF (T$=”4″) THEN 400
100 IF (T$=”Q”) THEN PRINT TAB(10);”<<>>”;
:SET F9:END
110 GOTO 50
150 CLS:PRINT “N1=”;N1;:INPUT N1
160 PRINT “E1=”;E1;:INPUT E1
170 PRINT “AZIMUTH1=”;AZI1;:INPUT AZI1
180 FANG=AZI1:GOSUB 3000:DAZI1=DANG
190 PRINT “N2=”;N2;:INPUT N2
200 PRINT “E2=”;E2;:INPUT E2
210 PRINT “AZIMUTH2=”;AZI2;:INPUT AZI2
220 FANG=AZI2:GOSUB 3000:DAZI2=DANG
230 DELTY=N2-N1:DELTX=E2-E1
240 GOSUB 3500:AZI12=Y:DIST12=X
250 NI=(TAN(DAZI2)*N2-TAN(DAZI1)*N1+E1-E2)/(TAN(DAZI2)-
TAN(DAZI1))
260 EI=(NI-N1)*TAN(DAZI1)+E1
290 PRINT “NI= “;NI;:PRINT
300 PRINT “EI= “;EI
310 GOTO 30
400 ‘Four points
410 CLS:PRINT “N1= “;N1;:INPUT N1
420 PRINT “E1= “;E1;:INPUT E1
430 PRINT “N2= “;N2;:INPUT N2
440 PRINT “E2= “;E2;:INPUT E2
450 PRINT “N3= “;N3;:INPUT N3
460 PRINT “E3= “;E3;:INPUT E3
470 PRINT “N4= “;N4;:INPUT N4
480 PRINT “E4= “;E4;:INPUT E4
490 DELTY=N2-N1:DELTX=E2-E1:GOSUB 3500:AZI12=Y:DIST12=X
500 DELTY=N4-N3:DELTX=E4-E3:GOSUB 3500:AZI34=Y:DIST34=X
510 DELTY=N3-N1:DELTX=E3-E1:GOSUB 3500:AZI=Y
512 IF (AZI12=90 OR AZI12=270) THEN NI=N1 ELSE 515
513 EI=(NI-N3)*TAN(AZI34)+E3
515 IF (AZI34=90 OR AZI34=270) THEN NI=N3 ELSE 520
520 IF NOT((AZI12=90 OR AZI12=270) OR (AZI34=90 OR
AZI34=270)) THEN NI=(TAN(AZI34)*N3-TAN(AZI12)*N1+E1-
E3)/(TAN(AZI34)-TAN(AZI12)) ELSE 540
530 EI=(NI-N1)*TAN(AZI12)+E1
540 PRINT “NI= “;NI;:PRINT
550 PRINT “EI= “;EI
560 GOTO 30
600 ‘AZI # DIST
610 CLS:PRINT “N1= “;N1;:INPUT N1
620 PRINT “E1= “;E1;:INPUT E1
630 PRINT “AZIMUTH= “;INAZI;:INPUT INAZI:FANG=INAZI:GOSUB
3000:DAZI1=DANG
635 CANG=0:PAZI=DAZI1:GOSUB 4500:DAZI2=NAZI
640 PRINT “N2= “;N2;:INPUT N2
645 PRINT “E2= “;E2;:INPUT E2
650 PRINT “DIST= “;DIST;:INPUT DIST
660 DELTY=N2-N1:DELTX=E2-E1:GOSUB 3500:DIST12=X:AZI12=Y
670 PHI1=AZI12-DAZI1:IF PHI1 675 PHI2=AZI12-DAZI2:IF PHI2 Int No.1
900 CANG=ANG2:PAZI=AZI12:GOSUB 4500:AZI2I2=NAZI
910 PNI=N2+DIST2*COS(AZI2I1)
920 PEI=E2+DIST2*SIN(AZI2I1)
930 MNI=N2+DIST2*COS(AZI2I2)
940 MEI=E2+DIST2*SIN(AZI2I2)
950 CLS:BEEP:PRINT “NI(1)= “;PNI;:PRINT
960 PRINT “EI(1)= “;PEI
970 CLS:PRINT “NI(2)= “;MNI;:PRINT
980 PRINT “EI(2)= “;MEI
990 GOTO 30
3000 ‘Convert input angle to degree
3010 DD=FIX(FANG)
3020 TEMP=FRAC(FANG)*100
3030 MM=FIX(TEMP)
3040 SS=FRAC(TEMP)*100
3050 DANG=DEG(DD,MM,SS)
3060 RETURN
3500 ‘Find Azimuth
3510 X=POL(DELTY,DELTX)
3520 IF Y180 THEN 4550 ELSE 4530
4530 TEMP=TEMP+180
4540 GOTO 4590
4550 IF TEMP>540 THEN 4580
4560 TEMP=TEMP-180
4570 GOTO 4590
4580 TEMP=TEMP-540
4590 NAZI=TEMP
4600 RETURN

ถึงยามต้องพรากจากกัน

  • ผ่านไปหลายปีสำหรับการทำงานในภาคสนาม ไปไหนมาไหนก็หอบเครื่องคิดเลขรุ่นนี้ไปใช้งาน ก่อนจะจากกันผมทิ้งเครื่องคิดเลขไว้ที่ตู้คอนเทนเนอร์หน้าไซต์งาน เพราะลูกน้องเอาไปใช้บ้าง ไม่ได้เอากลับที่พักตอนยามค่ำคืน คืนนั้นไฟฟ้าลัดวงจรที่ตู้คอนเทนเนอร์ ทราบข่าวอีกทีก็วอดแล้ว ยืนคอตกไปพักใหญ่ๆ เพราะว่าทรัพย์สินทางปัญญาหายวับไปกับพระเพลิง แต่ด้วยความไม่ประมาท ผมลอก Source code ไว้ที่เครื่องคอมพิวเตอร์ไว้บ้างแต่ไม่ทั้งหมด
  • ถึงแม้จะรู้ว่าใดๆในโลกนี้ล้วนอนิจจัง แต่ประสบการณ์เรื่องนี้ทำให้การใช้งานคอมพิวเตอร์ผมจะแบ็คอัพข้อมูลที่ทำงานได้สอง สามที่เสมอ และที่สำคัญคือผมชอบโปรแกรมมิ่ง จะเก็บ source code ไว้อย่างดี เก็บไว้หลายๆที่เช่นกัน แต่อย่างไรก็ตามต้องขยัน back up ครับ ถ้ามันพังในวันนี้ยังเอาของเมื่อวานมาใช้งานได้

สิบปีที่เครื่องคิดเลขห่างหาย

  • ตั้งแต่นั้นเป็นต้นมาผมก็ไม่เคยซื้อเครื่องคิดเลขอีกเลยประมาณสิบปี เพราะว่าเป็นยุคเวลาของโน๊ตบุ๊ค เอะอะจะคำนวณอะไรๆก็ต้องที่โน๊ตบุ๊ค แต่ชีวิตเหมือนขาดอะไรไป จนมาถึงยุคมือถือจอสัมผัสยิ่งแล้วไปกันใหญ่ เพราะมีโปรแกรมจำลองเครื่องคิดเลขมาให้ใช้งาน วันนั้นไปเดินห้างเห็นเครื่องคิดเลข FX 5800P วางอยู่ที่ชั้นขายของ ป้ายบอกลดราคาเหลือ 1990 บาท จากราคาเดิม 2890 บาทคิดอยู่ในใจว่ามันลดราคากระหน่ำแท้ๆ ราคานี้ไม่รวมสายลิ๊งค์ ดูสเป็คบอกว่าเขียนโปรแกรมได้ มีเม็มโมรีมาให้ 28500 ไบต์ ผมนึกถึง Fx-880P ทันที ตัดสินใจซื้อมาลอง ที่ไหนได้มาถึงบ้านเปิดดูในเน็ตเห็นราคาขายออนไลน์ราคา 1790 บาทได้สายลิ๊งค์ด้วย มันถูกว่ากว่าที่ผมซื้อหลายร้อยบาท เอาละวะ ภูมิใจที่ได้ใช้ของแพงกว่า
Casio FX5800P
Casio FX5800P

FX 5800P กับอารมณ์กระชากใจเหมือนกลับไปเรียน (Back to old school) อีกครั้ง

  • เอาละครับ ได้ใช้เครื่องคิดเลขจริงๆสักที คือหลายสิบปีที่ผ่านมาเหมือนรออะไรบางอย่าง แต่ไม่รู้ว่ามันคืออะไร เจออีกทีใช่เลย เวลากดคีย์เครื่องคิดเลขมันมีการตอบสนองได้อารมณ์เหมือนได้กลับไปใช้และเรียนมหาวิทยาลัยอีกครั้ง ผมเริ่มอ่านคู่มือเพื่อจะเขียนโปรแกรม ใช้เวลาไม่มากนักเพราะคุ้นๆอยู่ คู่มืออะไรก็หาง่ายในยุคนี้ ดาวน์โหลดมาอ่านได้สบายๆ
  • ช่วงเริ่มต้นกับเครื่องคิดเลข ผมเริ่มโปรแกรมง่ายๆก่อนเช่นจำพวก Cogo เช่นหามุม ระยะทางเมื่อกำหนดค่าพิกัดของจุดสามจุด ความรู้สึกแรกคือชอบเครื่องคิดเลขรุ่นนี้พอสมควร
  • เคยคิดเล่นๆว่าเครื่องคิดเลขพวกนี้ ยกเว้น FX-880P ที่เทพไปแล้ว จะสามารถเขียนโปรแกรมระดับ Advance ได้ไหมเช่นแปลงพิกัดไปมาระหว่าง UTM และ ค่าพิกัดภุมิศาตร์ (Geographic) หาระยะทางระหว่างสองจุดบน Ellipsoid หรือหาระยะทาง Geodesic distance

ข้อจำกัดด้านโปรแกรมมิ่งแต่ฟ้าปิดกั้นดินไม่ได้

  • จั่วหัวให้เว่อร์ซะยังงั้น ปัญหาจริงๆที่คนจะเขียนโปรแกรมพวกนี้คืออย่างแรกคือสูตร ไม่รู้จะใช้เวอร์ชั่นไหนดี อย่างที่สองคือเครื่องคิดเลขรุ่นพวกนี้มีตัวแปรจำกัด จาก A ถึง Z นับได้ 26 ตัว น้อยซะจริงๆ ถ้าสูตรมีการใช้ตัวแปรมากกว่านี้จะทำอย่างไร ปัญหานี้ยังมีทางออก เครื่องคิดเลขรุ่นนี้เตรียมตัวแปรอนุกรมให้คือ Z เราสามารถใช้งาน Z[1],Z[2],Z[3],…. ได้มากเท่าที่เมมโมรีเครื่องคิดเลขยังเหลือพอ
  • ข้อจำกัดอีกอย่างคือไม่สามารถนำไฟล์ข้อมูลให้โปรแกรมได้ ดังนั้นโปรแกรมเครื่องคิดเลขจึงจะต้องไม่ซับซ้อนมาก ถ้ามากกว่านี้ต้องใช้เครื่องคอมพิวเตอร์จะดีที่สุด

โปรแกรมแปลงพิกัดจากค่าพิกัดภูมิศาสตร์ ไปยัง ค่าพิกัดระบบพิกัดฉาก UTM

  • ก่อนจะไปต่อเรื่องโปรแกรมมิ่ง มาเรียกน้ำย่อยกันก่อน มาดูรูปการคำนวนกันก่อน ต้องการแปลงค่าพิกัด lat=14°27’44.71″ long=100°58’27.02″ ไปยังค่าพิกัด UTM

input_geo

  • แปลงพิกัดเป็น UTM ได้ค่า N=1599784.382 E=712796.211

output_utm1

  • และคำนวน zone  ของ UTM มาให้ด้วยจุดพิกัดนี้อยู่ในโซน 47

output_utm2

แปลงพิกัดบน WGS84

  • สำหรับโปรแกรมจะแปลงพิกัดบนพื้นฐาน WGS84 เท่านั้น ไม่มีการแปลงข้ามพื้นหลักฐานเพราะมันจะซับซ้อนยุ่งยากเกินกว่าเครื่องคิดเลข เมมโมรีน้อยๆจะทำได้

โปรแกรมคำนวณ

  • ข้อจำกัดอีกอย่างของเครื่องคิดเลขรุ่นนี้คือสามาถใช้สายลิ๊งค์โอนโปรแกรมจากเครื่องไปหาเครื่องอื่นเท่านั้น แต่ถ้าโอนโปรแกรมเข้าเครื่องคอมพิวเตอร์ผมพยายามหาในเน็ตตามฟอรั่ม มีคนพยายามจะแกะโปรโตคอลแต่ไม่น่าจะสำเร็จ ที่นี้จะเอาโปรแกรมมาแสดงบนคอมพิวเตอร์ได้ยังไง ก็ต้องนั่งแกะโปรแกรมทีละเม็ด เขียนด้วยเวิร์ดเพราะต้องการสัญลักษณ์ให้ตรงกับที่แสดงในเครื่องคิดเลขมากที่สุด สัญลักษณ์ส่วนใหญ่หาได้ในฟอนต์ Symbol, Wingdings
  • เมื่อแกะเสร็จแล้วก็ดูทวนอีกทีว่าพิมพ์ได้ตรงกับในเครื่องคิดเลขไหม ไม่มีอะไรตกหล่นก็ export มาเป็นภาพเอามาแปะ ผมพยายามแยกสีให้ดูง่ายตรงไหนเป็นฟังก์ชันใช้สีแดงเข้ม ตัวแปรสีน้ำเงิน เงื่อนไขโปรแกรมใช้สีเขียว ลองดูครับ
Geo2UTM
Geo2UTM
Geo2UTM(continued)
Geo2UTM(continued)
  • สำหรับเครื่องคิดเลขแล้ว ก็ไม่ถือว่าโปรแกรมใหญ่มากนัก แต่สังเกตดูตัวแปร ตั้งแต่ตัว A ถึงตัว Z ใช้แทบหมด ผมคงไม่อธิบายตัวโปรแกรมนะครบ จะยืดเยื้อ สำหรับคนที่เคยเขียนโปรแกรมเครื่องคิดเลขคาสิโอ้รุ่นเหล่านี้ มองแป๊ปเดียวก็โอเคแล้ว

วิธีใช้งาน

  • วิธีใช้งานก็ง่ายครับตามสไตล์เครื่องคิดเลข กดเรียกโปรแกรมก่อน Shift+Prog ที่เครื่องคิดเลขผมเลือก “GEO2UTM

20170108_122613

ตัวอย่างที่ 1

  • โปรแกรมจะถามค่าพิกัดแลตติจูดและค่าลองจิจูด ตัวอย่างการใช้นี้กำหนดให้ latitude = 26°12’3.6128″N longitude = 50°36’25.1928″E ที่เครื่องคิดเลขกดคีย์ “EXE” โปรแกรมจะถามค่าแลตติจูด ป้อนไปให้ทศนิยมครบ ตามรูปแรก  (เครื่องคิดเลขเวลาแสดงค่าพิกัดที่เราป้อนไปแล้ว จะแสดงแค่ทศนิยมสองตำแหน่ง) และป้อนค่าลองจิจูดให้ตามรูปถัดไป

20170108_123456

20170108_123609

ผลลัพธ์การคำนวณ

  • โปรแกรมจะคำนวณค่าพิกัดฉาก UTM ให้พร้อมบอกหมายเลขโซนของยูทีเอ็มมาด้วยและบอกว่าเป็นโซนด้านเหนือหรือด้านใต้ของเส้นศูนย์สูตร ในที่นี้จุดค่าพิกัดนี้แถวประเทศบาเรนห์ โซน 39N (เหนือ)

20170108_130154 20170108_130205

  • เทียบกับค่าที่คำนวณด้วยโปรแกรม Surveyor Pocket Tools ตรงกันครับ หมายเหตุนิดหนึ่งว่าการคำนวณการแปลงพิกัดในโปรแกรม  Surveyor Pocket Tools ใช้ไลบรารีของ Proj4 ผ่านทาง pyproj
Surveyor Pocket Tools
Surveyor Pocket Tools

ตัวอย่างที่ 2

  • มาลองดูกัน ถ้าผู้อ่านเกิดจับพลัดจับผลูไปทำงานต่างประเทศที่อยู่ใต้เส้นศูนย์สูตร ก็ยังสามารถใช้ได้ กำหนด แลตติจูด = 16d9’7.048″S ลองจิจูด = 33d33’49.779″E ค่าพิกัดอยู่ที่ประเทศโมซัมบิค ทวีปอาฟริกา
  • ป้อนค่าพิกัดเข้าดังนี้ ค่าแลตติจูดอยู่ใต้เส้นศูนย์สูตรให้ติดเครื่องหมายลบข้างหน้า ส่วนค่าพิกัดลองจิจูดก็ป้อนปกติ

20170108_132745
20170108_132903

  • ผลลัพธ์การแปลงพิกัดได้ดังนี้ จุดอยู่ที่โซน 36S ใต้เส้นศูนย์สูตร

20170108_132916
20170108_132923

  • เปรียบเทียบผลการคำนวณกับ Surveyor Pocket Tools ตรงกัน

surveyor-pocket-tools_2017-01-08_13-31-09

  • ก็พอหอมปากหอมคอครับ ตอนหน้ามาดูโปรแกรมแปลงค่าพิกัดฉากยูทีเอ็ม (UTM) ไปยังค่าพิกัดภูมิศาสตร์บ้าง ติดตามกันตอนต่อไปครับ

สนุกกับโปรแกรมเครื่องคิดเลขสำหรับงานสำรวจ ตอนที่ 2 โปรแกรมแปลงพิกัด “UTM2Geo” บนเครื่องคิดเลข Casio FX 5800P

โปรแกรม “UTM2Geo” สำหรับเครื่องคิดเลข Casio FX 5800P

  • สวัสดีครับผู้อ่านทุกท่าน พบกันตอนนี้เป็นตอนที่ 2 แล้วครับ ตอนแรกนำเสนอโปรแกรม “Geo2UTM” แปลงพิกัดจากคาพิกัดภูมิศาสตร์ (แลตติจูด/ลองจิจูด) ไปเป็นค่าพิกัดบนระบบพิกัดฉาก UTM มาตอนนี้กลับกันครับ เราจะมาเขียนโปรแกรมที่แปลงพิกัดจากระบบพิกัดฉากยูทีเอ็มไปเป็นค่าพิกัดภูมิศาสตร์
  • มาเขียนบทความที่นี่ให้กับ kns-engineering เป็นการเฉพาะกิจ สำหรับเพื่อนพี่น้องชาวสำรวจแวะไปเยี่ยมเยียนผมได้ที่บล็อกประจำ priabroy.name

เกร็ดความรู้เล็กน้อยสำหรับช่างสำรวจ

  • ก่อนจะเข้าไปว่าเรื่องโปรแกรมมิ่งบนเครื่องคิดเลข ขอซักซ้อมความรู้เซอร์เวย์สักเล็กน้อย มีสองอย่างคือมุม อะซิมัท (Azimuth) และฟังก์ชัน Atan2 สองอย่างนี้เกี่ยวพันกับเรื่องโปรแกรมในด้านการสำรวจเสียส่วนใหญ่ (แต่งานคำนวณแปลงพิกัดนี้ไม่ได้ใช้)
  • มุมอะซิมัทเป็นมุมที่แสดงทิศทางในด้านงานสำรวจของเรา เป็นมุมที่กวาดจากทิศเหนือตามเข็มนาฬิกา สมมติว่ามีจุด A มีค่าพิกัด (x,y,) = (500,500) และจุด B มีค่าพิกัด (x,y) = (586.603,550) คำถามพื้นฐานก็คือ
    • อะซิมัทจากจุด A ไปจุด B เท่าไหร่
    • อะซิมัทจากจุด B ไปจุด A เท่าไหร
    • ระยะทางหรือระยะราบเท่าไหร่ระหว่างจุดทั้งสอง

  • ระยะราบหาได้จากสูตรตรีโกณมิติสามัญพื้นฐาน ระยะราบ = √((500-586.603)² + (500-550)²) = 100 เมตร
  • แล้วอะซิมัทละคำนวณอย่างไร การคำนวณจะมาเกี่ยวพันกับเครื่องคิดเลขของเราอย่างแนบแน่น ฟังก์ชัน Atan2 เรียกกันบนเครื่องคอมพิวเตอร์ (ใครเขียนโปรแกรมบนเครื่องคอมพิวเตอร์จะรู้จักฟังก์ชั่นนี้ดี มีทุกภาษา) เที่ยบเท่าบนเครื่องคิดเลขก็คือ Pol() คือฟังก์ชั่นการที่จะมาช่วยย่นการคำนวณนี้
  • Atan2 ดั้งเดิมจะคำนวณหามุมที่กวาดจากแกน X ทวนเข็มนาฬิกา ดังนั้นในคู่มือคาสิโอ จะเขียนฟังก์ชั่น Pol() แบบนี้ครับ ในวงเล็บ x มาก่อน  y สังเกตว่ามุม θ กวาดจากแกน X

  • แต่สำหรับมุมที่ต้องการสำหรับงานสำรวจคือมุมอะซิมัท(Azimuth) คือมุมคือกวาดจากแกน Y ลงมาตามเข็มนาฬิกา จะทำอย่างไร
  • เทคนิคเวลาใช้งานสลับเอาไว้ค่า Y มาก่อนและ X ตามหลังครับ >> Pol(y,x) จะได้มุมอะซิมัท ผมเขียนรูปใหม่ดังนี้

  • ถ้ามีเครื่องคิดเลขก็ลองกดดูเลย Pol((550-500),(586.603-500))  ผลการคำนวณเครื่องคิดเลขเอาค่าระยะทางไปเก็บไว้ในตัว “I” และมุมอะซิมัทไว้ในตัว “J” ลอง RCL (recall) มาดูจะได้ I = 100.000 และ J = 60.000 มุม 60 ก็คืออะซิมัทจาก A ไป B นั่นเอง
  • จากโจทย์ข้างบนก็ตอบได้นะครับ อะซิมัทย้อน (Backward) ก็ให้เอาอะซิมัทไป (Forward) ± 180 ถ้ามากกว่า 180 ให้เอา 180 ไปลบ ถ้าน้อยกว่า 180 ให้เอา 180 ไปบวก ดังนั้นอะซิมัทจาก B มา A จะได้ = 60 + 180 = 240 
  • มุมอะซิมัทเป็นมุมมหัศจรรย์สำหรับช่างสำรวจ มหัศจรรย์ยังไงตอนหน้ามาว่ากันต่อ ตอนนี้ไปต่อเรื่องโปรแกรมกัน

โค๊ดโปรแกรม UTM2Geo

  • เหมือนเดิมนั่งแกะโปรแกรมทีละเม็ดแล้วเขียนลงในเวิร์ด แล้วตรวจทานสอง สามเที่ยว น้ำตาไหลพราก  ผมก็คิดเหมือนกันว่าผู้อ่านไปคีย์โปรแกรมตามที่ผมลิสต์ออกมาคงไม่ใช่เรื่องสนุก เอาอย่างนี้ ถ้าไม่อยากคีย์โปรแกรมเอง ติดต่อน้องๆที่ท้ายบล็อกผมจะลงชื่อไว้เพื่อขอดูดโปรแกรม มีข้อแม้เล็กๆน้อยๆว่าเครื่องต้องเป็นเครื่อง Casio Fx-5800P รุ่นเดียวกันและพกสายลิ๊งค์มาด้วยตัวเอง

ต้วอย่างที่ 1

  • สำหรับท่านที่ต้องการป้อนโปรแกรมลงในเครื่องคิดเลข วิธีป้อนโปรแกรมเลงเครื่องคิดเลขคงไม่ต้องสาธยายนะครับเพราะคงคุ้นเคยกันอยู่ มาดูการใช้งาน สมมติว่ามีหมุดที่มีค่าพิกัดฉาก UTM ดังนี้
    • Northing = 1,615,517.540 Easting = 395,698.272 Zone 47N อยู่บนพื้นหลักฐาน WGS84
  • เรียกโปรแกรม ด้วยการกด “Prog” แล้วเลือกโปรแกรม “UTM2Geo

  • ป้อนค่าพิกัด Norhting, Easting

  • ต่อไปโปรแกรมจะถาม Zone No ลำดับโซนหมายเลขที่เท่าไหร่ ถ้าโซนอยู่ด้านเหนือเส้นศูนย์สูตรให้ป้อนค่าบวกเช่น 47N ก็ป้อนตัวเลขธรรมดาไปครับ แต่ถ้าอยู่ใต้เส้นศูนย์สูตรป้อนเป็นตัวเลขลบ ในที่นี้ป้อนตัวเลข 47 เข้าไปดังรูป

ผลลัพธ์ของโปรแกรม

  • มาดูผลลัพธ์กันครับ

  • ผมเอาค่าพิกัดฉากนี้ไปป้อนในโปรแกรมแปลงพิกัด UTM-Geo Converter ที่อยู่ใน Surveyor Pocket Tools ก็ได้ตรงกันครับ แต่ทศนิยมในเครื่องคิดเลขที่ฟิลิปดาได้เต็มที่สองตำแหน่ง (น่าเสียดาย ทศนิยมของฟิลิปดาตำแหน่งที่สอง เทียบเป็นหน่วยเมตริกแล้วได้แค่ระดับหลักสิบเซนติเมตร) ถ้าต้องการหลักมิลลิเมตรก็ต้องบนคอมพิวเตอร์แล้วครับ แต่สำหรับเครื่องคิดเลขคิดมาได้ขนาดนี้ผมก็โอเคแล้วครับ

 ตัวอย่างที่ 2

  • สมมติว่าผู้อ่านมีโอกาสไปทำงานต่างประเทศไกลๆ มาลองค่าพิกัดที่อยู่โซนตะวันตกและอยู่ใต้เส้นศุนย์สูตรกันดูครับ สมมติว่าค่าพิกัดนี้อยู่ในบราซิลครับ ชีวิตจริงไม่เคยไปถึงทวีปอเมริกาครับ ไกลสุดแค่ยุโรปกับทวีปอาฟริกา
    •  Northing = 7,721,526.876 Easting = 505,464.207 Zone 24S บนพื้นหลักฐาน WGS84
  • พร้อมแล้วเรียกโปรแกรมป้อนตัวเลขกันเลย

 

  •  ต่อไปหมายเลขโซนให้ป้อนตัวเลขเป็นลบ -24 เพราะอยู่ใต้เส้นศูนย์สูตร

  • มาดูผลลัพธ์กัน จะได้ค่า Latitude = -20°36’19.28″ ค่าเป็นลบแสดงว่าอยู่ใต้เส้นศูนย์สูตร ได้ค่า Longitude = -38°56’51.220″ ได้ค่าเป็นลบแสดงว่าอยู่ทางด้านตะวันตกของตำบลกรีนนิช ของอังกฤษ (ตำบลกรีนนิชค่าลองจิจูด = 0)

  • เปรียบเทียบกับ Surveyor Pocket Tools เท่ากัน

  • ทรรศนะส่วนตัวผมมีโปรแกรมแปลงพิกัดติดเครื่องคิดเลขรุ่นนี้ FX – 5800P ทำให้เครื่องคิดเลขดูเทพขึ้นมาทันตาเห็น 🙂 สองโปรแกรมนี้กินเม็มไปจิ๊บๆครับ
  • ตอนหน้ามาว่าเรื่องยากขึ้นไปอีกนิด การคำนวณระยะทางที่สั้นที่สุดบนทรงรี (Geodesic distance) สูตรลากกันยาวเฟื้อย ตัวแปรบนเครื่องคิดเลขใช้หมดเกลี้ยงต้องไปดึงตัวอนุกรมมาช่วยด้วย ถ้าพิมพ์โปรแกรมตามผมต้องร้องว่า เจ็บกว่านี้มีอีกไหม
  • ตอนหน้ามาว่ากัน แต่ผมบอกก่อนว่าในฐานะช่างสำรวจ เรื่อง  geodesic distance บางครั้งเราใช้มันอย่างไม่รู้ตัว และไม่ใช่เรื่องไกลตัว ขอฝากน้องๆไว้ครับ มีความรู้ก็ใส่ตัวก็ใช่ว่าจะต้องไปเหนื่อยแบกหามตามโบราณที่ว่าไว้

ติดต่อขอลิ๊งค์โปรแกรม

  • ติดต่อคุณนันทวุฒิ อึ้งตระกูลได้ครับสำหรับลิ๊งค์ดูดโปรแกรม ผมจะฝากโปรแกรมไว้ที่นั่น แต่ถ้าใจไม่ร้อน รอบทความอีกสักสอง สามตอนจบก่อนก็ดี (เพราะบางโปรแกรมที่จะนำเสนอยังไม่ได้เขียนลงเครื่องคิดเลข)
  • สังคมจะน่าอยู่ ถ้าแบ่งปันกันและเอื้อเฟื้อเผื่อแผ่สำหรับคนที่มีโอกาสน้อยกว่า พบกันตอนใหม่ครับ

การเขียนโปรแกรมคำนวณการแปลงค่าพิกัดระหว่าง UTM Grid และ Geographic (Lat/Long) ด้วย Lazarus และ Delphi (ตอนที่ 3)

  • จากตอนที่ 2 จะเห็นโค๊ดที่ผม post 2 unit คือ GeoEllipsoids.pas และ GeoCompute.pas ถ้าสนใจก็ copy ไปวางที่ Text Editor แล้ว Save ชื่อไฟล์ตามที่ผม เราจะเริ่มต้นสร้าง New Project บน Lazarus จากนั้น Add file “GeoEllipsoids.pas” และ “GeoCompute.pas” เข้ามาในโปรเจคใหม่ที่สร้างขึ้น
  • ผมทิ้งท้ายเรื่องสร้างโปรเจคใหม่ในตอนที่ 1 ก่อนจะกล่าวถึง 2 unit ในตอนที่ 2 กลับมาที่ Lazarus อีกครั้ง
Lazarus 0.9.29 (beta) บน windows
Lazarus 0.9.29 (beta) บน windows
  • จะเห็นฟอร์ม unit1 นั้นเปล่าๆ ยังไม่ทีอะไร และ Project Inspector ด้านขวามือก็ยังว่างขึ้นเพียงแต่ project1.lpr และ unit1.pas เ่ท่านั้น ตรง component palette ยังมี component อีกมากมายของ Lazarus ที่ยังไม่ได้ติดตั้งเข้ามา ซึ่งจะอยู่ในไดเรคทอรี components ของ Lazarus สนใจตัวไหนก็ ใช้เมนูหลัก Packages > Open package file (.lpk) … ทำการ compile และ install ในตอนนี้ยังไม่ได้ใช้ตัวไหนเป็นพิเศษ จะใช้ standardd เท่านั้น
  • เมื่อสร้างโปรเจคใหม่แล้ว จุดที่สำคัญที่ต้องตั้งเลยก็คือ OS และ CPU Platform ที่เมนูหลักเลือก Project > Compiler Options… จะเห็น dialog ให้ตั้งค่าดังรูปด้านล่าง
ตั้งค่า OS และ Platform
ตั้งค่า OS และ Platform
  • ที่ Project Inspector ให้คลิกที่เครื่องหมาย + เพื่อ Add file…
Add file to project
Add file to project
  • Add file ที่ผม post ไว้ตอนที่ 2 เข้ามา 2 ไฟล์ดังรูปด้านล่าง
Project Inspector หลังจากเพิ่มไฟล์
Project Inspector หลังจากเพิ่มไฟล์
  • ที่ unit1 ทำการเพิ่ม object จาก standard palette เข้าไปตั้ง caption และชื่อดังรูปด้านล่าง
Create object on form.
Create object on form.
  • โค๊ดข้างล่างเป็นของ unit1 เขียน event เชื่อมต่อกับ ojbect บนฟอร์มแล้ว

[sourcecode language=”delphi”]
unit Unit1;

{$mode objfpc}{$H+}

interface

uses
Classes, SysUtils, FileUtil, LResources, Forms, Controls, Graphics, Dialogs,
StdCtrls, ExtCtrls, GeoEllipsoids, GeoCompute;

type

{ TForm1 }

TForm1 = class(TForm)
cmdEx1: TButton;
cmdEx2: TButton;
cmdCompute: TButton;
cmdClose: TButton;
cmbEllipsoids: TComboBox;
cmbHem: TComboBox;
Label7: TLabel;
txtZoneNo: TEdit;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
txtE: TEdit;
txtLat: TEdit;
Label1: TLabel;
Label2: TLabel;
radComputeType: TRadioGroup;
txtLong: TEdit;
txtN: TEdit;
procedure cmdCloseClick(Sender: TObject);
procedure cmdComputeClick(Sender: TObject);
procedure cmdEx1Click(Sender: TObject);
procedure cmdEx2Click(Sender: TObject);
procedure FormCreate(Sender: TObject);
procedure FormDestroy(Sender: TObject);
procedure radComputeTypeClick(Sender: TObject);
private
{ private declarations }
FEllipses : TEllipsoidList;
public
{ public declarations }
end;

var
Form1: TForm1;
function Coordinate(const X, Y : double) : TCoordinate;
implementation

function Coordinate(const X, Y : double) : TCoordinate;
begin
result.X := X;
result.Y := Y;
end;

{ TForm1 }

procedure TForm1.FormCreate(Sender: TObject);
var
i : integer;
begin
FEllipses := TEllipsoidList.Create;
for i := 0 to FEllipses.Count – 1 do
begin
cmbEllipsoids.Items.Add(TEllipsoid(FEllipses.Objects[i]).EllipsoidName);
end;
cmbEllipsoids.ItemIndex := FEllipses.Count – 1; //Default with WGS84
//start with UTM to Geographic computation.
txtN.Enabled := true;
txtE.Enabled := true;
txtZoneNo.Enabled := true;
cmbHem.Enabled := true;

txtLat.Readonly := true;
txtLong.Readonly := true;
end;

procedure TForm1.cmdCloseClick(Sender: TObject);
begin
Close;
end;

procedure TForm1.cmdComputeClick(Sender: TObject);
var
utm2geo : TUTM2Geo;
geo2utm : TGeo2UTM;
ell : TEllipsoid;
x, y : double;
utmproj : TUTMProjection;
geocoor, utmcoor : TCoordinate;
begin
ell := FEllipses.FindEx(cmbEllipsoids.Text);

if (radComputeType.ItemIndex = 0) then //UTM to Geographic
begin
x := strtofloat(txtE.Text);
y := strtofloat(txtN.Text);
utmcoor := Coordinate(x, y);
utmproj.ZoneNo := strtoint(txtZoneNo.Text);
if (cmbHem.ItemIndex = 0)then
utmproj.LatHem := hemNorth
else
utmproj.LatHem := hemSouth;

try
utm2geo := TUTM2Geo.Create;
utm2geo.Ellipsoid := ell;
utm2geo.UTMProjection := utmproj;
utm2Geo.UTMCoordinate := utmcoor;
utm2Geo.Compute;
txtLat.Text := format(‘%11.8f’, [utm2geo.GeoCoordinate.Y]);
txtLong.Text := format(‘%11.8f’, [utm2geo.GeoCoordinate.X]);
finally
utm2geo.Free;
end;
end
else  //Geographic to UTM
begin
x := strtofloat(txtLong.Text);
y := strtofloat(txtLat.Text);
geocoor := Coordinate(x, y);

try
geo2utm := TGeo2UTM.Create;
geo2utm.Ellipsoid := ell;
geo2utm.GeoCoordinate := geocoor;
geo2utm.Compute;
utmproj := geo2utm.UTMProjection;
txtN.Text := format(‘%11.3f’, [geo2utm.UTMCoordinate.Y]);
txtE.Text := format(‘%11.3f’, [geo2utm.UTMCoordinate.X]);
txtZoneNo.Text := format(‘%3d’, [utmproj.ZoneNo]);
if (utmproj.LatHem = hemNorth) then
cmbHem.ItemIndex := 0
else
cmbHem.ItemIndex := 1;
finally
geo2utm.Free;
end;
end
end;

procedure TForm1.cmdEx1Click(Sender: TObject);
begin
radComputeType.ItemIndex := 0; //UTM to Geographic
txtN.Text := ‘3885802.712’;
txtE.Text := ‘468327.864’;
txtZoneNo.Text := ’11’; //UTM Zone = 11
cmbHem.ItemIndex := 0; //North

txtLat.Text := ”;
txtLong.Text := ”;
end;

procedure TForm1.cmdEx2Click(Sender: TObject);
begin
radComputeType.ItemIndex := 1; //Geographic to UTM
txtLat.Text := ‘-37.6528211388889’;
txtLong.Text := ‘143.926495527778’;

txtN.Text := ”;
txtE.Text := ”;
txtZoneNo.Text := ”;
end;

procedure TForm1.FormDestroy(Sender: TObject);
begin
FEllipses.Free;
end;

procedure TForm1.radComputeTypeClick(Sender: TObject);
begin
if (radComputeType.ItemIndex = 0) then // UTM to Geographic
begin
cmdCompute.Caption := ‘&lt;==’;
txtN.Enabled := true;
txtE.Enabled := true;
txtZoneNo.Enabled := true;
cmbHem.Enabled := true;
txtN.SetFocus;
txtLat.Readonly := true;
txtLong.Readonly := true;
end
else
begin
cmdCompute.Caption := ‘==&gt;’;   //Geographic to UTM
txtN.Readonly := true;
txtE.Readonly := true;
txtZoneNo.Readonly := true;
cmbHem.Readonly := true;
txtLat.Enabled := true;
txtLong.Enabled := true;
txtLat.SetFocus;
end;
end;

initialization
{$I unit1.lrs}

end.
[/sourcecode]

  • ดูที่ procedure TForm1.FormCreate(Sender: TObject); เป็น event ของ form1 ผมดึง Ellipsoids (FEllipsoids) ที่สร้างไว้เป็น object มา iterate ด้วย for loop ลงไปที่ combo box ชื่อ cmbEllipsoids เวลาผู้ใช้คลิกเลือก Ellipsoid จาก combo box จะได้เพียงแต่ text ออกมาเราจะใช้ text ค้นกลับไปที่ FEllipses( TEllipsoidList) เพื่อดึง object ที่แสดงรูปทรงของ Ellipsoid
  • เมื่อผู้ใช้เลือกคลิกที่ radio group “Computation type” จะเลือกว่าจะคำนวณจาก UTM ไป Geographic หรือในทางกลับกัน อีเวนท์ของ radio group นี้อยู่ที่ procedure TForm1.radComputeTypeClick(Sender: TObject); จะได้ใช้ if clause เลือกการคำนวณถูกว่าจะใช้ class ไหน ระหว่าง TUTM2Geo หรือ TGeo2UTM
  • เพื่อให้ดูง่ายยิ่งขึ้นผมเพิ่มตัวอย่างเพื่อจะใช้คำนวณดู โดยกำหนดอยู่ที่ button ชื่อ cmdEx1 และ cmdEx2 พร้อมทั้งเขียนอีเวนท์ให้ 2 procedure คือ procedure TForm1.cmdEx1Click(Sender: TObject); และ procedure TForm1.cmdEx2Click(Sender: TObject); ตามลำดับ

[sourcecode language=”delphi”]
procedure TForm1.cmdEx1Click(Sender: TObject);
begin
radComputeType.ItemIndex := 0; //UTM to Geographic
txtN.Text := ‘3885802.712’;
txtE.Text := ‘468327.864’;
txtZoneNo.Text := ’11’; //UTM Zone = 11
cmbHem.ItemIndex := 0; //North

txtLat.Text := ”;
txtLong.Text := ”;
end;

procedure TForm1.cmdEx2Click(Sender: TObject);
begin
radComputeType.ItemIndex := 1; //Geographic to UTM
txtLat.Text := ‘-37.6528211388889’;
txtLong.Text := ‘143.926495527778’;

txtN.Text := ”;
txtE.Text := ”;
txtZoneNo.Text := ”;
end;
[/sourcecode]

  • สุดท้าย button ที่เป็นรูปลูกศรชี้ไปซ้ายเมื่อเลือกคำนวณจาก UTM ไป Geographic และชี้ไปด้านขวาเมื่อเลือก Geographic ไป UTM เมื่อป้อนค่าพิกัดที่เหมาะสมแล้วทำการคลิกที่ button ก็จะไปปลุกอีเวนท์ ดังโค๊ดด้านล่าง

[sourcecode language=”delphi”]
procedure TForm1.cmdComputeClick(Sender: TObject);
var
utm2geo : TUTM2Geo;
geo2utm : TGeo2UTM;
ell : TEllipsoid;
x, y : double;
utmproj : TUTMProjection;
geocoor, utmcoor : TCoordinate;
begin
ell := FEllipses.FindEx(cmbEllipsoids.Text);

if (radComputeType.ItemIndex = 0) then //UTM to Geographic
begin
x := strtofloat(txtE.Text);
y := strtofloat(txtN.Text);
utmcoor := Coordinate(x, y);
utmproj.ZoneNo := strtoint(txtZoneNo.Text);
if (cmbHem.ItemIndex = 0)then
utmproj.LatHem := hemNorth
else
utmproj.LatHem := hemSouth;

try
utm2geo := TUTM2Geo.Create;
utm2geo.Ellipsoid := ell;
utm2geo.UTMProjection := utmproj;
utm2Geo.UTMCoordinate := utmcoor;
utm2Geo.Compute;
txtLat.Text := format(‘%11.8f’, [utm2geo.GeoCoordinate.Y]);
txtLong.Text := format(‘%11.8f’, [utm2geo.GeoCoordinate.X]);
finally
utm2geo.Free;
end;
end
else  //Geographic to UTM
begin
x := strtofloat(txtLong.Text);
y := strtofloat(txtLat.Text);
geocoor := Coordinate(x, y);

try
geo2utm := TGeo2UTM.Create;
geo2utm.Ellipsoid := ell;
geo2utm.GeoCoordinate := geocoor;
geo2utm.Compute;
utmproj := geo2utm.UTMProjection;
txtN.Text := format(‘%11.3f’, [geo2utm.UTMCoordinate.Y]);
txtE.Text := format(‘%11.3f’, [geo2utm.UTMCoordinate.X]);
txtZoneNo.Text := format(‘%3d’, [utmproj.ZoneNo]);
if (utmproj.LatHem = hemNorth) then
cmbHem.ItemIndex := 0
else
cmbHem.ItemIndex := 1;
finally
geo2utm.Free;
end;
end
end;
[/sourcecode]

  • จากโค๊ดด้านบนซึ่งจะทำหน้าที่คำนวณเมื่อผู้ใช้คลิกที่ปุ่มคำนวณรูปลูกศร ผม declare local variable ไว้ 2 ตัวแปร คือ utm2geo และ geo2utm ตามแต่จะเลือกคำนวณ
  • มีตัวแปรที่สำคัญคือ ell : TEllipsoid เมื่อเริ่มต้นจะดึง object ของ TEllipsoid โดยใช้ชื่อที่ผู้ใช้คลิกที่ combo box “Ellipsoid” เป็นตัว index ในการค้นหา ดังโค๊ด ell := FEllipses.FindEx(cmbEllipsoids.Text);
  • ตัวแปร utmcoor และ geocoor ดึงค่าพิกัดที่ผู้ใช้ป้อนลง text edit
  • คงไม่ยากที่จะเข้าใจเพราะ class ที่ทำหน้าที่คำนวณ แยกออกไปต่างหากแล้ว ดูผลลัพธ์ที่ได้จากการคลิกปุ่ม Example 1 และ Example 2 ตามลำดับ

 

UTM to Geographic.
UTM to Geographic.
Geographic to UTM.
Geographic to UTM.
  • ผมตั้งใจจะรวมไฟล์โปรเจคนี้เป็น zip ไฟล์ แต่บริการแชร์ไฟล์ของ WordPress คือ Box.net ยังปิดบริการชั่วคราวครับ ตอนหน้าถ้าเป็น programming ด้วย Lazarus ผมอยากจะเขียนเรื่อง การคำนวณหาระยะทางและอะซิมัทบนรูปทรงรี (Geodesic distance & Azimuth on Ellipsoid) เมื่อ input เป็นค่าพิกัด Lat/Long ทั้ง 2 จุด (เป็นค่าพิกัด UTM ก็ได้แต่ต้องคำนวณเป็น Geographic ด้วย class ที่ผมเสนอในตอนนี้นั่นเอง)  ผู้อ่านเคยสงสัยบ้างไหมครับว่าเครื่อง GPS เช่น Garmin เมื่อเราป้อน Waypoint ก็คือจุดที่มีค่าพิกัด ผู้ใช้สามารถเลือกป้อนเป็นค่า Latitude,Longitude ก็ได้ แต่ถ้าป้อนเป็นยูทีเอ็มจะสังเกตเห็นว่าเครื่อง GPS จะบังคับให้เราป้อนโซนของยูทีเอ็มเช่น 47P, 47Q อะไรประมาณนี้ (ตัวเลข P,Q เรียกว่า UTM Sub Zone เริ่มจากตัว C ที่ใกล้ขั้วโลกใต้ขึ้นไปหาตัว X ที่ขั้วโลกเหนือ) เครื่อง GPS จะฝังฟังก์ชั่นที่กล่าวนี้เพื่อหาระยะทางและอะซิมัทบน Ellipsoid และก็เช่นเคยสูตรพวกนี้นักคณิตศาสตร์ได้คิดกันมาก่อนแล้วก็มีท่านอื่นปรับปรุงมาเรื่อยๆให้ได้ความถูกต้องเพิ่มขึ้น
UTM sub zone.
UTM sub zone.
  • คำถามคือทำไมไม่หาระยะทางและอะซิมัทจากค่าพิกัดยูทีเอ็มโดยตรง คำตอบคือหาได้ แต่ต้องอยู่บนโซนเดียวกัน แต่ค่าที่คำนวณจะผิดได้ ถ้าไม่มี UTM Zone No กำกับค่าพิกัดมาให้ เครื่อง GPS จึงบังคับให้เราป้อนค่าโชนดังที่ได้กล่าวไปแล้ว

Download sourcecode

  • sourcecode ของโปรแกรมตอนนี้สามารถ download ได้ที่นี่ UTM2LatLongProj.zip

การเขียนโปรแกรมคำนวณการแปลงค่าพิกัดระหว่าง UTM Grid และ Geographic (Lat/Long) ด้วย Lazarus และ Delphi (ตอนที่ 2)

สูตรที่ใช้ในการคำนวณ

  • ผมอ้างอิงมาจาก http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm เขียนโดย Steve Dutch กล่าวถึงสูตรที่นำมาใช้คำนวณแปลงค่าพิกัด โดยอ้างถึง U.S. Army Corps of Engineer และ USGS (U.S. Geological Survey Professional) ถ้าเป็นการแปลงจาก Lat/Long ไป UTM เอกสารของ U.S. Army Corps of Engineer ทำได้ดีกว่าชัดเจนกว่า แต่ถ้าเป็นการแปลงพิกัดจาก UTM ไป Lat/Long USGS ทำได้ดีกว่า
  • เราจะใช้สูตรที่ Steve Dutch ได้อ้างถึงและแสดงไว้ในเอกสาร นำมาโค๊ดเป็น Lazarus (Delphi) และที่น่าสนใจคือมีไฟล์ของ Excel สำหรับการคำนวณไว้ให้ดาวน์โหลดได้ด้วย ทำไว้ได้ดีมากครับ อยู่ตรงนี้ http://www.uwgb.edu/dutchs/UsefulData/UTMConversions1.xls

UTM Zone

UTM Zone No
UTM Zone No

GeoEllipsoids Unit

  • เราจะเริ่มด้วย Unit “GeoEllipsoids” (เก็บอยู่ในไฟล์ geoellipsoids.pas ที่ Lazarus พยายามให้เราตั้งชื่อเป็น lower case) เป็นยูนิตที่ declare class รูปแบบของ Ellipsoid ชื่อ TEllipsoid และมีสร้าง collection สำหรับเก็บ Ellipsoid ไว้ด้วยในยูนิตเดียวกัน
      [sourcecode language=”delphi”]

 

      unit GeoEllipsoids;

 

      {$mode objfpc}{$H+}

 

      interface

 

      uses

 

      Classes, Contnrs, SysUtils;

 

      type

 

      TEllipsoid = class

 

      private

 

      FEllipName : string;

 

      FMajorAxis : double;

 

      FInvFlat : double;

 

      public

 

      property EllipsoidName : string read FEllipName write FEllipName;

 

      property MajorAxis : double read FMajorAxis write FMajorAxis;

 

      property InvFlattening : double read FInvFlat write FInvFlat;

 

      constructor Create ;

 

      destructor Destroy ; override ;

 

      end;

TEllipsoidList = class (TStringList)
private
protected
function GetEllipsoid(idx : Integer ) : TEllipsoid ;
public
constructor Create ;
destructor Destroy ; override ;
procedure AddEx (const Name : string; const MajorAxisA : double; const InvFlat : double);
function FindEx(const Name : string): TEllipsoid ;
end;

implementation
//TEllipsoid
constructor TEllipsoid.Create;
begin
inherited create;
end;

destructor TEllipsoid.Destroy;
begin
inherited Destroy;
end;

//TEllipsoidList – Collection of TEllipsoid.
constructor TEllipsoidList.Create ;
begin
inherited Create;
Sorted := True ;

AddEx(‘WGS84’, 6378137, 298.257223563);
AddEx(‘GRS80’, 6378137, 298.257222101);
AddEx(‘WGS72’, 6378135, 298.260);
AddEx(‘Australian 1965’, 6378160, 298.250);
AddEx(‘Krasovsky 1940’, 6378245, 298.3);
AddEx(‘International (1924)’, 6378388, 297);
AddEx(‘Clake 1880’, 6378249.1, 293.460);
AddEx(‘Clarke 1866’, 6378206.4, 294.980);
AddEx(‘Airy 1830’, 6377563.4, 299.320);
AddEx(‘Bessel 1841’, 6377397.2, 299.150);
AddEx(‘Everest 1830’, 6377276.345, 300.8017);
AddEx(‘Hayford 1909’, 6378388.0, 296.999362081575);
AddEx(‘North American 1927’, 6378206.4, 294.978698213898);
AddEx(‘NAD 83’, 6378137, 298.257223563);
end ;

destructor TEllipsoidList.Destroy ;
var
i : Integer ;
begin
for i:= Count -1 downto 0 do
begin
if Assigned( Objects[i] ) then
TEllipsoidList( Objects[i] ).Free ;
end ;
inherited Destroy;
end ;

procedure TEllipsoidList.AddEx(const Name : string ; const MajorAxisA : double; const InvFlat : double);
var
e : TEllipsoid ;
pos : Integer ;
begin
if Find(Name, pos ) then
begin
if Assigned(Objects[pos]) then
TEllipsoid(Objects[pos]).Free ;
end ;
pos := Add(Name);
e := TEllipsoid.Create ;
with e do
begin
EllipsoidName := Name;
MajorAxis := MajorAxisA;
InvFlattening := InvFlat;
end ;
Objects[pos] := e;
end;

function TEllipsoidList.FindEx(const Name : String): TEllipsoid ;
var
pos : Integer ;
begin
if Find( Name, pos ) then
Result := TEllipsoid(Objects[pos])
else
Result := nil ;
end ;

function TEllipsoidList.GetEllipsoid(idx : Integer) : TEllipsoid ;
begin
Result := TEllipsoid( Objects[ idx ] );
end ;

end.
[/sourcecode]

  • จากโค๊ตด้านบนจะเห็นผม Declare class ชื่อ TEllipsoid มี property เก็บค่าของรูปทรงรี เนื่องจากค่าของรูปทรงรีที่คนใช้กันเป็นค่ามาตรฐานผมจึงนำ object ของ class “TEllipsoid” มาสร้างเก็บไว้ใน class ที่สืบทอดจาก TStringList ให้ชื่อว่า TEllipsoidList ทำไมต้อง derive จาก TStringList เนื่องจาก TStringList ออกแบบมาให้เก็บชุดของตัวอักษรแต่สามารถเก็บ Object ได้  มี Index ที่สามารถค้นหาได้ด้วยเมธอด Find ความหมายง่ายๆคือเราเอา Object ของ class “TEllipsoid” มายัดใส่ใน list ของ TEllipsoidList(derive มาจาก TStringList) โดยอาศัยชื่อของทรงรี (EllipsoidName) มาเป็นตัว index เอาไว้ค้นหา Object
  • ผมเพิ่มเมธอด FindEx สำหรับค้นหาและดึง object ของทรงรีออกมา และยังเพิ่ม AddEx สำหรับเผื่อไว้ในกรณีต้องการเพิ่มทรงรีอื่นในภายหลังได้

GeoCompute Unit

  • ต่อไปเป็น Unit ที่ 2 ชื่อ “GeoCompute” (เก็บอยู่ในไฟล์ geocompute.pas)

[sourcecode language=”delphi”]
unit GeoCompute;

{$mode objfpc}{$H+}

interface

uses
Classes, SysUtils, GeoEllipsoids, Math;

type
TZoneHemisphere = (hemEast, hemWest, hemNorth, hemSouth);

TUTMProjection = record
K0 : double;
FE : double;
FN : double;
ZoneNo : integer;
CM : double;
LatHem : TZoneHemisphere;
Longhem : TZoneHemisphere;
end;

TCoordinate = record
X : double;
Y : double;
end;

TGeo2UTM = class
private
FEllipsoid : TEllipsoid;
FGeoCoor : TCoordinate;
FUTMCoor : TCoordinate;
FUTMProj : TUTMProjection;
procedure SetGeoCoordinate (const geocoor : TCoordinate);
procedure SetEllipsoid(const ellipse : TEllipsoid);
public
property Ellipsoid : TEllipsoid write SetEllipsoid;
property GeoCoordinate : TCoordinate write SetGeoCoordinate; //input
property UTMCoordinate : TCoordinate read FUTMCoor;  //output
property UTMProjection : TUTMProjection read FUTMProj; //output
procedure Compute;
constructor Create;
destructor Destroy; override;
end;

//On UTM Grid coordianates we need to know the zone no, latitude hemisphere and
//longitude hemisphere also. Then we can convert utm coordinates to geographic.
TUTM2Geo = class
private
FEllipsoid : TEllipsoid;
FGeoCoor : TCoordinate;
FUTMCoor : TCoordinate;
FUTMProj : TUTMProjection;
procedure SetUTMProjection (const utmproj : TUTMProjection);
procedure SetUTMCoordinate (const utmcoor : TCoordinate);
procedure SetEllipsoid(const Ellipsoid : TEllipsoid);
public
property UTMProjection : TUTMProjection read FUTMProj write SetUTMProjection;
property Ellipsoid : TEllipsoid write SetEllipsoid;
property UTMCoordinate : TCoordinate write SetUTMCoordinate;  //input
property GeoCoordinate : TCoordinate read FGeoCoor; //output
procedure Compute;
constructor Create;
destructor Destroy; override;
end;

implementation

//TGeo2UTM
procedure TGeo2UTM.SetEllipsoid(const ellipse : TEllipsoid);
begin
FEllipsoid.EllipsoidName := ellipse.EllipsoidName;
FEllipsoid.MajorAxis := ellipse.MajorAxis;
FEllipsoid.InvFlattening := ellipse.InvFlattening;
end;

procedure TGeo2UTM.SetGeoCoordinate(const geocoor : TCoordinate);
var
dLat, dLong : double;
begin
FGeoCoor := geocoor;
dLat := FGeoCoor.Y;
dLong := FGeoCoor.X;

//set hemisphere
if (dLat &gt;=0) then
FUTMProj.LatHem := hemNorth
else
FUTMProj.LatHem := hemSouth;

if (dLong &gt;=0) then
FUTMProj.LongHem := hemEast
else
FUTMProj.LongHem := hemWest;

//compute zone
if (FUTMProj.LongHem = hemWest) then
FUTMProj.ZoneNo := trunc((180 + dLong) / 6) + 1
else if (FUTMProj.LongHem = hemEast) then
FUTMProj.ZoneNo := trunc(dLong/6) + 31;

//compute zone
if (FUTMProj.LatHem = hemNorth) then
FUTMProj.FN := 0
else if (FUTMProj.LatHem = hemSouth) then
FUTMProj.FN := 10000000;

FUTMProj.CM := 6 * FUTMProj.ZoneNo – 183;
end;

procedure TGeo2UTM.Compute;
var
K0, FE, FN : extended;
lat, long : extended; //in radian.
a, b, f, e, et2, rm, n, rho, nu, p, cm : extended;
M, A0, B0, C0, D0, E0 : extended;
Ki, Kii, Kiii, Kiv, Kv, A6 : extended;

begin
//Assign to new variable for easy looking at mathematic formula.
K0 := FUTMProj.K0;
FE := FUTMProj.FE;
FN := FUTMProj.FN;
a := FEllipsoid.MajorAxis;
f := 1 / FEllipsoid.InvFlattening;
long := PI/180 * FGeoCoor.X; //convert latitude to radian.
lat := PI/180 * FGeoCoor.Y;
cm := PI/180 * FUTMProj.CM;
p := long – cm; //already in radian.

b := a * (1 – f); //Minor axis;

// =.08 approximately. This is the eccentricity of the earth’s elliptical crosssection.
e := sqrt(1 – b*b / (a*a));

//= .007 approximately. The quantity e’ only occurs in even powers so it need only be calculated as e’2
et2 := power((e * a/b), 2);

n := (a – b) / (a + b);

//This is the radius of curvature of the earth in the meridian plane.
rho := a * (1 – e*e) / power (1 – power(e*sin(lat),2), 3/2);

//This is the radius of curvature of the earth perpendicular to the meridian plane.
nu := a / sqrt(1 – sqr(e*sin(lat)));

//=================== compute Meridian Arc===================================
A0 := a * (1 – n + (5.0/4.0)*(n*n – n*n*n) + (81.0/64.0)*(power(n,4) – power(n,5)));
B0 := (3*a*n/2) * (1 – n + (7/8)*(n*n – n*n*n) + (55/64)* power(n,4) – power(n,5));
C0 := (15*a*n*n/16) * (1 – n + (3/4)*(n*n – n*n*n));
D0 := (35*a*n*n*n/48) * (1 – n + (11/16)*(n*n – n*n*n));
E0 := (315*a*power(n,4)/512) * (1-n);
M := A0*lat – B0*sin(2*lat) + C0*sin(4*lat) – D0*sin(6*lat) + E0*sin(8*lat);
//===========================================================================

//===================Converting Lat/Long to UTM Grid Coordinate==============
Ki := M * K0;
Kii := K0 * nu * sin(2*lat) / 4;
Kiii := (K0 * nu * sin(lat) * power(cos(lat),3) / 24) * ((5 – tan(lat)*tan(lat)
+ 9*et2 * cos(lat)*cos(lat) + 4 * et2*et2 * power(cos(lat),4)));
Kiv := K0 * nu * cos(lat);
Kv := (K0 * nu * power(cos(lat),3)/6) * (1 – tan(lat)*tan(lat) + et2*cos(lat)*cos(lat));
//Now we ‘v got Grid UTM Coordinates.
FUTMCoor.Y := FN + Ki + Kii*p*p + Kiii*power(p,4);
FUTMCoor.X := FE + Kiv*p + Kv*p*p*p;
end;

constructor TGeo2UTM.Create;
begin
inherited create;
FUTMProj.K0 := 0.9996;
FUTMProj.FE := 500000;
FEllipsoid := TEllipsoid.Create;
end;

destructor TGeo2UTM.Destroy;
begin
FEllipsoid.Free;
inherited Destroy;
end;

//TUTM2Geo
procedure TUTM2Geo.SetEllipsoid(const Ellipsoid : TEllipsoid);
begin
FEllipsoid.EllipsoidName := Ellipsoid.EllipsoidName;
FEllipsoid.MajorAxis := Ellipsoid.MajorAxis;
FEllipsoid.InvFlattening := Ellipsoid.InvFlattening;
end;

procedure TUTM2Geo.SetUTMProjection(const utmproj : TUTMProjection);
begin
FUTMProj.LongHem := utmproj.LongHem;
FUTMProj.LatHem := utmproj.LatHem;
if (utmproj.LatHem = hemNorth) then
FUTMProj.FN := 0
else if (utmproj.LatHem = hemSouth) then
FUTMProj.FN := 10000000;

FUTMProj.CM := 6 * utmproj.ZoneNo – 183;

if (FUTMProj.CM &lt; 0) then
FUTMProj.Longhem := hemWest
else
FUTMProj.Longhem := hemEast;

end;

procedure TUTM2Geo.SetUTMCoordinate(const utmcoor : TCoordinate);
var
dN, dE : double;
begin
FUTMCoor := utmcoor;
dN := FGeoCoor.Y;
dE := FGeoCoor.X;
end;

procedure TUTM2Geo.Compute;
var
K0, FE, FN, a, b, f, e, et2, cm : double;
mu, M, e1 : double;
J1, J2, J3, J4 : double;
C1, T1, R1, N1, D : double;
fp, Q1, Q2, Q3, Q4, Q5, Q6, Q7 : double;
lat, long : double;
x, y : double;
begin
//Assign to new variable for easy looking at mathematic formula.
K0 := FUTMProj.K0;
FE := FUTMProj.FE;
FN := FUTMProj.FN;
a := FEllipsoid.MajorAxis;
f := 1 / FEllipsoid.InvFlattening;
x := (FE – FUTMCoor.X);
if (FUTMProj.LatHem = hemNorth) then
y := FUTMCoor.Y
else if (FUTMProj.LatHem = hemSouth) then
y := (10000000 – FUTMCoor.Y);  //temporary

cm := PI/180 * FUTMProj.CM;
b := a * (1 – f); //Minor axis;

// =.08 approximately. This is the eccentricity of the earth’s elliptical crosssection.
e := sqrt(1 – b*b / (a*a));

//= .007 approximately. The quantity e’ only occurs in even powers so it need only be calculated as e’2
et2 := power((e * a/b), 2);

M := y / K0; // M is Meridian Arc;
//Compute footprint Latitude.
mu := M/(a*(1 – e*e/4 – 3/64*power(e,4) – 5/256*power(e,6)));
e1 := (1 – sqrt(1 – e*e))/(1 + sqrt(1 – e*e));
J1 := (3*e1/2 – 27*e1*e1*e1/32);
J2 := 21*e1*e1/16 – 55*power(e1,4)/32;
J3 := 151*e1*e1*e1/96;
J4 := 1097*power(e1,4)/512;
//Footprint Latitude fp = mu + J1sin(2mu) + J2sin(4mu) + J3sin(6mu) + J4sin(8mu)
fp := mu + J1*sin(2*mu) + J2*sin(4*mu) + J3*sin(6*mu) + J4 * sin(8*mu);

C1 := et2*cos(fp)*cos(fp);
T1 := tan(fp)*tan(fp);
R1 := a*(1-e*e)/power(1-e*e*sin(fp)*sin(fp), 1.5);
N1 := a / sqrt(1 -e*e*sin(fp)*sin(fp));
D := x / (N1 * K0);
//Compute Latitude
Q1 := N1 * tan(fp)/R1;
Q2 := D*D/2;
Q3 := (5 + 3*T1 + 10*C1 – 4*C1*C1-9*et2)*power(D,4)/24;
Q4 := (61 + 90*T1 + 298*C1 + 45*T1*T1 – 3*C1*C1 – 252*et2)* power(D,6)/720;

Q5 := D;
Q6 := (1 + 2*T1 + C1) *D*D*D/6;
Q7 := (5 – 2*C1 + 28*T1 – 3*C1*C1 + 8*et2 + 24*T1*T1)*power(D,5)/120;

//lat
lat := 180/PI *(fp – Q1*(Q2 + Q3 + Q4));

//long = long0 + (Q5 – Q6 + Q7)/cos(fp)
long := 180/PI *(CM – (Q5 – Q6 + Q7)/cos(fp));
FGeoCoor.Y := lat;
FGeoCoor.X := long;
end;

constructor TUTM2Geo.Create;
begin
inherited create;
FEllipsoid := TEllipsoid.Create;
FUTMProj.K0 := 0.9996;
FUTMProj.FE := 500000;
end;

destructor TUTM2Geo.Destroy;
begin
FEllipsoid.Free;
inherited Destroy;
end;

end.
[/sourcecode]

  • จากโค๊ดด้านบนผม declare type ที่เป็น record ไว้ 2 แบบ คือ
    1. TUTMProjection เพื่อช่วยกำหนดตำแหน่งของ UTM Zone เช่น Zone No., CM (Central Meridian), FE (False Easting เท่ากับ 500000 เสมอ), FN(False Northing = 0 ถ้าอยู่ด้านเหนือเส้นศูนย์สูตร และเท่ากับ 10000000 ถ้าอยู่ด้านใต้เส้นศูนย์สูตร), K0 (Scale Factor = 0.9996)
    2. TCoordinate เพื่อเก็บค่าพิกัดซึ่งจะมี 2 ค่าคือ Latitude, Longitude หรือ X,Y (UTM)
  • ถัดลงมาอีกผม declare class ไว้ 2 class คือ TUTM2Geo และ TGeo2UTM มาดูรายละเอียด
    1. TGeo2UTM ทำหน้าที่คำนวณค่าพิกัดจาก Geographic ไปเป็น UTM
      • property สำคัญคือ Ellipsoid เราต้อง set โดยการส่ง object ของ class TEllipsoid ที่ผู้ใช้เลือกมาให้
      • property อีกตัวคือ GeoCoordinate เป็นค่าพิกัด Lat, Long ที่เราส่งไปให้เพื่อจะแปลงค่าพิกัดเป็นพิกัดยูทีเอ็ม
      • Compute เป็น method สำคัญที่ทำหน้าที่คำนวณแปลงค่าพิกัด ในเมธอดนี้จะมีสูตรที่ผมอ้างอิงถึงดังที่กล่าวไปแล้ว
      • ส่วนผลลัพธ์ส่งออกมาทาง property UTMCoordinate จะเป็นค่าพิกัดยูทีเอ็มที่เราต้องการ
      • และ property UTMProjection ก็เป็นผลลัพธ์อีกตัวที่ได้จากการคำนวณ ข้อดีของค่าพิกัดในรูป Geographic (Latitude/Longitude) สามารถคำนวณหาได้ว่าอยู่ Zone ไหนของยูทีเอ็ม, หา Central Meridian ได้
    2. TUTM2Geo ทำหน้าที่คำนวณค่าพิกัดจาก UTM ไปเป็น Geographic
      • Ellipsoid property เช่นเดียวกันกับ class TGeo2UTM
      • UTMProjection property ค่าพารามิเตอร์ของยูทีเอ็มที่เราส่งไปให้ เพื่อทำการคำนวณแปลงพิกัดจาก UTM ไปยัง Geographic จำเป็นต้องรู้พารามิเตอร์ของยูทีเอ็ม เช่น Zone No., Hemisphere อยู่เหนือหรือใต้เส้นศูนย์สูตร (ค่าพิกัดยูทีเอ็มค่าพิกัดเดียวกัน อาจจะมีซ้ำกันหลายๆโซน ดังนั้นถ้าบอกพิกัดยูทีเอ็มต้องบอก Zone No. และ Hemisphere ด้วย)
      • UTMCoordinate property เป็นค่าพิกัดยูทีเอ็มที่ต้องการแปลงค่าพิกัดไปยัง Geographic
      • Compute method เช่นเดียวกันเมื่อ input คือ Ellipsoid, UTMprojection และ UTMCoordinate พร้อมก็เรียกใช้ method นี้ได้เลยเพื่อทำการคำนวณ
      • GeoCoordinate property เป็น output คือค่าพิกัด Latitude และ Longitude ที่เราต้องการนั่นเอง
  • ทิ้งท้ายไว้ตรงนี้ครับ ตอนหน้าเป็นตอนสุดท้าย 🙂

การเขียนโปรแกรมคำนวณการแปลงค่าพิกัดระหว่าง UTM Grid และ Geographic (Lat/Long) ด้วย Lazarus และ Delphi (ตอนที่ 1)

ตัวอย่างโปรแกรมแปลงค่าพิกัด GeoCalc และ GeoTrans

  • ผู้ใช้งานด้าน GIS และคนที่ทำงานด้านสำรวจ (Surveying) คงหลีกหนีการแปลงพิกัดระหว่าง UTM Grid และ Geographic (Latitude, Longitude) ได้ยาก ทั้งทางตรงและทางอ้อม ทางตรงได้แก่ ใช้โปรแกรมแปลงพิกัดเช่น GeoCalc (ฟรี แต่ไม่ Opensource) สามารถดาวน์โหลดได้ที่ http://www.geocomp.com.au/geocalc/gcalc420.exe โปรแกรมนี้สามารถคำนวณค่าพิกัดข้าม Datum ไปอีก Datum ได้ (แต่ต้องระวังอีกนิดเพราะโปรแกรม Geocalc มี mistake อยู่ค่า Semi-major axis (a) ของทรงรี Everest 1830 ผิดต้องแก้ไขจาก 6377267.345 เป็น 6377276.345)
  • อีกโปรแกรมคือ GeoTrans ฟรีและ Opensource มี source code ทั้ง Java และ C สามารถ download ได้ที่ http://earth-info.nga.mil/GandG/geotrans/geotrans3.7/master.zip มีให้เลือก download ทั้ง Windows และ Unix (ผมยังไม่เคยนำไปใช้ใน Linux)  GeoTrans เป็นโปรแกรมที่สามารถใช้แปลงค่าพิกัดระหว่าง Datum สนับสนุน Datum & Map Projection ได้หมดทั่วโลก ตัว source code มีคนนำไปแปลงและใช้เช่น แปลงเป็น Delphi ใช้ในโปรแกรม TatukGIS (ไม่ฟรี) ซึ่งเป็น GIS Development Tool Kit มี ActiveX, .NET Winform และ Native VCL ให้เลือกใช้และใช้ได้บนวินโดส์เท่านั้น
  • ส่วนการใช้โปรแกรมแปลงค่าพิกัดที่เป็น Tool ติดมากับโปรแกรมอื่นๆ เช่นบน ERDAS Imagine, Terramodel,  Global Mapper หรือกระทั่งบน Web ก็มีให้เห็นมากพอสมควร

Library สำเร็จรูป “GeographicLib”

  • ถ้าเข้าไปดู source code ของ GeoTrans จะเห็นสูตรที่ใช้แปลงค่าพิกัดเหล่านี้ ยาวเหยียด แยกตามเส้นโครงแผนที่ (Map Projection) ถ้าเกิดจะมีใครพัฒนาโปรแกรมที่เกี่ยวข้อง GIS จาก Ground to top คงต้องศึกษา Library ของ GeoTrans ถ้าไม่คิดจะใช้ Library ของผู้อื่น พูดถึง Library ที่เป็น Opensource ก็ต้อง GeographicLib เขียนโดย Charles Karney ดาวน์โหลดได้ที่ Geographic.zip ตัว GeographicLib เขียนด้วย C++ ตัวนี้ปรับปรุงมาจาก GeoTrans  แต่มีฟังก์ชั่นมากกว่า GeoTrans เสียอีก เช่นการคำนวณเรื่อง Geoid หันมาดู source code ที่เขียนด้วย Delphi ผมหาไม่พบ ต้องยอมรับว่าเปอร์เซนต์ผู้ใช้ Delphi นั้นน้อยน่าใจหาย หลังจาก Borland ล้มหายตายจาก ผมไม่ทราบว่า Developer ของ Delphi จะมีเหลือมากน้อยเท่าใด ส่วน Lazarus ไม่ต้องพูดถึงว่าน้อยขนาดไหน


Datum Transformation Diagram
Datum Transformation Diagram
  • จากไดอะแกรมรูปข้างบนการแปลงค่าพิกัดกริดเช่นจาก UTM Grid (N,E,Z) บน WGS84 => UTM Grid (N,E,Z) บน Indian 1975 datum ได้สองเส้นทางหลัก

    • เส้นทางที่ 1 WGS84( N,E,Z)  ==> WGS84(Lat,Long,Height)  == Molodensky Transformation ==> Indian1975(Lat,Long,Height) ==> Indian1975(N,E,Z)
    • เส้นทางที่ 2 WGS84( N,E,Z)  ==> WGS84(Lat,Long,Height) ==>WGS84(X,Y,Z) ==7,5,3 Parameter Transformation==> Indian1975(X,Y,Z) ==> Indian1975(Lat,Long,Height) ==> Indian1975(N,E,Z)
  • เส้นทางที่ 1 จะสั้นกว่า จากค่าพิกัด Latitude,Long itude คำนวณโดยใช้ Molodensky Transformation หาค่าพิกัด Latitude,Longitude บนอีก datum ได้ ส่วนเส้นทางที่ 2 จะยาวกว่า จากค่าพิกัด Latitude,Longitude จะคำนวณไปหาค่าพิกัดในระบบ Cartesian แล้วใช้ค่าพารามิเตอร์ 3,5 หรือ 7 ค่าในการ Transformation ไปยังระบบ Cartesian บนอีก datum แล้วคำนวณขึ้นไปหา Latitude,Longitude จนถึงค่าพิกัดกริดยูทีเอ็ม ถ้าต้องการ
  • ซึ่งถ้าจะโปรแกรมกันจริงๆ ตามไดอะแกรมด้านบน ก็ใหญ่ขนาด GeographicLib ถ้าจะต้อง support ทุก datum ทุก Map Projection ที่กล่าวไปแล้ว แต่ในตอนนี้ผมขอบีบขอบเขตให้เล็กลง เราจะทำการคำนวณแปลงค่าพิกัดระหว่าง Lat/Long และ UTM Grid บน datum เดียวกันเท่านั้น เช่นคำนวณ Latitude,Longitude ของ WGS 84 ไปเป็นค่าพิกัดยูทีเอ็มกริด (Northing, Easting) บน WGS84 หรือในทางกลับกัน

สิ่งที่ต้องทราบก่อนจะโปรแกรม

  • ขอนำกลับมาที่ datum ก่อนจะไปต่อ datum เป็นพื้นหลักฐานที่แต่ละประเทศต่างก็ใช้แตกต่างกันไป เช่นประเทศในย่านเอเชียตะวันออกเฉียงใต้เช่นพม่า ไทย เขมร ลาว ใช้ Indian 1975 datum มี Spheroid (Ellipsoid) ทรงรีที่เราเรียกว่า Everest 1830 ทรงรีที่เราใช้แทนสัณฐานของโลกจะยุบที่ขั้วโลกเหนือใต้ โป่งออกตรงเส้นศูนย์สูตร สัณฐานของทรงรีจะกำหนดด้วยค่าSemi- major Axis (a) ความยาวของแกนหลัก ค่าการยุบตัว Flattening (f) และความยาวแกนรอง Semi-minor Axis (b) ค่า 3 ค่านี้จะสัมพันธ์กันหมดคือ b = a(1-f) สำหรับ Everest 1830 มีค่า a=6377276.345 , f = 1/300.8017
  • มาถึงยุคของ GPS พื้นหลักฐานของ GPS ในปัจจุบันคือ WGS84 ที่หลายประเทศรวมทั้งไทยเรา หันมาใช้ datum นี้กันมากขึ้น สัณฐานทรงรี WGS94 ค่า a=6378137 ค่า f=1/298.257 223 563

 รูปทรงของทรงรี (ภาพจาก Wikipedia)

  • ทรงรีที่เราใช้กันในโลกนี้มีไม่มาก ไม่น่าปวดหัว มีตัวหลักๆประมาณ 10 กว่าทรงรี
  • ถ้าจะเขียนโปรแกรมให้ใช้ได้ World wide ทั่วโลก สิ่งที่จะปวดหัวก็คือเส้นโครงแผนที่ Map Projection และ Transformation Parameters เพราะมีมากแต่ละประเทศจะใช้ค่าพารามิเตอร์จำนวน 3, 5, 7 เพื่อ Transformation ค่าพิกัดแตกต่างกันไป ตัวอย่างเช่นประเทศแถวตะวันออกกลาง จะใช้ Ellipsoid “International 1924” ใช้เส้นโครงแผนที่ UTM เหมือนกัน แต่จะต่างกันตรงที่ใช้ 3 Parameter Transformation แตกต่างกันไป แต่ถ้าใช้ในประเทศไทยเรา เส้นโครงแผนที่ที่ใช้ก็มีเพียงแต่ Transverse Mercator (TM) ซึ่งเส้นโครงแผนที่ TM ถ้าแบ่งโซนกันแน่นอนก็คือระบบพิกัด Universal Transverse Mercator (UTM) ที่เราคุ้นเคยกันอยู่

เริ่มสร้าง New Project ด้วย Lazarus

  • เปิดโปรแกรม Lazarus คลิกที่เมนู Project > New Project… จะเห็น Dialog เลือก Application คลิกที่ OK
สร้าง New Project
สร้าง New Project
  • จะเห็นฟอร์มเปล่าๆ ดังรูปข้างล่าง ในตอนนี้ผมรัน Lazarus บนวินโดส์
Lazarus รันบนวินโดส์
สร้างโปรเจคใหม่ (Lazarus บนวินโดส์)
  • ก่อนจะไปต่อตอนที่ 2 ผมขอแสดงโปรแกรมขนาดเล็ก ที่เราจะเขียน เพื่อคำนวณค่าพิกัดระหว่าง UTM และ Geographic (Lat/Long) ดังรูปด้านล่าง
โปรแกรมแปลงพิกัดระหว่าง UTM และ Geographic (Lat/Long)
โปรแกรมแปลงพิกัดระหว่าง UTM และ Geographic (Lat/Long)
  • ทิ้งท้ายไว้ตรงนี้ครับ ติดตามตอนต่อไป