นำเข้าขอบเขตแปลงที่ดินจากระบบแผนที่ราชการเข้า SketchUp

พอดีผมต้องขึ้นโมเดลอาคารบนที่ดินแปลงหนึ่ง ใช้ Add Location ของ SketchUp ดึงภาพถ่ายดาวเทียมกับระดับพื้นดินมาแล้วเรียบร้อย แต่สิ่งที่ขาดคือเส้นแบ่งแปลง — ภาพถ่ายมันเห็นแค่แนวรั้วกับแนวต้นไม้ ซึ่งไม่ตรงกับขอบเขตจริงเสมอไป จะวางอาคารให้ถอยร่นตามระยะก็กะไม่ถูก

ข้อมูลขอบเขตแปลงมีอยู่ในระบบแผนที่แปลงที่ดินของราชการอยู่แล้ว และเบื้องหลังเว็บมันเปิดเป็น WFS ซึ่งเป็นมาตรฐานเปิดของ OGC ตัวเดียวกับที่ระบบ GIS ทั่วโลกใช้ ผมเลยลองดึงมาวางในโมเดลดู ปรากฏว่าทำได้ครับ ตรงกับภาพถ่ายพอดี บทความนี้เป็นบันทึกทางเดินทั้งหมด

เข้าระบบตามปกติก่อน แล้วค่อยเรียก WFS

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

กับดักใหญ่ที่สุด: ชั้นข้อมูลนี้ไม่รับ bbox เปล่าๆ

ผมเริ่มแบบที่ทุกคนเริ่ม คือยิง WFS ขอฟีเจอร์ในกรอบพิกัดที่สนใจ ตามตำราเลย ผลที่ได้คือ ศูนย์ฟีเจอร์ ไม่มี error ไม่มีอะไร ตอบมาเป็น FeatureCollection ว่างเปล่าสวยงาม ทำให้หลงคิดว่าดึงผิดชั้นหรือพิกัดผิด นั่งไล่อยู่นานพอสมควร

สาเหตุคือชั้นข้อมูลแปลงที่ดินไม่ได้เป็นตารางตรงๆ แต่เป็น parameterized view — คือมุมมองที่ต้องส่งค่าพารามิเตอร์เข้าไปด้วย ฝั่ง GeoServer ใช้ตัวแปร viewparams ถ้าไม่ส่งไป เงื่อนไขในคิวรีข้างในจะไม่ตรงอะไรเลย ก็เลยได้ผลลัพธ์ว่างทุกครั้งแม้จะระบุกรอบพิกัดถูกต้อง

ค่าที่ต้องส่งคือ รหัสระวางแผนที่ ที่ครอบพื้นที่นั้น:

https://<host>/geoserver/<workspace>/wfs
  ?service=WFS&version=1.0.0&request=GetFeature
  &typeName=<workspace>:V_PARCEL47
  &viewparams=utmmap:503624412
  &outputFormat=application/json
  &srsName=EPSG:4326

รหัสระวางประกอบจากป้ายที่เห็นบนแผนที่ โดยแปลงเลขโรมันตรงกลางเป็นตัวเลข:

ROMAN = {"I": "1", "II": "2", "III": "3", "IV": "4"}

def sheet_code(label: str) -> str:
    """'5036 II 4412'  ->  '503624412'"""
    a, roman, b = label.split()
    return f"{a}{ROMAN[roman]}{b}"

print(sheet_code("5036 II 4412"))   # 503624412  (ค่าตัวอย่าง)

แล้วจะรู้ได้ยังไงว่าพื้นที่ที่สนใจอยู่ในระวางไหน — ถามจากชั้นดัชนีระวางครับ ชั้นดัชนีตัวนี้รับ bbox ได้ตามปกติ ไม่ต้องมี viewparams ส่งกรอบพิกัดของโมเดลไป มันคืนรายการระวางที่ทับกรอบนั้นมาให้ พร้อมป้ายชื่อในคุณสมบัติของแต่ละฟีเจอร์

ข้อควรระวังเรื่องลำดับพิกัดใน bbox — WFS เวอร์ชัน 1.0.0 เรียงเป็น lon,lat (x ก่อน) ส่วนเวอร์ชัน 1.1.0 ขึ้นไปกลับด้านเป็น lat,lon ตามที่ EPSG นิยามไว้ ถ้าใส่สลับ ผลลัพธ์คือได้ค่าว่างอีกเหมือนกัน แล้วจะไปโทษ viewparams ผิดตัว ผมล็อกไว้ที่ 1.0.0 แล้วเรียง lon,lat ตลอด จะได้ไม่ต้องมานั่งคิดทุกครั้ง

แปลงพิกัดโดยไม่ต้องลงไลบรารีเพิ่ม

เรขาคณิตต้นทางเก็บเป็น UTM โซน 47 หน่วยเมตร ซึ่งเอาไปใช้กับ SketchUp ตรงๆ ไม่ได้ เพราะฝั่งโมเดลคิดเป็นละติจูด/ลองจิจูด ทางที่คนทั่วไปทำคือลง pyproj แล้วแปลงเอง แต่เครื่องผมตอนนั้นไม่มี pyproj และไม่อยากลงเพิ่มกลางงาน

