MODIS のデータを R で扱ってみる（地表面温度 LST）

MODIS（MODerate resolution Imaging Spectroradiometer）のデータを R で扱う機会があったのでメモ。
なお、{MODIS} という R のパッケージがあり、取得から加工までできてしまうみたいだが今回これは使わず、データの取得はNASA のサイト EARTH DATA から DL、hdf への変換は cygwin ということにしておきます。

MODIS の 地表面温度（LST: Land Surface Temperature）に関するプロダクト MOD11A1 (LST_Day_1km) を可視化を含む分析において扱いやすくするために、

1. hdf 形式から GeoTiff 形式に変換し、
2. 座標を WGS1984 に変換し、
3. 物理量（ケルビンとか温度（セルシウス度））に変換する

コードをメモします。ポイント 1, 2 については MODIS Tool とか HEGとかがあるが、それを R で代替しようとしたものです。全部 library(MODIS) でやればいいのに。

library(gdalUtils)
library(raster)
library(rgdal)
####

## 1. GeoTiff への変換
setwd("C:/hoge/piyo/MOD11A1_hdf") # cygwin で変換した .hdf を保存したフォルダ 
for(i in 1:length(list.files())){
  gdalUtils::gdalinfo(list.files()[i])
  dat <- gdalUtils::get_subdatasets(list.files()[i])
  gdalUtils::gdal_translate(dat[1],  # 1 が LST_Day_1km
                 dst_dataset=paste("C:/hoge/piyo/MOD11A1_tif/LST_Day_", i, ".tif", sep=""), of="GTiff")
  print(i)
}

## 2. 座標を WGS1984 に変換
setwd("C:/hoge/piyo/MOD11A1_tif")
for(i in 1:length(list.files())) {
    gdalUtils::gdalwarp(paste("LST_Day_", i, ".tif", sep=""),
                        paste("C:/hoge/piyo/MOD11A1_tif_prj/LST_Day_", i, ".tif", sep=""),
                        t_srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs", # to WGS1984 
                        overwrite=T)
  print(i)
}

## 3. 物理量に変換する（ついでに関心の対象範囲に clip/crop する）
setwd("C:/hoge/piyo/MOD11A1_tif_prj")
TargetArea <- readOGR("C:/hoge/piyo/TargetArea.shp") # たとえば、国土数値情報行政区域から関東地方のshp
#plot(TargetArea)

for(i in 1:length(list.files())){
  dat2 <- raster(list.files()[i])
  values(dat2) <- values(dat2)*0.02 - 273.15
  # 0.02 は MOD11A1 の scale値；　273.15 はケルビンからセルシウス度への変換式
  dat3 <- crop(dat2, TargetArea)
  writeRaster(dat3, filename=paste("C:/hoge/piyo/MOD11A1_fin/LST_Day_", i, ".tif", sep=""), 
              format="GTiff", overwrite=TRUE)
  print(i)
}

scale 値というのは、観測値（the digital values）を物理量（ここではケルビン）に変換するための値です（参考： https://lpdaac.usgs.gov/products/mod11a1v006/ ）。ここではさらに温度とするために -273.15 しています。

MOD11A1 の LST_Day_1km は、広域に毎日とれるが、雲があったりで見えるところ／見ないところが毎日違ったりする。

関東地方の2020年のある日の地表面温度分布（ブランクは NA）

（もう一度）全部 library(MODIS) でやればいいのに。