Convert Points to Line, Polyline and Polygon

20120105

ไปเก็บข้อมูลโปรไฟล์ถนนวงแหวนตะวันออกมา ข้อมูลได้จาก GPS + IMU
ภายหลังประมวลผลแล้วได้ข้อมูลเป็น จุด ออกมา 164,669 จุด
เมื่อนำไปเปิดใน  QGIS แล้วเล่นทำเอาเครื่องอืดไปเลยทีเดียว
กระนั้นเลยแปลงจาก จุดเรือนแสนเหล่านั้นให้เป็นเส้นดูดีกว่า อาจจะเร็วขึ้นบ้าง
วิธีการนั้นสามารถทำได้หลายวิธี แต่ในที่นี้จะเสนอวิธีที่ใช้ซอฟต์แวร์ QGIS ที่เป็นซอฟต์แวร์ฟรี
ไม่เสียค่าลิขสิทธิ์ และมีความสามารถเพียงพอต่อการใช้งานครับ ดังนี้

QGIS มี Plugins ที่ชื่อ Points2One ที่ทำหน้าที่จัดการเรื่องนี้อยู่ ติดตั้งมันซะ

เมนู Plugins/Fetch Python Plugins…
เลือก Points2One แล้ว คลิก Install plugin

เมื่อติดตั้งเสร็จแล้วจะมี ไอคอนในทูลบาร์เป้นตัวหนังสือสีแดง P2 เมื่อเรียกใช้งานจะเห็นดังภาพ

เลือกรูปแบบว่าต้องการให้เป็น รูปปิด หรือเส้น
กำหนดชื่อ shapefile สำหรับบันทึกผล แล้วคลิก OK

เมื่อเลือกให้ไฟล์ที่แปลงเสร็จแล้วแสดงใน QGIS
ผลที่ได้จะเห็นว่า เส้นอยู่ในแนวเดียวกันกับข้อมูลจุด
เสร็จแล้ว เห็นไหม ไม่ยากเลย

Memory card ติดในเครื่องเอาออกอย่างไรดี

20111224

 

เครื่องโน๊ตบุ๊คสมัยนี้มีช่องสำหรับอ่านการ์ดจากกล้องถ่ายรูป บางทีก็โฆษณาว่าอ่านได้ 4 in 1 หรือ 10 in 1 บางครั้งก็ทำให้ผู้ใช้สับสนว่าอ่านได้หรือไม่ได้ ทีนี้บางทีเกิดพลาดเสียบเข้าไปแล้วตัวแผ่นการ์ดดันไปค้างอยู่ข้างใน ทำยังไงก็ไม่ออก เขี่ยก็แล้วแคะก็แล้วไม่ยอมออกมา ทำไงดีหรือว่าต้องไปเสียตังค์ให้ช่างพันธุ์ทิพย์รื้อเอาออกให้ดี

ผมมีทางออก เรามาดูกัน

การ์ด MS duo pro คาอยู่ในช่องอ่านของโน๊ตบุ๊ค

หากใช้วิธีการอื่นไม่ได้ผล ลองหาอุปกรณ์ดังนี้ครับ สก๊อตเทป หรือมาร์สกิ้งเทป กับคัดเตอร์หรือวัสดุเล็กๆ แบนๆ

อุปกรณ์มี เทป กับคัตเตอร์

ตัดเทปออกมาช่วงหนึ่ง ตัดความกว้างให้เล็กกว่าแผ่นการ์ดที่ติดอยู่ข้างในเล็กน้อย

วางเทปให้เลยช่องไปเล็กน้อยแล้วใช้ใบมีดคัตเตอร์สอดดันส่วนที่เลยให้ตัวเทปแปะติดกับการ์ดเข้าไปข้างในช่อง

 

จากนั้นก็ค่อยๆดึงเทปออกมาอาจต้องดึงแบบยกๆ เล็กน้อยนะครับ ถ้าไม่ออก็ทำใหม่หรือหาเทปกาวที่เหนียวกว่ามาทำครับ

เห็นไหมออกมาแล้ว

หวังว่าคงช่วยได้นะครับ

PostGIS: สำเนาฐานข้อมูลพร้อมแปลงระบบพิกัด

 

20111021

วัตถุประสงค์

