ค้นหาข้อมูลในบล็อก (Search in this blog)

Wednesday, September 30, 2009

gdalwarp แบบกำหนดค่าพิกัดจุดควบคุม (ของฟรีและดีก็มีในโลก)

วันนี้ขอเสนอ การใช้งาน gdalwarp แบบใช้ค่าพิกัดจุดควบคุม ตัวอย่างของการทำงานนี้คือ วิธีการแปลงระบบพิกัดจากข้อมูลภาพที่ไม่มีระบบพิกัดให้มีระบบพิกัดตามต้องการ เช่น แปลงข้อมูลภาพแผนที่ภูมิประเทศที่ผ่านการสแกนให้มีระบบพิกัด หรือข้อมูลภาพถ่ายดาวเทียมที่ไม่มีค่าพิกัดให้มีค่าพิกัด เป็นต้น
ขั้นตอนการทำงานโดยย่อมีดังนี้ครับ

1. กำหนดค่าพิกัดของจุดควบคุม
ค่าพิกัดของจุดควบคุมนั้นเป็นค่าที่ได้จาก การรังวัดบนภาพและการรังวัดค่าพิกัดบนพื้นดิน (อาจจะใช้การวัดจากแหล่งข้อมูลอื่น เช่น วัดจากข้อมูลภาพถ่ายทางอากาศ) รูปแบบข้อมูลที่จัดเก็บมีดังนี้ครับ
Image Coordinate (x, y) และ Ground Coordinate (E, N) เช่น 1 1 98.85231 18.93516 ซึ่งหมายถึง ข้อมูลภาพที่มีพิกัด (1,1) นั้นจะมีพิกัดพื้นดิน (98.85231,18.93516 ) เป็นต้น สำหรับการทดลองครั้งนี้ผมได้ลองใช้จุดควบคุมทั้งสิ้น 4 จุดดังนี้
จุดที่ 1 มีค่าพิกัด 1 1 98.85231 18.93516
จุดที่ 2 มีค่าพิกัด 12000 1 99.10726 18.89724
จุดที่ 3 มีค่าพิกัด 12000 12000 99.05081 18.67978
จุดที่ 4 มีค่าพิกัด 1 12000 98.79618 18.71761
2. เพิ่มค่าพิกัดของจุดควบคุมใส่ข้อมูลภาพ
เป็นการเพิ่มข้อมูลเกี่ยวกับจุดควบคุมเข้าสู่ข้อมูลภาพ การทำงานขั้นตอนนี้สามารถทำได้โดยการใช้งานผ่าน gdal_translate ดังตัวอย่างดังนี้
gdal_translate -gcp [x y E N] ... input_file output_file
ตัวอย่งเช่น
>>gdal_translate -gcp 1 1 98.85231 18.93516 -gcp 12000 1 99.10726 18.89724 -gcp 12000 12000 99.05081 18.67978 -gcp 1 12000 98.79618 18.71761 sample.tif sample_gcp.tif
3. การแปลงระบบพิกัด
เป็นขั้นตอนที่ใช้ข้อมูล(ค่าพิกัด)ที่อยู่ในข้อมูลภาพมาทำการแปลงระบบพิกัดการใช้งานสามารถทำงานได้ผ่านโปรแกรม gdalwarp ดังตัวอย่างดังนี้
>>gdawarp -r cubic -t_srs epsg:4326 sample_gcp.tif sample4326.tif

ข้อพึงระวัง
  • โปรแกรมยังไม่สามารถแสดงเศษเหลือ (Residual) ของจุดควบคุมได้ทำให้เราไม่ทราบว่าจุดควบคุมเหล่านั้นมีความคลาดเคลื่อนมากน้อยเพียงใด

Friday, September 25, 2009

การแก้ปัญหาระบบสมการเชิงเส้นด้วย Octave และตัวอย่างสำหรับการแปลงแบบสัมพรรค (ของฟรีและดีก็มีในโลก)

