備忘録1992-2015年の一年ごとの土地被覆を300m解像度でしめしたESA-CCIデータ（約1.1GB）を、土地被覆クラスをメタクラスにまとめ、解像度さげ、パーセント表示する作業。 はじめはgdal_calc.pyとgal_warpでいけるかとおもってためしてみたが一向に終わらず…

Githubにて。 https://github.com/naru-T/Rscripts/blob/master/raster_calc_function.mdもうちょっとmarkdownの書き方べんきょうせねば。

this does not necessary to prepare large amount of RAM.```{r} na_func x[x return(as.vector(x)) }stk system.time(ans summary(ans) ``` user system elapsed 10.585 0.254 10.898 #### calculation for temporal interpolating NA values in raster st…

this approach depends on RAM. it is only faster when using smaller raster data. Do not try this way if you use big raster data.```{r} system.time(stk[stk summary(stk) ``````{r} system.time(stk[stk summary(stk) ``` user system elapsed 0.984…

```{r} library(raster) r r1 r2 r3 r4 r5 r6 stk ```

this does not necessary to prepare large amount of RAM.```{r} na_func x[x return(as.vector(x)) }stk system.time(ans summary(ans) ``` user system elapsed 10.585 0.254 10.898 #### calculation for temporal interpolating NA values in raster st…

this approach depends on RAM. it is only faster when using smaller raster data. Do not try this way if you use big raster data.```{r} system.time(stk[stk summary(stk) ``````{r} system.time(stk[stk summary(stk) ``` user system elapsed 0.984…

```{r} library(raster) r r1 r2 r3 r4 r5 r6 stk ```

こちらのネタです。マルチパイプもいけます https://twitter.com/R_Programming/status/713217963801710594 `%s%` これだけ。> "A" %s% "B"%s% "C" [1] "ABC"

Ubuntu, imagemajick上にて. convert -delay 50 `ls -v ABC_*.png` output.gif delayのあとの数字は1/100sec `ls -v ABC_*.png` ワイルドカードは数字を順番にならべてくれない（1, 10,11, 12, 2, 3,...)ので、 ls -vでならべる

library(ggplot2) library(reshape2) tmp tmp.melt tmp.ggmap ベクターデータをggplot用に変換する） tmp.ggmap PC1 value)) + geom_polygon() + coord_map() + facet_wrap(~variable) + ##こいつが複数描く scale_fill_gradient2(low="blue", high="red", m…

シェルと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のみをピ…

色々とパッケージがあるそうですが、、raster libraryを使えば簡単だということがわかりました。ただしncdf4パッケージが必要です。 rstr.brick NAMEはgdalinfoで確認できる。 以上！

levelplotと似てるとこもあるけど違うとこも。 よくごっちゃになる。。 #タイトルの追加 spplot( raster, main = list(label="title", cex = 3)) #色の変更 spplot( raster, col.regions=colorRampPalette(brewer.pal(9, 'Reds')))#値の区切り spplot( raste…

色の変更 →library(RColorBrewer)を使って col.regions=colorRampPalette(brewer.pal(9, "RdBu"))などを追加 colorRampPaletteをつけないと変な挙動になることあり（なんで？？）デフォルトであるLongtitude, Latitudeの表示を消す → xlab="", ylab="" を追…

以前過去記事（時系列ラスターデータの処理について）ではでndvitsを用いたけどもかなり古いし、ライブラリがアーカイブされてしまったので、新たな取組を。 １．複数枚のtifファイルを複数バンドのもつ一つのtifにしてstackとして読み込む。 過去記事（大量…

#STACK is a raster stackras=STACK1 iii fun summary(lm(x ~ year))$coefficients} A=apply(STACK[iii], 1, fun) ras[iii]=A[1,] STACK[iii].intercept=ras ras[iii]=A[2,] STACK[iii].slope=ras ras[iii]=A[8,] STACK[iii].pval=ras

Raster １バンドのみで構成されるメッシュデータ r=raster("パス/◯◯.tif") Raster stack 複数バンドで構成されるメッシュデータ。brickと違い、データがすべてメモリ上に格納されるわけではないので、大規模データはstackを利用するほうがパンクしない。 s1=…

ここを参照 https://oku.edu.mie-u.ac.jp/~okumura/stat/colors.html #色を指定 #YnGnBlの連続色をかく。Hex codeをチェック brewer.pal(9, "YnGnBl")#色の定義 cols = colorRamp(c("#FFFFD9", "#EDF8B1", "#C7E9B4", "#7FCDBB", "#41B6C4", "#1D91C0", "#22…

今までやってこなかったので。 GRASSのi.groupで簡単に。i.group group=NAME in=image1,image2,image3,..image1,image2,image3,..のところはg.mlistで片付ける。 g.mlist pattern="AAA*" exclude="AAA1" separator=, 一枚のラスタ‐に大量のバンドが入ってる…

QGISでDissolveをするとテーブル内がわけわからないことになる。 （ランダムに行が選択されている？） GRASSでやるとテーブル内が消える。 ArcGISでやると和をとるか平均をとるかなど選択できる。が、選択が異様に面倒くさい。

LandsatのDLにて。 wget –nc –r –l 3 –A tar.gz http://XXX.XXX.XXX wget: Download data from http://XXX.XXX.XXX nc: Skip existing file which already download in the target directory. r: recursive. l 3: Directory level. Download only 3 level fr…

ルートのパーティションとホームのパーティションをddコマンドでimg化してバックアップし、その後リストアするメモ。１．DVDかUSBでUBUNTUを立ち上げる。Try UBUNTUを選択 ２．バックアップするデバイスをマウント 2-1. マウント先の作成 sudo mkdir /media/…

http://m884.jp/RzTutorial_Ja.htmlGUIの起動 library(Rz)

ここを参考に。 元shpデータ内テーブルと結合したいテーブル内にマッチングさせるIDが必要です。まず、もとshpをインポート A_shp=readShapePoly("A.shp", IDvar="OBJECTID")テーブルの確認 A_shp@data summary(A_shp)結合したいテーブルデータをインポート …

ここを参考にしています。 QGISでもPoint Sampling Toolがあるけども、ラスターの数が多いとしんどいですので、Rでやってみました。library(rgdal) library(maptools) library(raster) #ラスタースタックの作成 #実際は１００枚ほどのラスターをスタックrasS…

こちらを参考にしました。HOME=$HOME/【インストール先ディレクトリ】 /usr/bin/dropbox start -iこんだけでオッケっす。

簡単なことでも忘れてしまうのでメモ。tar.sh - #!/bin/shfor i in `ls | grep AA` dotar czvf ${i}.tar.gz ${i} rm ${i} done -

この記事はFOSS4G Advent Calendar 2013に参加しています。 普段は自分用メモブログなので今回も雑になる可能性があります。。 GISはいうまでもなく空間データ処理にフォーカスがあたってますが、時系列ものにはめっぽう弱い印象です。唯一救いがあるのがGRA…

http://hiroshiykw.blogspot.jp/2008/01/blog-post_24.html困っていた１の次に１０がくる問題。↑記事により解決です。ls -v GIMMS_*{?,??}*.bin | uniq とすれば番号順に並ぶ。なぜかダブって表示されるのでuniqで対応。