• ต้องการทำสำเนาฐานข้อมูลปริภูมิแต่เปลี่ยนระบบพิกัดของข้อมูลปริภูมิ
• ต้องการแปลงระบบพิกัดจาก Geographic/WGS84 ไปเป็น UTM Zone47/WGS84
• กำหนดฐานข้อมูล ชื่อ river มีตาราง river4326 และกำหนดให้ฐานข้อมูลใหม่ให้ชื่อ river32647
• ข้อมูล river4326 ประกอบไปด้วย Column [ID],[name],[length],[basin],[subbasin],[the_geom]

ที่จริงในการเรียกใช้งาน สามารถแปลงระบบพิกัดแบบ on-the-fly ได้เลยและแนะนำให้ทำอย่างนั้นเมื่อใช้งานจริง แต่ในที่นี้ลองว่าถ้าต้องการทำแบบนี้ทำได้อย่างไรบ้าง

วิธีการ

•  SQL Query: SELECT * FROM river4326 WHERE 1=0
  เพื่อ ให้แสดงเฉพาะชื่อ Column จะได้
  | [ID] | [name] | [length] | [basin] | [subbasin] | [the_geom] |
• สร้าง ตารางขึ้นใหม่ในฐานข้อมูลนั้น ให้ชื่อ river_32647 โดยใช้คำสั่ง SQL
  CREATE TABLE river_32647 as SELECT ID, name, length, basin, subbasin, ST_Transform(ST_SetSRID(the_geom, 4326), 32647) as the_geomFROM river4326
  เป็นการสร้างตาราง river_32647  ขึ้นมาใหม่โดยให้มี column เหมือน river4326 และเปลี่ยน the_geom จากระบบพิกัด Geographic ไปเป็น UTM
• เลือกตาราง river_32647 ทำการ backup ไว้ให้ชื่อไฟล์ river_32647.backup
• สร้างฐานข้อมูลให้ เลือกtemplate PostGIS ให้ชื่อฐานข้อมูลใหม่ว่า river32647
• restore ไฟล์ river_32647 เข้ามาในฐานข้อมูล river32647

นี่เป็นเพียงวิธีการหนึ่งที่ทำได้นะครับ ผมเชื่อว่ามีวิธีที่ง่ายกว่านี้แน่ แต่ยังหาไม่เจอ ท่านในมีคำแนะนำช่วยชี้แนะเพิ่มเติมด้วยครับ
อย่างที่บอกถ้าใช้งานจริงแนะนำให้แปลงระบบพิกัด on-the-fly ดีกว่าง่ายกว่าเยอะเลยและเก็บฐานข้อมูลแค่ชุดเดียวด้วย

Libre Office/Open Office: How to check day in month

20111023

When you want to know number of day in month especially number of day in February. it is 28 or 29
if you use formula in cell like “=Day(B3)” when B3 contain 02/02/00 (02 February 2000) you will get “2″.
the February 2000 has 29 day try this formula

“=DAY(EOMONTH(B3,0))” you will get “29″

have fun with Open Office ( it also work in MS Excel)

Libre Office/Open Office: Macro: convert text to string / String to text

20101022

การดึงค่าจากตารางมาสู่ตัวแปรของภาษาเบสิคใน Libre/Open Office มีแค่สองอย่างคือ Value ที่เป็นตัวเลขกับ String ที่เป็นตัวอักษร บางทีได้ค่ามาเป็นตัวอักษรต้องการแปลงเป็นตัวเลขทำได้ง่ายๆ ดังนี้

A = “123.456″
myValue = Val(A)

เมื่อได้ค่ามาเป็นตัวเลขแล้วจะทำการเปลี่ยนเป็นชนิดตัวเลขเป็นอย่างอื่นภายหลังได้

หรือสามารถใช้ฟังก์ชั่นในการเปลี่ยนได้เลยจากค่าในตารางได้ดังนี้
CStr(Var)
converts any data type into a string.
CInt(Var)
converts any data types into an integer value.
CLng(Var)
converts any data types into a long value.
CSng(Var)
converts any data types into a single value.
CDbl(Var)
converts any data types into a double value.
CBool(Var)
converts any data types into a Boolean value.
CDate(Var)
converts any data types into a date value.

อ้างอิง http://wiki.services.openoffice.org/wiki/Documentation/BASIC_Guide/Conversion_Functions_(Runtime_Library)

English: R.I.P.

20111013

มีข่าวว่า Dennis Ritchie ผู้สร้างภาษา C ปูชนียบุคคลอีกหนึ่งคนในวงการคอมพิวเตอร์ ได้เสียชีวิตลงแล้วเมื่อสุดสัปดาห์ที่ผ่านมา
จากนั้นเข้าไปอ่านก็มีหลายคนให้คอมเมนท์ที่มีคำว่า R.I.P. ก็เลยอยากรู้ว่าหมายความว่าอย่างไหร