วันนี้ขอแนะนำวิธีการแก้ปัญหาระบบสมการแบบใช้เมทริกซ์เข้ามาช่วย (ลืมเรื่องสมการเชิงเส้นและเมทริกซ์กันรึยังครับ แฮะๆ) ตัวอย่างผมมีระบบสมการเชิงเส้นแบบสองตัวแปร (ตัวไม่รู้ค่า 2 ตัว) ดังนี้
2x + y = 19 ----- (1)
x + y = 14 ----- (2)
ถ้าเป็นสมัยมัธยมต้นเราก็จะใช้วิธีย้ายข้างสมการ ดังนี้ครับ
เริ่มต้นจากย้ายข้าง สมการที่ (1) จะได้ y = 19 + 2x แล้วแทนค่า y ในสมการที่ (2) จะได้ค่า x เท่ากับ 5
หลังจากนั้นให้แทนค่า x ในสมการที่ (1) หรือ (2) ก็ได้แล้วจะได้ค่า y = 9
ดังนั้นจากระบบสมการดังกล่าวเราจะได้ค่า x และ y เท่ากับ 5 และ 9 ตามลำดับ

จากตัวอย่างดังกล่าวจะเห็นได้ว่าขั้นตอนไม่ค่อยซับซ้อนเท่าไหร่นะครับ คราวนี้เรามาดูตัวอย่างที่มีการใช้งานจริงในด้านภูมิสารสนเทศนะครับ
สมมุติผมต้องการแปลงระบบพิกัดของข้อมูลภาพ (Image Georeferencing) โดยอาศัยการแปลงแบบสัมพรรคแบบสองมิติ (2D Affine Transformation) ซึ่งก็คือสมการโพลีโนเมียลกำลังหนึ่ง (1st order polynomial) นั่นเอง ก่อนอื่นผมต้องคำนวณก่อนว่าหากใช้สมการโพลีโนเมียลกำลังหนึ่งนั้นต้องรังวัดจุดควบคุมภาคพื้นดินและรังวัดบนภาพทั้งหมดกี่จุด (วิธีการคำนวณดูได้จากกระทู้ก่อนหน้านี้ ที่่นี่ ) เพื่อความสะดวกและง่ายต่อการอธิบายผมขออนุญาตใช้จำนวนสมการเท่ากับจำนวนตัวไม่รู้ค่า (Unknown) หรือแค่ให้พอหาคำตอบได้เท่านั้นแต่ไม่เหมาะกับการทำงานจริง (การทำงานจริงควรรังวัดให้เกินพอนะครับ จำนวนจุดยิ่งเยอะและค่า residual น้อยๆยิ่งดีครับ)
สรุป นะครับ หากต้องการแปลงระบบพิกัดด้วยการแปลงสัมพรรคแบบสองมิติ นั้นต้องใช้ จุดควบคุม 3 จุด (ขั้นต่ำ)
สำหรับครั้งนี้ผมสมมุติว่าผมรังวัดค่าพิกัดภาพ (Image Coordinate) และพิกัดพื้นดิน (Ground Coornidate) ได้ดังนี้ครับ

โดย (ix,iy) คือพิกัดที่วัดได้บนภาพหรือพิกัดภาพนั่นเอง ส่วน (gX,gY) คือค่าพิกัดพื้นดิน (ระบบใดๆก็ได้ ตัวอย่างนี้ผมขอใช้ระบบพิกัดภาพเทียบกับระบบพิกัดกริด UTM 47N) และจากระบบสมการโพลีโนเมียลกำลังหนึ่งหรือสมการการแปลงแบบสัมพรรคซึ่งมีรูปสมการดังนี้

Xi = a (xi) + b (yi) + c
Yi = d (xi) + e (yi) + f

จากสมข้างต้นจะเห็นได้ว่าถ้าเราต้องการทราบพิกัดของจุดภาพใดๆ(จุดใดก็ได้บนภาพ) ว่ามีค่าพิกัดพื้นดิน (UTM 47N) เป็นเท่าไหร่ เราก็สามารถทำได้เพียงใส่ค่าพิกัดภาพ (xi,yi) แทนลงไปในสมการก็จะได้ค่าพิกัดพื้นดิน (Xi,Yi) ทันที
ง่ายไหมครับ ดูผ่านๆหมือนง่ายนะครับ แต่จากสมการข้างต้นสิ่งที่เรายังไม่รู้คือค่าของสัมประสิทธิ์ (Coeficient) ต่างๆ ซึ่งก็คือค่า a,b,c,d,e,f นั่นเอง แล้วเราจะหาได้อย่างไรละครับ (ค่าสัมประสิทธิ์เหล่านี้จะได้จากการหาความสัมพันธ์ระหว่างระบบพิกัดทั้งของโดยใช้ค่าพิกัดของจุดร่วมหรือ Common points ที่เรารังวัดได้จากทั้งระบบพิกัดภาพและระบบพิกัดพื้นดิน)

