Rを用いた衛星データ(LST, NDVI)の空間統計分析

ここでは、衛星リモートセンシングデータから生成された地表面温度 (Land Surface Temperature: LST) と植生指数 (Normalized Difference Vegetation Index: NDVI) を使って、空間統計学の解析を通じて、LSTを推定する際の標準誤差が大きい地点を可視化する演習を説明しています。統計処理用にRを、またデータの前処理にシェルスクリプト、Python、可視化にQGISを使っています。

LST, NDVIの定義式については、「反射率、輝度温度、標高データを用いた土地被覆分類」に記載されています。

また本演習では、NASAが打ち上げて運用しているModerate Resolution Imaging Spectrometer (MODIS) というセンサで取得された画像を利用しています。空間分解能250 m〜1 kmで、同一地点を毎日取得できる低空間分解能の衛星センサです。画像の空間解像度は落ちますが、その代わりに毎日画像を取得できることで、広域の長期間にわたる時系列変化を調べるのに適しています。36バンドで得られたデータから、陸域だけでなく、海洋、大気、雪氷圏の解析に役立つ様々なプロダクトが生成されています。

    本演習内容:
  1. (参考)Cygwinのインストール
  2. R, QGISのインストール
  3. (参考)使用するデータ
  4. データの前処理
    1. 必要な範囲のLST, NDVIデータのコピー (シェルスクリプト)
    2. LST, NDVIデータの統合 (Python)
  5. 重回帰モデルのパラメータ推定
  6. QGIS上での結果の表示
  7. 課題
時空間モデリングによる地表面温度(LST)推定の標準誤差 この演習で最終的に得られる結果の一例:亀岡市のDOY=225(2012年1月1日をDOY=1)時における地表面温度(LST)推定の標準誤差

注意:以下で(参考)と書いてあるところは、演習中には自分で処理する必要はなく、既に用意されたデータ等を使用すればよい。今後、自分でさらに処理する時のために記述した。

1. (参考)Cygwinのインストール

MODISというセンサで取得されたデータから生成されるLST, NDVIのプロダクトは、各々8日間、16日間ごとに処理される。通常、WebサイトからLST, NDVIのプロダクトをダウンロードする際には対象期間に含まれるプロダクトを全てダウンロードする。

本演習では、LST, NDVIのプロダクトを一つの表にまとめる。LSTは8日間ごとにファイルが存在するが、処理で必要とされるのは16日ごとのファイルであり、一旦ダウンロードしたファイルのうちコピーするファイルを選別する必要が出てくる。手作業でコピーすることもできるが非効率的で、シェルスクリプトを使うと簡単に実行できる。そのシェルスクリプトを実行する環境を構築するために、cygwinをインストールする。

2. R, QGISのインストール

以下のサイトからインストーラーをダウンロードし、インストールする。

3.(参考)使用するデータ

演習では、演習内に提示するデータを利用して処理するが、自分たちでも直接LST, NDVIのプロダクトをダウンロードできる。

(参考)2022年3月修士課程修了の楠瀬智也氏が残してくれた、LST, NDVIプロダクトのダウンロードおよび処理の手引書

4. データの前処理

空間統計解析の前に、上記のデータを処理する必要がある。

4.1.(参考)必要な範囲のLST, NDVIデータのコピー (シェルスクリプト)

衛星画像(データ)の処理においてはファイル数が多数に上り、またファイル名に規則性があることから、自動的に一括処理するシェルスクリプトが多用される。元々、筆者は1990年代後半から2000年代半ばまでUNIX (Solaris, HP-UX) で処理することが多く、シェルスクリプトを多用していた。現在はWindowsでの作業が多く、cygwinを使うことがあるため、cygwinを前提に下記の通り説明している。一括して大量のファイルをコピーしたり、ファイル名を変更したりできれば、別にシェルスクリプトでなくても構わない。

(参考)シェルスクリプト講座: [基礎編] [制御文編] [実例編1] [実例編2]