R.I.P.  ย่อมาจาก Rest In Peace หมายถึงหลับให้สบาย ใช้กับคนตายเท่านั้น

R.I.P. ครับ Dennis Ritchie

ผมอาจจะไม่เคยใช้ภาษา C แต่ก็ใช้ภาษาที่อาศัยโครงสร้างภาษา C เช่น JavaScripts, Python ครับ

อ้างอิง Andrew Biggs Academy

Python: THEOS full scene Pan-Sharpening

20111003

 

ในโครงการ THEOS2Go ส่วนผลิตภัณฑ์เพิ่มค่า Pan-Sharpening หรือ Image Fusion ภาษาไทยเรียกการหลอมภาพ เป็นความสามารถที่ใช้โปรแกรมภาษา Python ในการหลอมภาพโดยอาศัยไลบรารี่ของ Gdal ซึ่งทำได้นานแล้วในโครงการนี้ใช้อัลกอริธึมในการหลอมชื่อ Brovy transform ถ้าถามผมก็เป็นการเทียบบัญยัติไตรยางค์แบบง่ายๆ

newR = (3 * P * R) / (R + G + B + 1.0)
newG = (3 * P * G) / (R + G + B + 1.0)
newB = (3 * P * B) / (R + G + B + 1.0)

where :
newR, newG and newB is new pixel value in Red, Green and Blue band
P =  Panchromatic band
R, G and B = Red, Green and Blue band

เมื่อทำการทดสอบกับภาพเต็มซีนของภาพข้อมูลดาวเทียมธีออส ปรากฎว่าเกิดอาการทำนองตัวแปรของอะเรย์ที่ไช้ในการระบุตำแหน่งของพิกเซลมันล้นเกินความสามารถของอินเด็กซ์ คล้ายๆ กับมันล้นแสดงเออเรอ แล้วก็หยุดทำ แต่ผลของการหลอมก็ยังแสดงเป็นภาพได้แต่ได้ประมาณ 60% ของภาพและเป็นผลที่ไม่ผิดอะไรซูมเข้าไปดูภาพก็หลอมได้ถูกต้อง พยายามแก้โดยการลองดูวิธีการใช้ตัวแปรสำหรับเป็นอินเด็กซ์ในการระบุตำแหน่งพิกเซลของภาพ ด้วยภูมิไม่ถึงจึงไม่สามารถแก้ไขได้
ปรึกษากับคุณ PK ได้คำแนะนำว่าวิธีการเขียนค่าลงไปของผมนั้นใช้วิธีเขียนทีละพิกเซลซึ่งเสียเวลาในการเขียนค่อนข้างมากและอาจเป็นสาเหตุที่ทำให้ตัวแปรของอินเด็กซ์ล้น
อาศัยเอกสารประกอบการสอนของคุณ Chris Garrard จากหน่วยงาน  RS/GS Laboratory ของมหาวิทยาลัย Utah university  ก็ได้รู้ ได้ศึกษาวิธีการอ่านเขียนข้อมูลแบบบล็อค ทำให้ความเร็วของการเขียนมากกว่าเดิมกว่าสองเท่า เยี่ยมเลยแต่ปัญหาหลักก็ยังไม่ถูกแก้ไข ยังคงทำได้ไม่เต็มซีนภาพหลังจากการทำสอบพบว่าสามารถทำได้ประมาณ 7000 x 7000 พิกเซลภาพจึงได้ทำการแบ่งทำทีละ 1/4 ภาพจึงทำให้สามารถเขียนภาพเต็มซีนได้สำเร็จ

ผลการหลอมภาพเต็มซีนของข้อมูลภาพธีออส

ทดสอบกับภาพข้อมูลดาวเทียม SPOT แบบไม่เต็มซีน
ทดสอบหลายครั้งจึงมั่นใจนำไปรวมกับชุดโปรแกรม THEOS2Go ก็สามารถใช้งานได้เขียนลงบน USB drive ก็ได้จึงได้นำไปประกอบการอบรมพอนำไปติดตั้งกับ QGIS 1.7 เครื่องอบรม แล้วทดสอบปรากฎว่า่ python.exe error ซะงั้น เชคแล้วเชคอีกพบว่า ระบบ import gdal ไม่ได้ มึนไปคืนนึงก็จับ gdal ไปยัดใน python ที่คาดว่ามันเรียกใช้ผลปรากฎว่าไช้งานได้ นำไปทดสอบกับคอมพิวเตอร์กรมที่ดินที่ไม่ได้ลง pypthon ไว้สามารถหลอมภาพได้ ทดสอบกับเครื่องของโปรแกรมเมอร์ของ สทอภ. ก็สามารถทำงานได้ก็น่าจะผ่านการทดสอบการใช้งานในเรื่องการทำงานแล้วนะ ส่วนความสวยงามเรื่องนี้คงต้องคุยกันยาว