วิธีการก็คือเราสามารถหาได้โดยอาศัยความสัมพันธ์ของพิกัดจุดภาพและพิกัดพื้นดิน (จากรังวัดจุด) จากตารางตัวอย่างข้างบนครับ ขั้นตอนต่อไป ผมจะทำการแทนค่าพิกัดของจุดร่วมที่รังวัดระบบบนภาพและระบบพิกัพื้นดิน ซึ่งจะได้ระบบสมการใหม่ดังนี้ (วัด 3 จุดจะได้ 6 สมการ)
499271.67=a(8445.75)+b(6611.75)+c ----(1)
2074132.01=d(8445.75)+e(6611.75)+f ----(2)
487786.57=a(3436.19)+b(8011.69)+c ----(3)
2073436.38=d(3436.19)+e(8011.69)+f ----(4)
497860.20=a(8433.13)+b(9457.38)+c ----(5)
2068443.25=d(8433.13)+e(9457.38)+f ----(6)

ระบบสมการข้างต้นจะเห็นว่าจำนวนสมการ (6 สมการ) เท่ากับจำนวนตัวไม่รู้ค่า (6 ตัวคือ a,b,c,d,e,f)
ดังนั้นสามารถหาคำตอบ!!!
แต่............................. ถ้าจะให้มานั่งย้ายข้างสมการไปมาหรือใช้ความรู้ชั้นมัธยมต้นคงจะไม่ไหวแล้วละครับ เพราะว่าตัวเลขมากขึ้นจำนวนสมการก็มากขึ้นตามไปด้วย แล้วเราจะทำยังไงดี ???

วันนี้ผมมีตัวอย่างวิธีการแก้ปัญหาระบบสมการโดยอาศัยเมทริกซ์เข้ามาช่วยมานำเสนอครับ
จากระบบสมการก่อหน้านี้ผมขอจัดระบบสมการใหม่อีกทีและแทนในเมทริกซ์ต่างๆในรูป AX = B ดังนี้ครับ
เมทริกซ์ A บรรจุสัมประสิทธิ์ของตัวไม่รู้ค่า มีขนาด 6x6
เมทริกซ์ X บรรจุตัวไม่รู้ค่า มีขนาด 6x1
เมทริกซ์ B บรรจุค่าพิกัดพื้นดิน มีขนาด 6x1
หลังจากนั้นเราจะใช้คุณสมบัติของเมทริกซ์ ซึ่งไม่สามารถย้ายข้างสมการแบบสมการทั่วๆไปได้ เช่น ย้ายจากฝั่งซ้าย(การคูณ)ไปฝั่งขวา(กลายเป็นการหาร) เป็นต้น
คุณสมบัติของเมทริกซ์นั้นจะย้ายข้างสมการได้แต่เปลี่ยนจากคูณฝั่งซ้ายไปเป็นการคูณด้วยอินเวอร์สที่ฝั่งขวา เช่นจากตัวอย่างข้างต้น AX=B จะเป็น X= Inverse (A)*B ซึ่งความยากคือการหาอินเวอร์สนี่แหละครับ (ซึ่งมีหลายวิธี เช่น Gauss-Jordan เป็นต้น) ครั้งนี้ผมขออนุญาตให้ทุกท่านใช้เครื่องมือฟรี (โปรแกรม) ที่นิยมอย่างมากในแวดวง Open Source ซึ่งคือเจ้า Octave นั่นเอง ซึ่งก็โหลดได้จาก ที่นี่ ครับ (โหลดแล้วติดตั้งกันเลยนะครับ จะได้ทดลองทำตามกันเลย)
หลังจากติดตั้งแล้วให้ทำตามขั้นตอนดังนี้
1) ป้อนค่าเมทริกซ์ A ดังนี้

>>a=[8445.75 6611.75 1 0 0 0;0 0 0 8445.75 6611.75 1;3436.19 8011.69 1 0 0 0;0 0 0 3436.19 8011.69 1;8433.13 9457.38 1 0 0 0;0 0 0 8433.13 9457.38 1]

