如何自动将多边形形状文件的许多字段转换为 R 中的栅格
How to automatically convert many fields of a polygon shapefile to raster in R
我有一个表示泰森多边形的 shapefile。
每个多边形都与 table 的许多值相关联。
thiessen <- readOGR(dsn = getwd(), layer = poly)
OGR data source with driver: ESRI Shapefile
Source: ".../raingauges/shp", layer: "thiessen_pol"
with 10 features
It has 5 fields
head(thiessen)
est est_name p001 p002 p003
0 2 borges 1 8 2
1 0 e018 2 4 3
2 5 starosa 5 15 1
3 6 delfim 4 2 2
4 1 e087 1 1 3
5 3 e010 0 1 0
列'est'和'est_name'与雨量计的ID和名称有关。以下几列对我很重要,代表第 1 天、第 2 天等的降水值(在示例中我只保留了三天,但实际上,我有 8 年的每日降水数据)。
我需要将多边形转换为光栅,但 table.
的每个字段(列 p001、p002 等)需要一个光栅
有一种简单的方法可以使用 R 中的函数 rasterize 将多边形转换为栅格。
r_p001 <- rasterize(thiessen, r, field = thiessen$p001)
plot(r_p001)
writeRaster(r_p001, filename=".../raingauges/shp/r_p001.tif")
问题是我需要手动设置 table 的字段(列),其中包含要转换为栅格的多边形值。因为我有大约 2900 天(每个雨量计有 2900 列降水值),所以不可能手动完成。
文档没有帮助阐明如何自动执行此过程,我在互联网上没有找到任何帮助我的东西。
有谁知道如何自动将每个字段转换为栅格并另存为 tif 格式?
这里有一个方法:
示例数据
library(raster)
r <- raster(ncols=36, nrows=18)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
att <- data.frame(id=1:3, var1=10:12, var2=c(6,9,6))
pols <- spPolygons(p1, p2, p3, attr=att)
重要的是要有一个唯一的字段如果你的数据没有它,像这样添加它
pols$id <- 1:nrow(pols)
栅格化
r <- rasterize(pols, r, field='id')
为所有其他变量创建层
x <- subs(r, data.frame(pols), by='id', which=2:ncol(pols), filename="rstr.grd")
x
#class : RasterBrick
#dimensions : 18, 36, 648, 2 (nrow, ncol, ncell, nlayers)
#resolution : 10, 10 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : rstr.grd
#names : var1, var2
#min values : 10, 6
#max values : 12, 9
另一种方法是保留一个具有光栅属性 Table 的图层,这样更快,但根据您的目的,可能是一种不太有用的方法:
r <- rasterize(pols, r, field='id')
f <- as.factor(r)
v <- levels(f)[[1]]
v <- cbind(v, data.frame(pols)[,-1])
levels(f) <- v
f
#class : RasterLayer
#dimensions : 18, 36, 648 (nrow, ncol, ncell)
#resolution : 10, 10 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names : layer
#values : 1, 3 (min, max)
#attributes :
# ID var1 var2
# 1 10 6
# 2 11 9
# 3 12 6
然后你可以这样做:
z <- deratify(f)
得到与第一个例子相同的结果
z
#class : RasterBrick
#dimensions : 18, 36, 648, 2 (nrow, ncol, ncell, nlayers)
#resolution : 10, 10 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names : var1, var2
#min values : 10, 6
#max values : 12, 9
我有一个表示泰森多边形的 shapefile。
每个多边形都与 table 的许多值相关联。
thiessen <- readOGR(dsn = getwd(), layer = poly)
OGR data source with driver: ESRI Shapefile
Source: ".../raingauges/shp", layer: "thiessen_pol"
with 10 features
It has 5 fields
head(thiessen)
est est_name p001 p002 p003
0 2 borges 1 8 2
1 0 e018 2 4 3
2 5 starosa 5 15 1
3 6 delfim 4 2 2
4 1 e087 1 1 3
5 3 e010 0 1 0
列'est'和'est_name'与雨量计的ID和名称有关。以下几列对我很重要,代表第 1 天、第 2 天等的降水值(在示例中我只保留了三天,但实际上,我有 8 年的每日降水数据)。
我需要将多边形转换为光栅,但 table.
的每个字段(列 p001、p002 等)需要一个光栅有一种简单的方法可以使用 R 中的函数 rasterize 将多边形转换为栅格。
r_p001 <- rasterize(thiessen, r, field = thiessen$p001)
plot(r_p001)
writeRaster(r_p001, filename=".../raingauges/shp/r_p001.tif")
问题是我需要手动设置 table 的字段(列),其中包含要转换为栅格的多边形值。因为我有大约 2900 天(每个雨量计有 2900 列降水值),所以不可能手动完成。
文档没有帮助阐明如何自动执行此过程,我在互联网上没有找到任何帮助我的东西。
有谁知道如何自动将每个字段转换为栅格并另存为 tif 格式?
这里有一个方法:
示例数据
library(raster)
r <- raster(ncols=36, nrows=18)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
att <- data.frame(id=1:3, var1=10:12, var2=c(6,9,6))
pols <- spPolygons(p1, p2, p3, attr=att)
重要的是要有一个唯一的字段如果你的数据没有它,像这样添加它
pols$id <- 1:nrow(pols)
栅格化
r <- rasterize(pols, r, field='id')
为所有其他变量创建层
x <- subs(r, data.frame(pols), by='id', which=2:ncol(pols), filename="rstr.grd")
x
#class : RasterBrick
#dimensions : 18, 36, 648, 2 (nrow, ncol, ncell, nlayers)
#resolution : 10, 10 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : rstr.grd
#names : var1, var2
#min values : 10, 6
#max values : 12, 9
另一种方法是保留一个具有光栅属性 Table 的图层,这样更快,但根据您的目的,可能是一种不太有用的方法:
r <- rasterize(pols, r, field='id')
f <- as.factor(r)
v <- levels(f)[[1]]
v <- cbind(v, data.frame(pols)[,-1])
levels(f) <- v
f
#class : RasterLayer
#dimensions : 18, 36, 648 (nrow, ncol, ncell)
#resolution : 10, 10 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names : layer
#values : 1, 3 (min, max)
#attributes :
# ID var1 var2
# 1 10 6
# 2 11 9
# 3 12 6
然后你可以这样做:
z <- deratify(f)
得到与第一个例子相同的结果
z
#class : RasterBrick
#dimensions : 18, 36, 648, 2 (nrow, ncol, ncell, nlayers)
#resolution : 10, 10 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory
#names : var1, var2
#min values : 10, 6
#max values : 12, 9