Libre/Open Office; Macro: Check Sheet Name is Exist?

20110707

“วันนี้คณะวิศวฯ จุฬาฯ รับปริญญา”

การเขียนมาโครใน Libre Office หรือ Open Office มีหลายภาษาในที่นี้ใช้ VBA เพราะเคยทำมาบ้างใน Excel คำสั่งทั่วไปเหมือนกันจะได้ไม่ต้องเรียนรู้ใหม่มาก ที่ต่างกันคือการเรียกฟังก์ชั่นภายใน หรือเรียกคอมโพเนนท์ของ  Libre/Open Office ที่สำคัญตามหายากหน่อยว่าเรียกว่าอะไร พอทำไว้ได้ก็ต้องเอามาแปะไว้ต้องการเมื่อไหร่จะได้มาหา และเป็นการแบ่งปันผู้ที่ยังไม่รู้ ให้เข้าถึงความรู้ได้เร็วขึ้น

Check Sheet Name is Exist?

this script for CALC
ใช้สำหรับโปรแกรม CALS เป็นการตรวจสอบว่า ชีทในชื่อที่กำหนดนี้มีหรือยัง ฟังดูเหมือนง่ายนะครับแค่ดูก็เห็นแล้วว่ามีหรือไม่ แต่ทำยังไงให้เครื่องรู้ลองหลายวิธี หาโค้ทเก่าที่ทำบน Excel มาปรับปรุงสุดท้ายก็ใช้ได้ มาดูกันเลย

