## 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')

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

# 先設定投影座標系統是 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)

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

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

# 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)

penghu

taiwan

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

# 請先下載 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'))

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

# 再看一下圖
plot(penghu_mask)

# 輸出澎湖的
raster::writeRaster(penghu_mask, filename = 'Penghu_20m_checked.tif', 
              format = 'GTiff', overwrite = TRUE, options = c('COMPRESS=LZW', 'TFW=YES'))


# 接下來清臺灣本島的
r_max = raster::getValues(taiwan)
which(r_max > 3952)
# check the values
taiwan[which(r_max > 3952)]
# 找到有問題的點位
raster::xyFromCell(taiwan, which(r_max > 3952))

# 這個效率不怎麼好
# taiwan[taiwan>3952] = NA 
# 用笨方法，其實可以用  lapply，但不知道為什麼沒法成功使用在 raster layer 中
taiwan[96676141] = NA; taiwan[96910164] = NA; taiwan[96920339] = NA

# approximately 50 mins
taiwan_mask = raster::mask(taiwan, taiwan_boundary_epsg3826)

plot(taiwan_mask)

# 輸出初步檢查過的臺灣本島 DEM
raster::writeRaster(taiwan_mask, filename = 'Taiwan_20m_checked.tif', 
              format = 'GTiff', overwrite = TRUE, options = c('COMPRESS=LZW', 'TFW=YES'))

taiwan_masked_val = raster::getValues(taiwan_mask)

hist(taiwan_masked_val)

taiwan_mask

# 使用 focal 檢查，weight 矩陣為 3x3。先計算標準差
# Penghu
penghu_focal = raster::focal(penghu_mask, w=matrix(1,3,3), fun = sd,  pad = T, padValue = 0)

# Taiwan  (計算需要大約 1.5 hr)
taiwan_focal = raster::focal(taiwan_mask, w=matrix(1,3,3), fun = sd,  pad = T, padValue = 0)

# 先輸出
raster::writeRaster(penghu_focal, filename = 'Penghu_20m_focal_sd.tif', 
              format = 'GTiff', overwrite = TRUE, options = c('COMPRESS=LZW', 'TFW=YES'))
raster::writeRaster(taiwan_focal, filename = 'Taiwan_20m_focal_sd.tif', 
              format = 'GTiff', overwrite = TRUE, options = c('COMPRESS=LZW', 'TFW=YES'))

taiwan_focal

penghu_focal

taiwan_focal_val = getValues(taiwan_focal)