GDAL commands

備忘録

1992-2015年の一年ごとの土地被覆を300m解像度でしめしたESA-CCIデータ(約1.1GB)を、土地被覆クラスをメタクラスにまとめ、解像度さげ、パーセント表示する作業。
はじめはgdal_calc.pyとgal_warpでいけるかとおもってためしてみたが一向に終わらず・・・
次にお得意のRのraster::calcでやってみるが途中エラーが出てめんどくさくなって断念(←このやり方ではおそらくできる。がRAMを使わずtmpフォルダで処理をするため遅い)。

ということで、最も早かった方法は、gdal_translateで元ファイルを分割(1°メッシュごと)し、gdal_calc.pyでクラス編集し、gal_warpで解像度を下げ、最後にgdalbuildvrtでvrtファイルを作成する。丸2日にて完成

try user-defined function approach


this does not necessary to prepare large amount of RAM.

```{r}
na_func <- function(x){
x[x <= 0.5] <- NA
return(as.vector(x))
}

stk <- stack(r1,r2,r3,r4,r5,r6)
system.time(ans <- calc(stk, na_func))
summary(ans)
```

user system elapsed
10.585 0.254 10.898



#### calculation for temporal interpolating NA values in raster stack.


very nice function by https://stat.ethz.ch/pipermail/r-sig-geo/2010-December/010216.html
It is always a good way to use calc function, especially when using big raster data.
this do not depend on big RAM space.

```{r}
require(zoo)


fill.NA <- function(x,na.rm = F){
#x[x==0] <- NA
if(length(which(is.na(x)==F))>2){
fill <- approxfun(x, y=NULL, method='linear', rule=2)
rec <- which(is.na(x))
x[rec] <- fill(rec)
}

return(x)
}

stk_fillNA <- calc(stk, fun=fill.NA)

summary(stk_fillNA)

```

try normal approach

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 <= 0.5 ] <- NA)
summary(stk)
```

```{r}
system.time(stk[stk <= 0.5 ] <- NA)
summary(stk)
```

user system elapsed
0.984 0.406 1.403


set up raster stack

```{r}
library(raster)
r <- raster(ncols=3000, nrows=3000)
r1 <- setValues(r, runif(ncell(r)))
r2 <- setValues(r, runif(ncell(r)))
r3 <- setValues(r, runif(ncell(r)))
r4 <- setValues(r, runif(ncell(r)))
r5 <- setValues(r, runif(ncell(r)))
r6 <- setValues(r, runif(ncell(r)))
stk <- stack(r1,r2,r3,r4,r5,r6)

```

try user-defined function approach


this does not necessary to prepare large amount of RAM.

```{r}
na_func <- function(x){
x[x <= 0.5] <- NA
return(as.vector(x))
}

stk <- stack(r1,r2,r3,r4,r5,r6)
system.time(ans <- calc(stk, na_func))
summary(ans)
```

user system elapsed
10.585 0.254 10.898



#### calculation for temporal interpolating NA values in raster stack.


very nice function by https://stat.ethz.ch/pipermail/r-sig-geo/2010-December/010216.html
It is always a good way to use calc function, especially when using big raster data.
this do not depend on big RAM space.

```{r}
require(zoo)


fill.NA <- function(x,na.rm = F){
#x[x==0] <- NA
if(length(which(is.na(x)==F))>2){
fill <- approxfun(x, y=NULL, method='linear', rule=2)
rec <- which(is.na(x))
x[rec] <- fill(rec)
}

return(x)
}

stk_fillNA <- calc(stk, fun=fill.NA)

summary(stk_fillNA)

```

try normal approach

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 <= 0.5 ] <- NA)
summary(stk)
```

```{r}
system.time(stk[stk <= 0.5 ] <- NA)
summary(stk)
```

user system elapsed
0.984 0.406 1.403