อ้างอิงตัวอย่างจากกระทู้นี้ (referent:http://www.ozgrid.com/forum/showthread.php?t=40992&page=1)

Function bWorksheetExists(WSName As String) As Boolean
    On Error Resume Next
    bWorksheetExists = (ThisComponent.Sheets.getByName(WSName).getname = WSName)
End Function

Sub ChectTab()
	inputTab = "Sheet2"'
	If Not bWorksheetExists(inputTab) Then
         'Sheet does not exist
         'do something
         Msgbox(inputTab & " does no exist")
    Else
         'Sheet exists
         Msgbox(inputTab & " exist")
    End If

End Sub

จะนำไปใช้เพิ่มเติมยังไงก็แล้วแต่จะพัฒนาต่อนะครับ ^_^

Python: Least Square Adjustment; Example 5

20110419

find parameter of function y; f(y)

f(y) = a[0] + a[1]*x + a[2]*x*x + a[3]*x*x*x

data is
x = 1, 2, 3, 4, 5  and x is constance
y = 1.74, 2.79, 4.33, 6.16, 8.51

n=5, n[0]=4, r=1

condition equation

y[i] + v[i] = a[0] + a[1]*x[i] + a[2]*x[i]*x[i] + a[3]*x[i]*x[i]*x[i]

then

v[i] = a[0] + a[1]*x[i] + a[2]*x[i]*x[i] + a[3]*x[i]*x[i]*x[i] – y[i]

Solve by Python

## Least Square Adjustment
## Example: 5
## Adjustment by Parameter
## Leveling

''' data x,y fix x
x = 1., 2., 3., 4., 5.
y = 1.74, 2.79, 4.33, 6.16, 8.51

f(y) = a_0 + a_1*x + a_2*x^2 + a_3*x^3
'''

## n= 5, n_0=4, r = n-n_0 = 1
''' condition equation
v_1 = a_0 + a_1*x_1 + a_2*x_1^2 + a_3*x_1^3 - y_1
v_2 = a_0 + a_1*x_2 + a_2*x_2^2 + a_3*x_2^3 - y_2
v_3 = a_0 + a_1*x_3 + a_2*x_3^2 + a_3*x_3^3 - y_3
v_4 = a_0 + a_1*x_4 + a_2*x_4^2 + a_3*x_4^3 - y_4
v_5 = a_0 + a_1*x_5 + a_2*x_5^2 + a_3*x_5^3 - y_5'''

import numpy as np

##set data
x = (1.,2.,3.,4.,5.)
y = (1.74, 2.79, 4.33, 6.16, 8.51)

##general equation for v
##v[i] = a_0 + a_1*x[i] + a_2*x[i]^2 + a_3*x[i]^3 - y[i]
## set data to v + B * Delta = f
## set all a = 1 for find index of a
a = (1.,1.,1.,1.)

B = []
for i in range(5):
 B.append([-a[0],-a[1]*x[i],-a[2]*x[i]*x[i],-a[3]*x[i]*x[i]*x[i]])

''' from v + B * Delta = f
y[i] + v[i] - a[0] - a[1]*x[i] - a[2]*x[i]*x[i] - a[3]*x[i]*x[i]*x[i] = 0
v[i] - a[0] - a[1]*x[i] - a[2]*x[i]*x[i] - a[3]*x[i]*x[i]*x[i] = -y[i]
v[i] + [-a[0] - a[1]*x[i] - a[2]*x[i]*x[i] - a[3]*x[i]*x[i]*x[i]] = -y[i]
'''
f = np.asmatrix(np.dot(y,-1)).reshape(len(y),1)

## weight matrix
W = np.matrix([[1./np.square(0.02),0.,0.,0.,0.],[0,1./np.square(0.02),0,0,0],[0,0,1./np.square(0.02),0,0],[0,0,0,1./np.square(0.04),0],[0,0,0,0,1./np.square(0.04)]])

B = np.asmatrix(B)
## N = B^T * W * B
N = B.T * W * B

## t = B^T * W * f
t = B.T * W * f

## delta = N^(-1) * t
invN = np.linalg.inv(N)
delta = invN * t

print "parameter"
print delta

result is
parameter
[[ 1.13834711]
 [ 0.35227273]
 [ 0.25132231]
 [-0.00549587]]

parameter is a[0], a[1], a[2] and a[3]

parameter
a[0] = 1.13834711
a[1] = 0.35227273
a[2] = 0.25132231
a[3] = -0.00549587
Ans

Python: Least Square Adjustment using parameter

20110418

Least square adjustment

Example 4 : Adjustment by parameter

Leveling problem as show in image

example4: Least square adjustment

the different of height show in table below

–from– —– — to — — Δh (m.) –
======|=======|==========
E                     A              42.107
B                     A              12.424
C                     B               42.251
C                     D                 8.464
D                     E                 4.138
D                     A              46.269
D                     B               33.802
======|=======|==========

Elevation of point E is 500.000 m. and weight of observation data is invert distance of each point. find all point elevation

solve:
n = 7, n_0 = 4, r = 7-4=3, u = n_0 = 4

set data to equation
v + BΔ = Wf

from least square let
N = B(T) * W * B
t = B(T) * W * f
and
Δ =inv(N) * t

calculate in Python

## Least Square Adjustment
## Example: 4
## Adjustment by Parameter
## Leveling

## Point E has elevation is 500.00 m. MSL
## n = 7, n0 = 4, r = 3

import numpy as np
## V + B * Delta = d - l * W

##Delta

Delta = np.array([["A"],["B"],["C"],["D"]])

##Matrix B

B = np.matrix([[-1.,0,0,0],[-1,1,0,0],[0,-1,1,0],[0,0,1.,-1],[0,0,0,1],[-1,0,0,1],[0,-1,0,1]])

print B

##Matrix of Weight

W = np.matrix([[1./8,0.,0.,0.,0.,0.,0.],[0,1./6,0.,0.,0.,0.,0.],[0.,0.,1./3,0.,0.,0.,0.],[0,0,0,1./6,0,0,0],[0,0,0,0,1./7,0,0],[0,0,0,0,0,1./3,0],[0,0,0,0,0,0,1./5]])

print W

##Matrix f

f = np.matrix([[-542.107],[-12.424],[-42.251],[-8.464],[495.862],[-46.269],[-33.802]])

print f

## N = B^T * W * B

N = B.T * W * B

## t = B^T * W * f

t = B.T * W * f

## delta = N^(-1) * t
invN = np.linalg.inv(N)
delta = invN * t

print "Leveling"
print delta

and output show this
Leveling
[[ 542.11261789]
 [ 529.66589837]
 [ 487.40762703]
 [ 495.85708435]]

เริ่มจากง่ายๆ ไปก่อนนะครับจะได้เรียนรู้การใช้ array ต่างๆ ของไพธอนและการทำงานของแมทริกซ์ใน numpy
คำตอบที่ได้คือค่าระดับของจุด A,B,C และ D ตามลำดับ