Notebook

檢查內政部 20m DEM 資料¶

最近內政部釋出了全國的 20m DEM 數值資料，但志展從原始資料產製出等高線圖時發現了一些數值上不合理的錯誤(請見參考資料[1])，所以這個文件用 R 來示範如何檢查這類的資料錯誤。

information compiled by Cheng-Tao Lin (mutolisp@gmail.com)

Part 1. 前言及概念¶

1.1 資料集來源¶

1.2 資料的結構¶

內政部釋出的資料是 raster 格式(Geotiff)，簡單來說是使用格柵式(grid)的排列方式把數值填入，簡單說資料結構的概念如下，在表頭的部份會記錄詮釋資料，例如座標系統、維度、格式等資訊。底下則會用像矩陣一樣排列數值：

   # 詮釋資料(座標系統、維度、格式)
   2 3 5 6 1
   5 4 1 3 5
   2 2 2 5 4
   1 2 1 3 1
   1 1 1 1 1

1.3 座標系統¶

跟空間有關的資料集中，座標系統是最重要的詮釋資料(metadata)，所以先看一下來源的 URL 有沒有標註？檢查之後沒有發現任何詮釋資料有關座標系統的描述，因此我們先下載資料後先使用 libGDAL 中的 gdalinfo 確認一下座標系統(後面會再用 R 的 code 詳細解釋如何使用 gdal 查詢)

臺灣地區 DEM (filename: 20m_dem.tif)