### สำหรับการป้อนค่าสมาชิกของเมทริกซ์ทนั้นให้ป้อนค่าตัวเลขที่ละแถวซึ่งจะแบ่งค่าของแต่ละคอลัมน์ด้วยการเว้นว่าง (Space) และเมื่อป้อนค่าสมาชิกในแต่ละแถวแล้วการขึ้บรรทัดใหม่นั้นให้ใช้เครื่องหมาย ;

2) ป้อนค่าของเมทริกซ์ B ดังนี้

>>b=[499271.67;2074132.01;487786.57;2073436.38;497860.20;2068443.25]

3) ทำการคำนวณค่าสัมประสิทธ์จากสมการ X=Inverse (A) * (B) โดยพิมพ์คำสั่งในโปรแกรมดังนี้
>>inv (a) * b

4)ผลลัพธ์ที่ได้หมายถึงค่าของสัมประสิทธิ์ของแต่ละตัว (a,b,c,d,e,f) ตามลำดับ

จากตัวอย่างการคำนวณทำให้เราได้ความสัมพันธ์ระหว่างพิกัดภาพและพิกัดพื้นดินดังนี้

Xi = 2.22330667 (xi) - 0.71668202 (yi) + 485888.70006051
Yi = -0.42032228 (xi) - 2.00098553 (yi) + 2090911.96301988

นั่นหมายความว่า จุดภาพที่ (1,1) นั้นค่าพิกัดพื้นดินเท่ากับ (485890.21,2090909.54) และจุดภาพที่ (1,2) มีค่าพิกัดพื้นดินเท่ากับ (485889.49,2090907.54) ... จนถึงจุดภาพสุดท้าย

Friday, September 18, 2009

gdaldem เครื่องมือสำหรับวิเคราะห์และการมองเห็นข้อมูล DEM (ของฟรีและดีก็มีในโลก)

วันนี้ขอแนะนำเครื่องมือสำหรับการวิเคราะห์ (Analyze) และการมองเห็น (Visualize) สำหรับข้อมูลแบบจำลองระดับสูง (DEM) ในตอนแรกนี้จะนำเสนอการคำนวณ ความลาดเอียง (Slope) อย่างง่ายโดยขั้นตอนมีดังนี้ครับ
ก่อนอื่นเราต้องเปลี่ยนที่อยู่ไปยังที่เก็บไฟล์แบบจำลองระดับสูง เช่น
>>cd c:\sample
และใช้คำสั่ง gdaldem ซึ่งมี syntax ดังนี้ครับ

gdaldem -option input_file out_file

ตัวอย่างการคำนวณความลาดเอียงมีดังนี้ครับ

>>gdaldem slope srtm32647.tif slope.tif

หากต้องการให้ผลการคำนวณอยู่ในแบบเปอร์เซนต์ก็ให้เพิ่ม option -p ดังนี้ครับ (ค่าเริ่มต้นหรือค่า default คือหน่วยองศาครับ)

>>gdaldem slope -p srtm32647.tif slope_p.tif

Wednesday, September 9, 2009

หากต้องการแปลงพิกัดด้วยโพลีโนเมียลแล้วต้องใช้จุดควบคุมจำนวนเท่าไหร่???

