ggplotで複数プロット（スケール一緒）を描く

library(ggplot2)
library(reshape2)
tmp <- shp1@data[,c(13,19,24)]　##描きたい値（複数）
tmp.melt <- melt(tmp, id="id")　　　##idを使ってmelt. idがない場合はcbindでくっつけておく
tmp.ggmap <- fortify(shp1, region="id") ##shp1（描きたいベクターデータをggplot用に変換する）
tmp.ggmap <- join( tmp.ggmap, tmp.melt, by ="id")　##くっつける
PC1 <- ggplot( tmp.ggmap, aes(x = long, y = lat, group=group, fill=value)) +
geom_polygon() +
coord_map() +
facet_wrap(~variable) +　　##こいつが複数描く
scale_fill_gradient2(low="blue", high="red", midpoint=0) +　##色
theme(legend.position="right")


参考: http://www.cartopedia.co.uk/blog/2013/06/02/facet-wrapping-multivariate-data-reshape-and-ggplot/

全球土地被覆データMCD12Q1の作成　HDFダウンロード〜GeoTiffの抽出・マージ（モザイク）〜投影法の変換まで

シェルとGDALのみでいきます。
1.HDFダウンロード
wget -r -l 1 -A hdf http://e4ftl01.cr.usgs.gov/MOTA/MCD12Q1.051/2009.01.01/

2.　各HDFファイルからTIFFに変換。今回はLand Cover Type 1 (2001): IGBP global vegetation classification schemeのみをピックアップします。

以下のシェル（mcd12q1_hdf_to_tif.sh）を実行する（元ネタはどこだったっけ・・・）
#!/bin/bash
#Set field separator to linefeeds rather than spaces
echo "#!/bin/bash" > batch.sh
chmod +x batch.sh

for FILE in *hdf
do
echo $FILE
LAYERS=$(gdalinfo $FILE | grep SUBDATASET | grep NAME | grep Land_Cover_Type_1 | sed 's/[A-Z0-9\_]*=//g')

for LAYER in $LAYERS
do
LAYER=$(echo $LAYER | sed 's/ //g')
NEWFILE=$(echo $LAYER | \
sed 's/HDF4_EOS:EOS_GRID:"//g' | sed 's/.hdf":Grid:/_/g' | sed 's/ //g' | sed 's/"//g' | sed 's/ //g'| sed 's/:/./g').tif
CMD="gdal_translate -of GTiff '${LAYER}' $NEWFILE"
echo $CMD >> batch.sh
done
done

exec ./batch.sh

３．できたGeoTiffファイルをマージする

http://invisibleroads.com/tutorials/gdal-raster-merge.html

このシェル（gdal_tif_merge.tif）を実行
#!/bin/bash
#Set field separator to linefeeds rather than spaces
echo "#!/bin/bash" > batch.sh
chmod +x batch.sh
#IFS=$'\n'

LIST=`ls | grep Land_Cover_Type_1.tif`
LIST2=`echo $LIST | sed 's/\n/ /g'`
gdal_merge.py $LIST2 -o MOD12Q1_global.tif



４，投影法の変換もGDALで。
WGS84（EPSG:4326）へ変換する。

gdalwarp -t_srs EPSG:4326 -of GTiff MOD12Q1_global.tif MOD12Q1_global_wgs84.tif

netCDFデータのインポート（R)

色々とパッケージがあるそうですが、、raster libraryを使えば簡単だということがわかりました。

ただしncdf4パッケージが必要です。
rstr.brick <- brick("AAA.nc", vername = "NAME")
NAMEはgdalinfoで確認できる。
以上！

spplotのTips

levelplotと似てるとこもあるけど違うとこも。
よくごっちゃになる。。


#タイトルの追加
spplot( raster, main = list(label="title", cex = 3))


#色の変更
spplot( raster, col.regions=colorRampPalette(brewer.pal(9, 'Reds')))

#値の区切り
spplot( raster, col.regions=colorRampPalette(brewer.pal(9, 'Reds')), at=c(0.5,1.5,2,5))


#凡例を消す
spplot( raster, colorkey=F)

#凡例の位置を下にする（ついでにレベルの文字サイズを大きくするなど）
colorkey=list(space="bottom", tick.number=2, tck=2, labels=list(cex=1.5))


spplot( raster, col.regions=colorRampPalette(brewer.pal(9, 'Reds')), at=c(0.5,1.5,2,5),
colorkey=list(space="bottom", tick.number=2, tck=2, labels=list(cex=1.5)))


#ベクターを重ねる
spplot(raster) + layer(sp.polygons(bnd))
もしくは
layer=list("sp.polypons", shpfile, fill="gray70", lwd=0)
spplot(raster,
sp.layout=list(layer))

#緯度経度をつける
scales=list(tick.number=2, draw=TRUE, tck=0.5, cex=1.2)

spplot( raster, col.regions=colorRampPalette(brewer.pal(9, 'Reds')), at=c(0.5,1.5,2,5),
colorkey=list(space="bottom", tick.number=2, tck=2, labels=list(cex=1.5)),
scales=list(tick.number=2, draw=TRUE, tck=0.5, cex=1.2))