เริ่มต้นใช้งาน GDAL/OGR

อย่างที่ทราบกันดีนะครับใน QGIS หรือ Desktop GIS หลายๆ ตัวนั้นได้เรียกใช้งานไลบรารี่จาก GDAL/OGR ในการประมวลผลหรือวิเคราะห์ข้อมูลใน core หลักของโปรแกรมตนเอง โดย GDAL จะเป็นชุดคำสั่งที่ทำงานกับข้อมูลราสเตอร์ และ OGR จะเป็นชุดคำสั่งที่ทำงานกับข้อมูลเวกเตอร์นั่นเองครับ

มาเริ่มต้นการใช้งาน GDAL/OGR อย่างง่ายกันก่อนนะครับ

GDAL Quick Start

1.เรียกดูข้อมูลราสเตอร์ด้วยชุดคำสั่ง gdalinfo

>> gdalinfo landsat-8.tif

gdalinfo landsat-8.tif
c:\Workspace\gdal_ogr\raster>gdalinfo landsat-8.tif
Driver: GTiff/GeoTIFF
Files: landsat-8.tif
Size is 527, 438
Coordinate System is:
PROJCS["WGS 84 / UTM zone 47N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",99],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32647"]]
Origin = (578355.000000000000000,1891425.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  578355.000, 1891425.000) ( 99d44'11.47"E, 17d 6'21.45"N)
Lower Left  (  578355.000, 1878285.000) ( 99d44' 9.79"E, 16d59'13.88"N)
Upper Right (  594165.000, 1891425.000) ( 99d53' 6.42"E, 17d 6'19.31"N)
Lower Right (  594165.000, 1878285.000) ( 99d53' 4.41"E, 16d59'11.76"N)
Center      (  586260.000, 1884855.000) ( 99d48'38.02"E, 17d 2'46.65"N)
Band 1 Block=527x1 Type=Float32, ColorInterp=Gray
Band 2 Block=527x1 Type=Float32, ColorInterp=Undefined
Band 3 Block=527x1 Type=Float32, ColorInterp=Undefined
Band 4 Block=527x1 Type=Float32, ColorInterp=Undefined
Band 5 Block=527x1 Type=Float32, ColorInterp=Undefined
Band 6 Block=527x1 Type=Float32, ColorInterp=Undefined
Band 7 Block=527x1 Type=Float32, ColorInterp=Undefined

2.การแปลงรูปแบบข้อมูลราสเตอร์ด้วย gdal_translate

>> gdal_translate -of HFA inputfile.tif outputfile.img

c:\Workspace\gdal_ogr\raster>gdal_translate -of HFA landsat-8.tif landsat-8.img Input file size is 527, 438 0...10...20...30...40...50...60...70...80...90...100 - done. c:\Workspace\gdal_ogr\raster>

ในกรณีต้องการแปลงพิกัดด้วยสามารถเพิ่ม -a_srs EPSG:4326

>> gdal_translate -a_srs EPSG:4326 -of HFA inputfile.tif outputfile.img

c:\Workspace\gdal_ogr\raster>gdal_translate -a_srs EPSG:4326 -of HFA landsat-8.tif landsat-8_4326.img Input file size is 527, 438 0...10...20...30...40...50...60...70...80...90...100 - done. c:\Workspace\gdal_ogr\raster>

3.การแปลงระบบพิกัดด้วย gdal_warp

>> gdalwarp -s_srs EPSG:32647 -t_srs EPSG:4326 -of GTiff inputfile.tif outputfile.tif

c:\Workspace\gdal_ogr\raster>gdalwarp -s_srs EPSG:32647 -t_srs EPSG:4326 -r near -of GTiff landsat-8.tif landsat_latlon.tif Creating output file that is 537P x 430L. Processing landsat-8.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done. c:\Workspace\gdal_ogr\raster>

4.การต่อภาพหรือ mosaic ข้อมูลภาพด้วย gdal_warp หรือ gdal_merge.py

>> gdalwarp -of GTiff n15_e099_1arc_v3.tif n15_e100_1arc_v3.tif warpmerged.tif

>> gdal_merge -of GTiff n15_e099_1arc_v3.tif n15_e100_1arc_v3.tif -o merged.tif

** นอกจากนี้ถ้ามีไฟล์จำนวนมากในการ mosaic สามารถให้อ่านจาก Textfile แทนได้เช่นกัน ดังตัวอย่าง

>> gdal_merge -ot UInt16 -of GTiff -optfile inputfiles.txt -o mosaic.tif

5.การสร้าง index ข้อมูลราสเตอร์ด้วย gdalindex

>> gdaltindex -tileindex location tileindex.shp *.tif

C:\Workspace\gdal_ogr\index>gdaltindex -tileindex location -f "GPKG" tileindex.gpkg n15_e099_1arc_v3.tif n15_e100_1arc_v3.tif n16_e099_1arc_v3.tif n16_e100_1arc_v3.tif Creating new index file...

** ถ้าต้องการเปลี่ยนเป็น ESRI Shapefile ดังตัวอย่าง

>> gdaltindex -tileindex location -f “ESRI Shapefile” tileindex.shp *.tif

OGR Quick Start

1.เรียกดูข้อมูลเวกเตอร์ด้วยชุดคำสั่ง ogrinfo

>> ogrinfo -ro geodb.gpkg

c:\Workspace\gdal_ogr\vector>ogrinfo -ro geodb.gpkg INFO: Open of `geodb.gpkg' using driver `GPKG' successful. 1: flood (3D Multi Polygon) 2: landuse (Multi Polygon) 3: urban (Multi Polygon) 4: trans (Multi Line String) 5: transmission (Multi Line String) 6: site (Multi Polygon) 7: river (Multi Polygon)

เรียกดูข้อมูลสรุปแต่ละเลเยอร์ด้วย -so

>> ogrinfo -ro -so geodb.gpkg landuse

c:\Workspace\gdal_ogr\vector>ogrinfo -ro -so geodb.gpkg landuse
INFO: Open of `geodb.gpkg'
using driver `GPKG' successful.
Layer name: landuse
Geometry: Multi Polygon
Feature Count: 820
Extent: (629633.000000, 1826540.000000) - (658385.000000, 1847940.000000)
Layer SRS WKT:
PROJCS["WGS 84 / UTM zone 47N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",99],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32647"]]
FID Column = fid
Geometry Column = geom
AREA: Real (0.0)
PERIMETER: Real (0.0)
LUCODE_52: String (16.0)
DES_TH52: String (192.0)
DES_EN52: String (100.0)
DES_52: String (100.0)
GLU_CODE: String (2.0)

** ถ้าต้องการเปลี่ยนเป็น ESRI Shapefile ดังตัวอย่าง

>> ogrinfo -ro amphoe.shp
>> ogrinfo -ro -so amphoe.shp amphoe

2.การแปลงรูปแบบเวกเตอร์ด้วย ogr2ogr

>> ogr2ogr -formats >> ogr2ogr -f GeoJSON landuse.geojson geodb.gpkg landuse

ogr2ogr -s_srs EPSG:32647 -t_srs EPSG:4326 -f GeoJSON landuse.geojson geodb.gpkg

** ถ้าข้อมูลเป็น ESRI Shapefile

>> ogr2ogr -s_srs EPSG:32647 -t_srs EPSG:4326 -f GeoJSON building.geojson building.shp

Reference:

[1] https://www.gdal.org/gdal_utilities.html

[2] https://www.gdal.org/ogr_utilities.html

[3] https://live.osgeo.org/en/quickstart/gdal_quickstart.html#get-to-know-gdal