2017-07-11 1 views
1

1つの州の土壌データのラスタを取得し、郡ごとにトリミングし、各郡のセル値を郡fipsコードに変更してから郡ラスタを再び状態ラスタにマージし直します。オーバーラップエクステントを使用してラスタの切り抜き、値の変更、マージ

ここでは、土壌ラスタ(セルの値として各土壌タイプに関連付けられたマップ単位キーとしてのデフォルト)と米国郡のポリゴンを読み込みます。私は1つの州のちょうどのポリゴンを選択し、それを土壌ラスタと同じコーディアンテシステムに変換してから土壌のラスタとポリゴンの2つの郡を選択します。

state_soils_raster <- raster("MapunitRaster_IL_10m.tif") 
us_county_polygons <- readOGR("cb_2016_us_county_500k/cb_2016_us_county_500k.shp") 

IL_county_polygons <- us_county_polygons[us_county_polygons$STATEFP == 17,] 
IL_county_polygons <- spTransform(IL_county_polygons, CRS = crs(state_soils_raster)) 

county1 <- "Douglas" 
county2 <- "Coles" 

county1_polygon <- IL_county_polygons[IL_county_polygons$NAME %in% county1,] 
county2_polygon <- IL_county_polygons[IL_county_polygons$NAME %in% county2,] 

county1_raster <- crop(state_soils_raster, county1_polygon) 
county2_raster <- crop(state_soils_raster, county2_polygon) 

Iは、単独で各郡をプロットする場合は、トリミングされた領域の範囲が矩形であり、郡自体の領域を越えて延びていることがわかります。 mukeyの値はどこにでもあるので(典型的には郡によってグループ化されていますが)、着色は狂っています。郡1は郡2の北に位置しています。

plot(county1_raster) 
plot(county1_polygon, add = T) 

enter image description here

plot(county2_raster) 
plot(county2_polygon, add = T) 

enter image description here

私があるような値のままにして2つの郡のラスタをバックマージする場合は、すべてが正常です。 2つのラスタの範囲が重なっていても、どのラスタmergeが引っ張られているかにかかわらず、セルの値は同じです。私は実際にラスターmergeがこの場合に引っ張られるのかどうかはわかりませんが、それは本当に重要ではありません。すべてがうまく組み合わさり、セルの値が正しい。

both_counties_raster <- merge(county1_raster, county2_raster) 
plot(both_counties_raster) 
plot(county1_polygon, add = T) 
plot(county2_polygon, add = T) 

enter image description here

しかし、私は何をしたい郡ラスタを再結合する前に、郡でセル値にを変更することです。

values(county1_raster) <- 1 
values(county2_raster) <- 2 
both_counties_raster_new <- merge(county1_raster, county2_raster) 

すべてがうまくマージしますが、私は今、新しい組み合わせラスタをプロットするときには、両方の郡ラスタmergeに含まれた細胞のためだけのラスタのうちの一つからセルの値を取ったことは明らかです。明らかにmergeはデフォルトで最初の入力ラスタに優先順位を付けます。

plot(both_counties_raster_new) 
plot(county1_polygon, add = T) 
plot(county2_polygon, add = T) 

enter image description here

私が探しているだけで、各郡の境界内にセルの値を変更してから、再び一緒にすべての郡をマージすることです。

私は(hereを説明した)10メートルのセルの分解能で、NAに郡の境界の外に何かを変えることができますraster::mask機能の認識しています、これは時間の非常識な量を取ります!

また、raster::rasterize関数を使用して郡の境界ポリゴンを、州の土壌ラスターの同じセルサイズと範囲を持つラスターに変換する代替アプローチを試しました。繰り返しますが、10mのセル分解能では、これは永遠にかかります。 1.5時間で私の8つのコアごとに1つの郡を処理することができました。そして、私は全国にやるべきだ!

誰かが私にそれを指摘すれば、それは素晴らしいだろうが、私は10mのラスター米国郡データセットに気づいていない。

土壌データはgSSURGOデータです - gSSURGOには多くのテーブルに郡属性があるかどうかもわかりません。それがあれば、私はそれを見つけることができません。それは簡単な解決策でもあります。

答えて

0

raster::cellFromPolygonで試したことがありますか?
ここに簡単な例を示します。

# Create a raster with zero values 
r <- raster(ncols=30, nrows=30, res = 1/3) 
values(r) <- 0 
# Create polygons 
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)) 
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)) 
pols <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), Polygons(list(Polygon(cds2)), 2))) 
plot(r) 
plot(pols, add = TRUE) 

r2 <- r 
# Find which cells are in which polygons 
cellpol <- cellFromPolygon(r, pols) 
# Not a really clean way to attribute values in the global environment... 
lapply(1:length(cellpol), function(x) values(r2)[cellpol[[x]]] <<- x) 
plot(r2) 
plot(pols, add = TRUE) 
関連する問題