まず、この作業で想定している、フォルダ(またはディレクトリ)とファイルの構造を下記に示す。ドライブDの下のフォルダ名を「rsgis」としているが、この名前は何でも構わない。

C:\Program Filesなど

D:\rsgis
   |__ Data 
        |__ Prefecture_LST (ダウンロードしてきたLST:8日ごと)
             |__ Aichi
                  |__ Agui-cho_2012001_LST.csv
                  |__ Agui-cho_2012009_LST.csv
                  |__ Agui-cho_2012017_LST.csv
             |__ Akita
                  |
             (他の都道府県のフォルダが存在)

        |__ Prefecture_NDVI(ダウンロードしてきたNDVI:16日ごと)
             |__ Aichi
                  |__ Agui-cho_2012001_NDVI.csv
                  |__ Agui-cho_2012017_NDVI.csv
                  |__ Agui-cho_2012033_NDVI.csv
             |__ Akita
                  |
             (他の都道府県のフォルダが存在)

        |__ Prefecture(上記2つからコピーしたファイルの保存用)

        |__ Prefecture_rst(4年間のLSTとNDVIを1つにまとめたファイルを保存)

Cygwinで実行する内容:打ち込むコマンドを赤字で示す。ここでは、全ての都道府県のデータをコピーする処理を示しているが、演習では特定の都道府県名を明記して処理する。。

$ pwd ← 現在のディレクトリを表示
/home/susaki
$ cd /cygdrive/d/rsgis ← データのあるディレクトリへ移動

$ ls ← ファイルやディレクトリを表示
(結果が表示される:省略)

$ ./copy.sh

copy.shの中身

#!/usr/bin/bash

# LSTデータ:NDVIに合わせて16日おきのデータをコピー
# NDVIデータ:丸ごとコピー

DIR_OUT=Data/Prefecture

DIR_IN1=Data/Prefecture_LST # 001, 009, 017と8日おきのデータ
DIR_IN2=Data/Prefecture_NDVI # 001, 017, 033と16日おきのデータ

# 対象年:開始、終了、間隔
year_st=2012
year_end=2015
year_intv=1

# 対象DOY (Day of year):開始、終了、間隔
doy_st=1
doy_end=353
doy_intv=16