ทางลัดคือให้ฝั่งเซิร์ฟเวอร์แปลงให้ ใส่ srsName=EPSG:4326 เข้าไปในคำขอ GeoServer จะรีโปรเจกต์ให้เสร็จก่อนส่งกลับ ได้ GeoJSON ที่เป็น lon/lat มาเลย ไม่ต้องแตะเรื่องระบบพิกัดเองสักบรรทัด

โซน 47 ครอบประเทศไทยส่วนใหญ่ ถ้างานอยู่แถบตะวันออกสุดต้องใช้ชั้นข้อมูลของโซน 48 แทน ตรงนี้ดูจากลองจิจูดของพื้นที่งานได้เลย เส้นแบ่งอยู่ที่ 102 องศาตะวันออก

ฝั่ง SketchUp: latlong_to_point คือทั้งหมด

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

สคริปต์ฝั่ง SketchUp ทั้งหมดมีเท่านี้ครับ:

require 'json'

model = Sketchup.active_model
ents  = model.active_entities

# ไฟล์นี้เตรียมไว้ก่อนแล้วจากฝั่ง Python: [[[lon,lat],[lon,lat],...], ...]
parcels = JSON.parse(File.read('D:/work/parcels.json'))

model.start_operation('Import parcels', true)

parcels.each do |ring|
  pts = ring.map do |lon, lat|
    p = model.latlong_to_point([lon, lat])
    Geom::Point3d.new(p.x, p.y, p.z + 5.cm)   # ยกขึ้นกันเส้นจมภาพถ่าย
  end
  pts << pts.first unless pts.first == pts.last

  created = ents.add_edges(pts)

  # วงปิดที่อยู่ระนาบเดียวกันจะเกิดหน้าขาวขึ้นมาเองและบังภาพดาวเทียม -> ลบทิ้ง
  created.each do |edge|
    next unless edge.valid?
    edge.faces.each { |f| f.erase! if f.valid? }
  end
end

model.commit_operation

มีสองบรรทัดที่แลกมาด้วยการลองผิดหลายรอบ

บรรทัดยก z ขึ้น 5 เซนติเมตร — ถ้าวางเส้นไว้ที่ระดับศูนย์พอดี มันจะอยู่ระนาบเดียวกับภาพถ่ายดาวเทียม แล้วเกิดอาการ z-fighting คือเส้นกะพริบเข้าๆ ออกๆ ตอนหมุนมุมมอง ยกขึ้นนิดเดียวก็หายสนิท

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

อีกเรื่องที่ควรรู้ไว้ — Add Location มีตัวเลือกดึงกลุ่มอาคารสามมิติมาด้วย ซึ่งจะได้กล่องสีขาวๆ โผล่มาในโมเดล อันนั้นเป็นอาคารจากชุดข้อมูลของบริการแผนที่ ไม่ใช่เส้นแปลงที่เราเพิ่งนำเข้า อย่าเผลอลบผิดตัว

ตัดขอบให้พอดีกรอบก่อนวาด

ระวางหนึ่งใบมีแปลงเป็นร้อย แต่โมเดลครอบแค่หย่อมเดียว ผมเลยตัดทิ้งตั้งแต่ฝั่ง Python ก่อนส่งเข้า SketchUp ไม่งั้นโมเดลจะหนักโดยไม่จำเป็น:

def clip(rings, bbox):
    """เก็บเฉพาะแปลงที่มีจุดใดจุดหนึ่งตกในกรอบโมเดล"""
    lon0, lat0, lon1, lat1 = bbox
    out = []
    for ring in rings:
        if any(lon0 <= lon <= lon1 and lat0 <= lat <= lat1 for lon, lat in ring):
            out.append(ring)
    return out

# กรอบเอามาจาก point_to_latlong ของมุมโมเดลใน SketchUp (ค่าตัวอย่าง)
bbox = (99.7980, 16.1960, 99.8060, 16.2030)

ผมเลือกเงื่อนไข "มีจุดใดจุดหนึ่งตกในกรอบ" แทนที่จะเป็น "ทุกจุดต้องอยู่ในกรอบ" เพราะแปลงที่คร่อมขอบโมเดลก็ยังอยากเห็น จะได้รู้ว่าเพื่อนบ้านอยู่ตรงไหน

ทำไมไม่ใช้ปลายทางที่คลิกทีเดียวได้เลย

ระบบมีปลายทางแบบคลิกบนแผนที่แล้วได้ข้อมูลแปลงกลับมา ตอนแรกผมนึกว่าจะจบง่ายๆ ตรงนั้น แต่พอไล่ดูจริงมันคนละงานกัน:

สรุปคือใช้คู่กันครับ ดึงรูปร่างจาก WFS แล้วถ้าอยากได้เลขที่แปลงกับเนื้อที่มาแปะเป็นข้อความในโมเดล ค่อยไปถามอีกปลายทางเป็นรายแปลง

สรุป

ทางเดินทั้งหมดสั้นกว่าที่คิดครับ — เอากรอบพิกัดจากโมเดลไปถามชั้นดัชนีว่าอยู่ระวางไหน เอารหัสระวางไปขอชั้นแปลงพร้อม viewparams และ srsName=EPSG:4326 ตัดขอบให้เหลือเฉพาะที่ต้องใช้ แล้ววาดด้วย latlong_to_point ปิดท้ายด้วยยก z ขึ้นเล็กน้อยและลบหน้าที่เกิดเอง

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

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