Driver: GTiff/GeoTIFF
Files: dem_20m.tif
Size is 10175, 19112
Coordinate System is:
PROJCS["TWD_1997_TM_Taiwan",
    GEOGCS["GCS_TWD_1997",
        DATUM["TWD_1997",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",121],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",250000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (148310.000000000000000,2801730.000000000000000)
Pixel Size = (20.000000000000000,-20.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  148310.000, 2801730.000) (119d59'23.86"E, 25d19'16.05"N)
Lower Left  (  148310.000, 2419490.000) (120d 0'57.80"E, 21d52'12.03"N)
Upper Right (  351810.000, 2801730.000) (122d 0'40.43"E, 25d19'16.02"N)
Lower Right (  351810.000, 2419490.000) (121d59' 6.38"E, 21d52'12.01"N)
Center      (  250060.000, 2610610.000) (121d 0' 2.12"E, 23d35'56.25"N)
Band 1 Block=128x128 Type=Float64, ColorInterp=Gray
  NoData Value=-1.79769313486231571e+308
  Metadata:
    RepresentationType=ATHEMATIC

澎湖地區 DEM (filename: Penghu_20m_dem.tif)

Driver: GTiff/GeoTIFF
Files: Penghu_20m.tif
Size is 2220, 3495
Coordinate System is:
PROJCS["TWD_1997_TM_Taiwan",
    GEOGCS["GCS_TWD_1997",
        DATUM["TWD_1997",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",121],
    PARAMETER["scale_factor",0.9999],
    PARAMETER["false_easting",250000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (76661.476474763854640,2633852.793932494707406)
Pixel Size = (20.000000000000000,-20.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   76661.476, 2633852.794) (119d17'56.69"E, 23d47'58.06"N)
Lower Left  (   76661.476, 2563952.794) (119d18'25.77"E, 23d10' 6.71"N)
Upper Right (  121061.476, 2633852.794) (119d44' 4.81"E, 23d48'13.14"N)
Lower Right (  121061.476, 2563952.794) (119d44'26.44"E, 23d10'21.34"N)
Center      (   98861.476, 2598902.794) (119d31'13.49"E, 23d29'10.38"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
  NoData Value=-3.40282306073709653e+38

所以兩個檔案應該都是 TWD1997 Zone 121 (EPSG:3826)，但澎湖的中央子午線(central meridian)也是 121º，有點奇怪，理論上澎湖的投影座標系統都是以 119ºE 當做中央子午線，後面我們再檢查確認。

參考資料¶

[1] [應用內政部20公尺網格數值地形模型資料](http://blog.nutsfactory.net/2016/09/14/taiwan-moi-20m-dtm/)
[2] [TWD1997 TM2 Zone 121, EPSG:3826](http://spatialreference.org/ref/epsg/3826/)
[3] [TWD1997 TM2 Zone 119, EPSG:3825](http://spatialreference.org/ref/epsg/3825/)

Part 2. 實做¶

預先安裝的程式¶

libgdal

如果是 Windows 的使用者可以安裝 OSGEO4W；若你是 Mac 使用者可以參考使用 homebrew 安裝 gdal-20；GNU/Linux 使用者，我想你應該知道怎麼安裝 XD

R，需要的 package 如下：
- raster
- rgdal
- sp
RStudio(選擇性安裝)

Step 2.1 設定工作目錄及看詮釋資料¶

第一步驟我們將載入相關的 packages (raster 及 rgdal)，設定工作目錄以及看資料集本身的詮釋資料，這個步驟很重要，有助於了解座標系統、空間中的範圍、網格的維度(列數、欄數)、沒有資料的空值等

In [1]:

## Preprocess of 20 m DEM published by Minister of Interior (MOI)
# original dataset url: http://data.gov.tw/node/35430
# 載入 raster 和 rgdal 兩個 package
require(raster)
require(rgdal)
require(sp)

setwd('/Users/psilotum/Documents/Dropbox/Databases/GISLayers/MOI_20m_DEM/')
# get information from gdalinfo
# 先看一下臺灣本島的 raster info
rgdal::GDALinfo('./raw/dem_20m.tif')

Loading required package: raster
Loading required package: sp
Loading required package: rgdal
rgdal: version: 1.1-10, (SVN revision 622)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.0, released 2016/04/25
 Path to GDAL shared files: /usr/local/Cellar/gdal-20/2.1.0/share/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-3 
Warning message:
: statistics not supported by this driver

rows        19112 
columns     10175 
bands       1 
lower left origin.x        148310 
lower left origin.y        2419490 
res.x       20 
res.y       20 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0
+ellps=GRS80 +units=m +no_defs 
file        ./raw/dem_20m.tif 
apparent band summary:
   GDType hasNoDataValue    NoDataValue blockSize1 blockSize2
1 Float64           TRUE -1.797693e+308        128        128
apparent band statistics:
         Bmin       Bmax Bmean Bsd
1 -4294967295 4294967295    NA  NA
Metadata:
AREA_OR_POINT=Area

In [2]:

# 再看澎湖的
rgdal::GDALinfo('./raw/Penghu_20m.tif')

Warning message:
: statistics not supported by this driver

rows        3495 
columns     2220 
bands       1 
lower left origin.x        76661.48 
lower left origin.y        2563953 
res.x       20 
res.y       20 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0
+ellps=GRS80 +units=m +no_defs 
file        ./raw/Penghu_20m.tif 
apparent band summary:
   GDType hasNoDataValue   NoDataValue blockSize1 blockSize2
1 Float32           TRUE -3.402823e+38        128        128
apparent band statistics:
         Bmin       Bmax Bmean Bsd
1 -4294967295 4294967295    NA  NA
Metadata:
AREA_OR_POINT=Area

看起來兩個 raster 的資料類型和 NoDataValue (沒有資料的數值)都不統一，臺灣的資料類型是 Float64、澎湖的則是 Float32 (話說 20m DEM 有必要使用到 Float64 嗎？ XD)。NoDataValue 則是分別為 -1.797693e+308 及 -3.402823e+38。

Step 2.2 載入資料¶

使用 rgdal 載入 GeoTiff 檔案，並用 raster 轉換成 RasterLayer 格式

In [3]:

# 先設定投影座標系統是 epsg:3826
taiwan_moi_20m = rgdal::readGDAL('./raw/dem_20m.tif', p4s = '+init=epsg:3826')
penghu_moi_20m = rgdal::readGDAL('./raw/Penghu_20m.tif', p4s = '+init=epsg:3826')

# convert SpatialGridDataFrame to RasterLayer object (raster library)
penghu = raster::raster(penghu_moi_20m)
taiwan = raster::raster(taiwan_moi_20m)

./raw/dem_20m.tif has GDAL driver GTiff 
and has 19112 rows and 10175 columns
./raw/Penghu_20m.tif has GDAL driver GTiff 
and has 3495 rows and 2220 columns

Step 2.3 檢查資料的錯誤¶

In [4]:

# 再看一次載入的 raster 資訊的摘要
# check the layer information
penghu
taiwan

class       : RasterLayer 
dimensions  : 3495, 2220, 7758900  (nrow, ncol, ncell)
resolution  : 20, 20  (x, y)
extent      : 76661.48, 121061.5, 2563953, 2633853  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3826 +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : in memory
names       : band1 
values      : -999, 70.5  (min, max)

class       : RasterLayer 
dimensions  : 19112, 10175, 194464600  (nrow, ncol, ncell)
resolution  : 20, 20  (x, y)
extent      : 148310, 351810, 2419490, 2801730  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3826 +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : in memory
names       : band1 
values      : -999, 6246.99  (min, max)

疑？values (代表 raster 數值) 顯示載入 raster 檔案的最大值和最小值，澎湖地區的圖層最大值 70.5 (單位應該是公尺)，最小值 -999；臺灣地區的最大值是 6246.99 ，最小值也是 -999。這和我們認知中臺灣的海拔高度範圍有極大的落差。最小值應該不可能為 -999，這已經接近海溝的深度了，合理懷疑可能是處理資料時原始軟體把 -999 當成「沒有資料的數值(No data value)」，一般來說如果使用 ErDAS 或 ArcGIS 的某些 raster 格式會把沒有資料的數值設定為 -9999 或 -999。

所以我們要先找出不合理的數值並試圖修正，臺灣的海拔範圍應該是 0–3952m (最高峰玉山山頂)，但考量到平地有湖泊池沼、海岸也有可能會有低於海平面(0 m)的點，因此先將不合理數值切成兩部分：

大於最大值的數值 ( x > 3952 )
小於海拔 0 公尺的數值

對於大於最大值的數值，實作則是用 hist() 來繪製資料數值的頻度分布，有助於快速觀察並發現問題。若小於海拔 0 公尺的，我們則使用遮罩(mask)的概念，先將臺灣國界外的落海區域先去除，再來詳細檢查國界內是否有小於海拔 0 公尺太多的數值。

In [5]:

# 先用 raster::getValues() 取得圖層的數值
penghu_val = raster::getValues(penghu)
taiwan_val = raster::getValues(taiwan)
# 再用 hist 繪圖
hist(penghu_val)
hist(taiwan_val)

看起來應該是 no data value，所以我們使用 raster::calc() 把 -999 處理成 NA

In [6]:

# calc 是用 function 來計算 raster 中的數值，所以先
# 建立一個 omit_na function 來處理 -999 的數值
omit_na = function(x){ 
  x[x == -999] = NA
  return(x)
}

penghu = raster::calc(penghu, fun = omit_na)
taiwan = raster::calc(taiwan, fun = omit_na)

In [7]:

penghu

class       : RasterLayer 
dimensions  : 3495, 2220, 7758900  (nrow, ncol, ncell)
resolution  : 20, 20  (x, y)
extent      : 76661.48, 121061.5, 2563953, 2633853  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3826 +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : in memory
names       : layer 
values      : -0.5, 70.5  (min, max)

In [8]:

taiwan

class       : RasterLayer 
dimensions  : 19112, 10175, 194464600  (nrow, ncol, ncell)
resolution  : 20, 20  (x, y)
extent      : 148310, 351810, 2419490, 2801730  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3826 +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : /private/var/folders/ys/z84vdxw10dg46yk5zzj991qh0000gn/T/RtmpwrzejI/raster/r_tmp_2016-09-16_233716_18031_76851.grd 
names       : layer 
values      : -33.29, 6246.99  (min, max)

In [9]:

# 畫個圖確認一下
plot(penghu)
plot(taiwan)

看起來外海還有許多有資料的數值，所以下個步驟就是載入內政部的直轄市、縣界圖當成遮罩，把這些外海的資料去除

Step 2.4 使用遮罩除去國界外的數值¶

使用 raster::mask(raster, shape_polygon) 來處理

In [10]:

# 請先下載 http://data.gov.tw/node/7442 並解壓縮至 County1041215A 的目錄中
# 政府開放資料中的 shapefile 投影座標系統是 EPSG:3824，為了和 raster 的投影座標系統相符
# 所以我們先載入後再轉換成 EPSG:3826 
taiwan_boundary = rgdal::readOGR('./County1041215A/County_MOI_1041215.shp', 
                                 layer = 'County_MOI_1041215', p4s='+init=epsg:3824')
taiwan_boundary_epsg3826 = sp::spTransform(taiwan_boundary, CRS('+init=epsg:3826'))

OGR data source with driver: ESRI Shapefile 
Source: "./County1041215A/County_MOI_1041215.shp", layer: "County_MOI_1041215"
with 22 features
It has 4 fields

Warning message:
In rgdal::readOGR("./County1041215A/County_MOI_1041215.shp", layer = "County_MOI_1041215", : p4s= argument given as: +init=epsg:3824
 and read as: +proj=longlat +ellps=GRS80 +no_defs 
 read string overridden by given p4s= argument value

In [11]:

# 澎湖的比較單純，所以先試做澎湖的
penghu_mask = raster::mask(penghu, taiwan_boundary_epsg3826)

In [12]:

# 再看一下圖
plot(penghu_mask)

In [13]:

# 輸出澎湖的
raster::writeRaster(penghu_mask, filename = 'Penghu_20m_checked.tif', 
              format = 