# LSTデータのコピー:NDVIに合わせて16日おきのデータをコピー
for dir in `ls ${DIR_IN1}`
do
	echo "++++++ ${dir} ++++++"
	dirpath_in=${DIR_IN1}/${dir}
	dirpath_out=${DIR_OUT}/${dir}
	echo ${dirpath_in}
	echo ${dirpath_out}
	
	if [ ! -e ${dirpath_out} ] # 出力先ディレクトリが存在しない場合
	then
		mkdir ${dirpath_out} # 出力先ディレクトリを作成
	fi
	
	year=${year_st}
	
	while [ ${year} -le ${year_end} ]
	do
		echo "***** ${year} *****"
		doy=${doy_st}
		
		while [ ${doy} -le ${doy_end} ]
		do
			doyid=`printf "%03d" ${doy}` # DOYの3桁表示
			yeardoy=${year}${doyid}
			echo "----- ${yeardoy} -----"
			cp ${dirpath_in}/*${yeardoy}*csv ${dirpath_out}
			doy=`expr ${doy} + ${doy_intv}`
		done
		year=`expr ${year} + ${year_intv}`
	done
done


# NDVIデータのコピー:全ファイルをそのままコピー
for dir in `ls ${DIR_IN2}`
do
	echo "++++++ ${dir} ++++++"
	dirpath_in=${DIR_IN2}/${dir}
	dirpath_out=${DIR_OUT}/${dir}
	echo ${dirpath_in}
	echo ${dirpath_out}
	
	if [ ! -e ${dirpath_out} ] # 出力先ディレクトリが存在しない場合
	then
		mkdir ${dirpath_out} # 出力先ディレクトリを作成
	fi

	cp ${dirpath_in}/*csv ${dirpath_out}
done

4.2. LST, NDVIデータの統合 (Python)

以下、簡単に処理の方針を示す。

  1. 市町村単位で関連するデータ(LST,NDVI)を収集する。
  2. ファイルごとに含まれる点数が異なるため、まず全ファイルを探索して、Lat, Lonの和集合を作成する。
  3. Lat, Lonのテーブルに対し、16日ずつDOYをずらしながらファイルが存在するか確認する。ファイルが存在する場合には、存在する地点に対しデータを格納し、存在しない地点にはN/Aを付与する。ファイルが存在しない場合(特にLST)、全点に対しN/Aを付与する。
#!/usr/bin/env python
# coding: utf-8

import pandas as pd
import numpy as np
import os
import glob
import csv

# DOYのリスト作成
year_list = list(range(2012, 2016))
doy_list = list(range(1, 365, 16))
doy_name = []
for year in year_list:
    for doy in doy_list:
        doy_tmp = f'{doy:03}' # 数字を3桁の文字として出力
        doy_name_tmp = str(year) + doy_tmp # 例:2012 + 001 -> 2012001
        doy_name.append(doy_name_tmp)

# ファイルが存在するフォルダのPATH
path0 = "./Data/Prefecture/Kyoto"
path0_out = "./Data/Prefecture_rst/Kyoto"

# プロダクトリスト(2種類)
prods = ['LST', 'NDVI']


files0 = os.listdir(path0)
f0 = files0[0:]
num0 = len(f0)

# 市町村名の抽出
areas = []
#for infile in f0:
for i in range(num0):
    # '2'を使って文字列を分割
    # 京都市のファイルが不規則に命名されているため、'_'では失敗する
    strtmp = f0[i].split('2') 
    areas.append(strtmp[0])

# 市町村名の重複削除
areas_unique = list(dict.fromkeys(areas))
num_areas = len(areas_unique)

# 市町村単位で2012001_LST(NDVI), ..., 2015361_LST(NDVI)のファイルを処理
for area_name in areas_unique:
    
    # Lat, Lonの和集合作成:開始
    pathtmp = path0 + "/" + area_name + "*.csv"
    print(pathtmp)
    area_file = glob.glob(pathtmp) # 該当する市町村のファイルのみ抽出
    print(area_file)
    
    num_area_file = len(area_file)
    
    ll = []
    for i in range(1):
        
        # 各CSVファイルのデータを読み込み
        with open(area_file[i], encoding='utf8', newline='') as f:
            csvdata = csv.reader(f)
            lines = [row for row in csvdata]
            num_lines = len(lines)
            for j in range(1, num_lines):
                #print(j, num_lines)
                ll_tmp = [lines[j][1], lines[j][2]]
                ll.append(ll_tmp)
        print(i, 'old ll = ', len(ll))
    
    for i in range(1, num_area_file):
        # 各CSVファイルのデータを読み込み
        with open(area_file[i], encoding='utf8', newline='') as f:
            csvdata = csv.reader(f)
            lines = [row for row in csvdata]
            num_lines = len(lines)
            for j in range(1, num_lines):
                ll_tmp = [lines[j][1], lines[j][2]]
                ll.append(ll_tmp)
        
        #重複削除
        print(i, 'old ll = ', len(ll))
        ll2 = []
        for line in ll:
            if line not in ll2:
                ll2.append(line)
        ll = ll2
        print(i, 'new ll = ', len(ll))
    
    
    # 該当するデータを格納。存在しない場合にはN/Aを格納
    header = ['Lon', 'Lat']
    data = ll
    num_lines_all = len(data)
    
    # Lat, Lonの和集合作成:終了
    
    
    # 2012001, ..., のように順番にファイルが存在するか点検
    # ファイルが存在しない場合、その日のデータには全てN/Aを格納
    
    flag_data = [0] * num_lines_all # データの有無のフラッグ
    
    # DOYのリスト順にファイルが存在すれば開く
    for yeardoy in doy_name:
        
        for prod in prods:
            print(yeardoy, prod)
            yrdyprod = yeardoy + prod
            flag_file = 0
            for afile in area_file:
                if afile.find(prod) > 0 and afile.find(yeardoy) > 0:
                    #print(afile)
                    flag_file = 1
                    break
            if flag_file == 1:
                
                # 各CSVファイルのデータを読み込み
                with open(afile, encoding='utf8', newline='') as f:
                    csvdata = csv.reader(f)
                    lines = [row for row in csvdata]
                    num_lines = len(lines)
                    
                    # "2012001LST"などをヘッダーに追加
                    header.append(yrdyprod)
                    
                    # フラッグの初期化
                    for k in range (num_lines_all):
                        flag_data[k] = 0
                    
                    # Lat, Lonが一致すればLST(NDVI)データの書き込み
                    for j in range(1, num_lines):
                        for k in range(num_lines_all):
                            if lines[j][1] == data[k][0] and 
                               lines[j][2] == data[k][1]:
                                data[k].append(lines[j][3]) # LST or NDVI
                                flag_data[k] = 1
                                break
                    
                    for k in range(num_lines_all):
                        if flag_data[k] == 0:
                            data[k].append('N/A')
                    
            else: # ファイルが存在しない場合
                # "2012001LST", "2012001NDVI"などをヘッダーに追加
                header.append(yrdyprod)
                for k in range(num_lines_all):
                    data[k].append('N/A')
            
    # ファイルへの出力
    file_path = path0_out + "/" + area_name + "LST_NDVI.csv"
    
    data_tmp = data
    data_tmp.insert(0, header)
    
    with open(file_path, "w", encoding="utf-8") as fout:
        writer = csv.writer(fout, lineterminator='\n')
        writer.writerows(data_tmp)

5. 重回帰モデルのパラメータ推定 (R)

以下の手順に従って入力していく。「QGIS, Rを用いた公示地価データの空間統計分析:静的な時空間モデリング」と同様に、低ランクガウス過程による時空間統計解析を行う。

下記の例では、i=15 (DOY=225)の結果を出力している。演習では、自動的に連続してi=1から92(2015年のDOY=353:2012/01/01を起点とすると1449)まで出力させたい。(Rではfor文を使うと連続出力ができなかった。Pythonでできるかもしれないので、各自で試してほしい)

#ライブラリの読み込み
install.packages("sp")
install.packages("mgcv")
install.packages("spacetime")
library(sp)
library(mgcv)
library(spacetime)

#データの読み込み
dir <- "Data/Prefecture_rst/Kyoto"
city_name <- "Kameoka-shi"

infile <- paste0(dir, "/", city_name, "_LST_NDVI.csv", sep = "")
dat <- read.csv(infile)

# 奇数:モデル作成用, 偶数:検証用
dat_tmp <- split(dat, rep(c("odd", "even")))
dat1 <- dat_tmp$odd # モデル作成用
dat2 <- dat_tmp$even # 検証用

# 各年度のデータを下記の並びに整理
# Lon, Lat, LST, NDVI, Time

dat1a <- data.frame(Lon=dat1$Lon, Lat=dat1$Lat, 
LST=dat1$X2012001LST, NDVI=dat1$X2012001NDVI, Time=1)

dat2a <- data.frame(Lon=dat2$Lon, Lat=dat2$Lat, 
LST=dat2$X2012001LST, NDVI=dat2$X2012001NDVI, Time=1)

year_list <- list(2012:2015)
doy_list <- seq(from = 1, to = 353, by = 16)

doy2012_list <- paste0("2012", sprintf("%03d", doy_list))
doy2013_list <- paste0("2013", sprintf("%03d", doy_list))
doy2014_list <- paste0("2014", sprintf("%03d", doy_list))
doy2015_list <- paste0("2015", sprintf("%03d", doy_list))

total_doy2012_list <- doy_list
total_doy2013_list <- 366 + doy_list
total_doy2014_list <- 366 + 365 + doy_list
total_doy2015_list <- 366 + 365*2 + doy_list

doy_all <- c(doy2012_list, doy2013_list, doy2014_list, doy2015_list)
tdoy_all <- c(total_doy2012_list, total_doy2013_list, 
total_doy2014_list, total_doy2015_list)



# 2012017, 2012033などのリストを生成
num_doy <- length(doy_all)

for (i in 2:num_doy) { # "2012001"は除外するため、i=2から開始

	# "X2012017LST", "X2012017NDVI"のような列名を作成
	lstid <- paste0("X", doy_all[i], "LST")
	ndviid <- paste0("X", doy_all[i], "NDVI")
	timeid <- tdoy_all[i]

	dat1b <- data.frame(Lon=dat1["Lon"], Lat=dat1["Lat"], 
      LST=dat1[lstid], NDVI=dat1[ndviid], Time=timeid)

	dat2b <- data.frame(Lon=dat2["Lon"], Lat=dat2["Lat"], 
      LST=dat2[lstid], NDVI=dat2[ndviid], Time=timeid)


	# 列名を"X2012017LST"から"LST"へ置換
	dat1b <- dplyr::rename(dat1b, LST=lstid)
	dat2b <- dplyr::rename(dat2b, LST=lstid)

	dat1b <- dplyr::rename(dat1b, NDVI=ndviid)
	dat2b <- dplyr::rename(dat2b, NDVI=ndviid)

	dat1a <- rbind(dat1a, dat1b)
	dat2a <- rbind(dat2a, dat2b)
}

dat1 <- dat1a
dat2 <- dat2a


# N/Aが含まれる行を削除
dat1tmp <- subset(dat1, !(is.na(dat1$LST)))
dat1 <- dat1tmp
dat1tmp <- subset(dat1, !(is.na(dat1$NDVI)))
dat1 <- dat1tmp

dat2tmp <- subset(dat2, !(is.na(dat2$LST)))
dat2 <- dat2tmp
dat2tmp <- subset(dat2, !(is.na(dat2$NDVI)))
dat2 <- dat2tmp

# LSTが文字列として認識されたいたので、数値へ変換
dat1$LST <- as.numeric(dat1$LST)
dat2$LST <- as.numeric(dat2$LST)

# 低ランクガウス過程による時空間統計解析の実行
f1 <- LST ~ te(Lon, Lat , Time, k=c(50,10),d=c(2,1)) + 
s(NDVI) + s(Lon) + s(Lat) + s(Time)
model1 <- gam(f1, data= dat1)

# 結果の出力
est_rst <- paste0(dir, "/", city_name, "_LST_NDVI_rst.txt", sep = "")
sink(est_rst)
summary(model1)
sink()


# 検証
pred_rst <- predict(model1, newdata = dat2, se.fit = TRUE, interval =
'confidence', level = 0.95)

pred_rst_table   <- data.frame(dat2[, c("Lon", "Lat", "Time", "LST")], 
pred=pred_rst$fit, se.pred=pred_rst$se.fit)

# 検証結果の出力
old.op <- options(max.print=2147483647) # 出力桁数の上限の変更
outfile <- paste0(dir, "/", city_name, "_LST_NDVI_val_rst.txt", sep = "")
sink(outfile)
pred_rst_table
sink()


# 全データに対して予測
dat_all <- rbind(dat1, dat2)
pred_rst <- predict(model1, newdata = dat_all, se.fit = TRUE, interval =
'confidence', level = 0.95)

est_table   <- data.frame(dat_all[, c("Lon", "Lat", "Time", "LST")], 
Pred=pred_rst$fit, se.pred=pred_rst$se.fit)
est_table$LST <- as.numeric(est_table$LST)

dif <- est_table$LST - est_table$Pred
est_table <- cbind(est_table, dif)

# i=15 (DOY=225)の時の結果を出力:手動で修正して出力する。
# for文で繰り返して出力させようとしたができなかったため。

i=15
est_table_tmp   <- est_table[est_table$Time == tdoy_all[i], ]
outfile <- paste0(dir, "/", city_name, "_LST_NDVI_est_", 
tdoy_all[i], ".txt", sep = "")
sink(outfile)
est_table_tmp
sink()

低ランクガウス過程による時空間統計解析の実行結果(亀岡市の例)

Family: gaussian 
Link function: identity 

Formula:
LST ~ te(Lon, Lat, Time, k = c(50, 10), d = c(2, 1)) + s(NDVI) + 
    s(Lon) + s(Lat) + s(Time)

Parametric coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 290.93436    0.01799   16174   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                     edf  Ref.df        F  p-value    
te(Lon,Lat,Time) 269.045 349.267   90.077  < 2e-16 ***
s(NDVI)            7.007   8.099   59.630  < 2e-16 ***
s(Lon)             8.103   8.358    4.168 2.44e-05 ***
s(Lat)             7.378   7.991    6.236  < 2e-16 ***
s(Time)            9.000   9.000 3775.343  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.952   Deviance explained = 95.3%
GCV = 4.3025  Scale est. = 4.2026    n = 12989

6. QGIS上での結果の表示

まず、プロジェクトに座標系を与える。新規プロジェクトを立ち上げてから、下記のように設定する。

「メニュー」→「プロジェクト」→「プロパティ」

次に、推定結果をQGIS上で表示する。一例として、上記の出力結果の「Kameoka-shi_LST_NDVI_est_225.txt」をQGIS上に取り込でいく。自動的に、各列の名前(フィールド名)を読み込んでくれる。

注意:Rで出力されたファイルはカンマ区切りのCSVファイルではなく、スペースの長さが一定ではない空白区切りのファイルである。そのため、以下の設定が必要になる。

下記のフィールド名でデータが読み込まれるはずである。

次に、背景画像と重ね合わせて表示させる。好きな背景を選んで構わない。

「メニューバー」→「レイヤ」→「XYZ接続」→「新規」

最後に、上記で作成した「Kameoka-shi_LST_NDVI_est_225.txt」を表示させ、標準誤差に応じて表示させる色を変えていく。

  「レイヤパネル」中の「Kameoka-shi_LST_NDVI_est_225」で右クリック
→「プロパティ」
→「シンボロジ」
亀岡市の時におけるLST推定の標準誤差:(左)DOY=145, (中)DOY=225, (右)DOY=305(2012年1月1日をDOY=1)

7. 課題

各班で以下の課題に取り組んで、レポートと発表資料を作成し、PandA上の班長のフォルダに提出する。提出後はメールにて連絡する。また講義内の発表会で成果を発表する。

  1. 今回の例では京都府内の複数都市を選定し、LSTを推定する際の標準誤差の分布図を作成した。各班で他の都道府県から3〜4つの市町村を選び(但し点数が少ないとRでの解析が正常にできないことに注意)、標準誤差の時空間分布の傾向や違いを考察せよ。
  2. 本来は平年に比べて異常にLSTが高い地点を抽出しようと試みたが、今回は標準誤差の表示に留まった。標準誤差自体は推定の良否に関連する指標であり、LSTの異常値を表すわけではない。各班で、標準誤差に代わる指標を使って、上記の課題と同様に3〜4つの市町村における異常値の時空間分布の傾向や違い、また標準誤差との違いを考察せよ。

[GRASSのインストール、標高データを用いた地滑り危険度マップの作成]
[植生指数 (NDVI) の計算、表示]
[標高データ (SRTM)の表示、植生指数 (NDVI) の3次元表示]
[反射率、輝度温度、標高データを用いた土地被覆分類]
[QGIS, Rを用いた公示地価データの空間統計分析:空間的自己回帰モデル]
[QGIS, Rを用いた公示地価データの空間統計分析:静的な時空間モデリング]
[Rを用いた衛星データ(LST, NDVI)の空間統計分析]
[Pythonを用いた空間統計分析のための衛星データ(LST, NDVI)の処理1(市町村別データの生成)]
[Pythonを用いた空間統計分析のための衛星データ(LST, NDVI)の処理2(人口データの重みを加味した市町村別データの生成)]

[講義「リモートセンシングと地理情報システム(GIS)」のトップページ]
[須ア純一のページへ]

須ア純一 京都大学大学院 工学研究科社会基盤工学専攻 空間情報学講座