หลายท่านคงเคยทำการแปลงข้อมูลภาพที่ยังไม่มีระบบพิกัด(โลก)ให้มีระบบพิกัดที่ต้องการ เช่น ระบบพิกัดภูมิศาสตร์หรือระบบพิกัดกริดยูทีเอ็ม หรือที่เรียกกันว่า Image registration หรือ Image warping กระบวนการนี้มีวัตถุประสงค์หลักคือ การทำให้ข้อมูลภาพนั้นมีระบบพิกัดที่ต้องการ ขั้นตอนการทำงานนั้นสามารถทำได้หลากหลายวิธี ซึ่งแต่ละวิธีนั้นมีข้อมูลหรือสิ่งที่ต้องการแตกต่างกัน (ในที่นี้เราจะพูดถึงเฉพาะการแปลงแบบ 2 มิติเท่านั้นนะครับ)
จุดควบคุมภาคพื้นดิน (Ground Control Point - GCP) นั้นจำเป็นอย่างยิ่งสำหรับกระบวนการนี้ เนื่องจากใช้สำหรับอ้างอิงหรือเป็นตัวกำหนดความสัมพันธ์ระหว่างสองระบบพิกัดใดๆ วิธีการได้มาของจุดควบคุมภาคพื้นดินนั้นสามารถทำได้หลายทางเช่น อ่านค่าพิกัดจากแผนที่ภูมิประเทศ แผนที่ภาพหรือการสำรวจภาคสนาม เป็นต้น
หลังจากได้ข้อมูลจุดควบคุมภาคพื้นแล้วสิ่งที่ต้องทำในขั้นต่อไปคือ การเลือกวิธีการแปลงพิกัด ซึ่งมีให้เลือกมากมายเช่น การแปลงแบบ Helmert, Affine หรือ Polynomial เป็นต้น ในส่วนของวิธีการแปลงและการคำนวนค่าของสัมประสิทธิ์ขอข้ามไปก่อนนะครับเนื่องจากมีรายละเอียดค่อนข้างเยอะ (ขี้เกียจพิมพ์นะครับ แฮะๆๆ)
วิธีการที่นิยมใช้กันอย่างแพร่หลายคือการแปลงโดยอาศัยโพลีโนเมียลอันดับต่างๆ ซึ่งแต่ละวิธีนั้นจะให้ความถูกต้อง ความน่าเชื่อถือรวมถึงความยุ่งยากในการคำนวนที่แตกต่างกัน
วันนี้ขอนำเสนอการคำนวนจำนวนจุดควบคุมภาคพื้นดินที่ต้องการ(อย่างน้อย)สำหรับโพลีโนเมียลอันดับต่างๆ การคำนวนหาจำนวนจุดควบคุมที่ต้องการสำหรับโพลีโนเมียลอันดับต่างๆสามารถคำนวนหาได้จากสูตรดังนี้

[(n+1)*(n+2)]/2

โดย n คือ จำนวนอันดับของโพลีโนเมียลที่ต้องการ เช่น หากต้องการแปลงระบบพิกัดด้วยโพลีโนเมียลอับดับที่ 2 (2nd order) จะใช้จำนวนจุดควบคุมอย่างน้อย [(2+1)*(2+2)]/2 หรือเท่ากับ 6 จุด เป็นต้น

ในทางปฏิบัติเรามักจะรังวัดจำนวนจุดควบคุมให้เกินพอเพื่อให้เข้าสู่ระบบ Overdetermined System แล้วจึงเลือกเอาเฉพาะจุดที่มีความถูกต้องและน่าเชื่อถือ (หรือขจัดจุดที่มีความคลาดเคลื่อนสูงออกจากการคำนวน) โดยอาศัยการปรับแก้แบบกำลังสองน้อยที่สุด (Least Square Adjustment)
ข้อดีของการเพิ่มจำนวนจุดควบคุมให้มากกว่าจำนวนจุดที่ต้องการหรือทำให้เข้าสู่ระบบ Overdetermined System คือ เราสามารถขจัดความคลาดเคลื่อนที่แฝงอยู่ในข้อมูลได้ เช่น อาจมีการรังวัดหรือป้อนค่าพิกัดของจุดควบคุมผิดพลาด (Mistake) หรือจุดควบคุมที่เลือกใช้มีความคลาดเคลื่อนสูงเกินกว่าจะยอมรับได้ (ใช้เป็นการตรวจสอบคุณภาพของจุดควบคุมเบื้องต้นได้ด้วยนะ เฮ่อๆๆ) เป็นต้น
ส่วนบางท่านที่ไม่ต้องการจำสูตรก็สามารถดูจำนวนจุดที่ต้องการสำหรับการแปลงระบบพิกัดด้วยโพลีโนเมียลอันดับต่างๆได้ดังนี้ครับ
1st ต้องการจุดควบคุมอย่างน้อย 3 จุด แนะนำให้ใช้อย่างน้อย 4 จุดสำหรับการทำงานจริง (หรือการแปลงแบบ Affine)
2nd ต้องการจุดควบคุมอย่างน้อย 6 จุด แนะนำให้ใช้อย่างน้อย 7 จุดสำหรับการทำงานจริง
3rd ต้องการจุดควบคุมอย่างน้อย 10 จุด แนะนำให้ใช้อย่างน้อย 11 จุดสำหรับการทำงานจริง
4th ต้องการจุดควบคุมอย่างน้อย 15 จุด แนะนำให้ใช้อย่างน้อย 16 จุดสำหรับการทำงานจริง
5th ต้องการจุดควบคุมอย่างน้อย 21 จุด แนะนำให้ใช้อย่างน้อย 22 จุดสำหรับการทำงานจริง