r - 创建方形网格并将其导出为 shapefile 或表格

Question

我想使用谷歌地球引擎来提取某些国家的数据。我需要方形网格形式的数据，所以我想为某个国家创建这些方形网格，将它们添加到 shapefile，然后将 shapefile 导入地球引擎。我已经找到了一些创建方形网格的代码（在 shapefile 中创建一个网格），但现在我有两个问题。

首先，我需要导出方形网格，以便将它们导入地球引擎。我对 shapefile 的替代品非常开放。

其次，后续代码适用于某些国家（如法国），但不适用于其他国家（如泰国）。

library(raster)
shp = getData(country = "FRA", level = 0)

shp = spTransform(shp, CRSobj = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
plot(shp)

cs = c(10000, 10000)
grdpts = makegrid(shp, cellsize = cs)

spgrd = SpatialPoints(grdpts, proj4string = CRS(proj4string(shp)))

spgrdWithin = SpatialPixels(spgrd[shp,])
plot(spgrdWithin, add = T)

在第 2 行中将“FRA”替换为“THA”会导致spTransform出现错误。

标签： rgeospatialshapefiler-rasterrasterizing

Answer 1

失败是因为您使用的是 utm zone 32。您需要根据国家/地区的经度使用该区域。你可以在这里看到它们

您可以使用自动查找区域ceiling((longitude+180)/6)

library(raster)
s <- getData(country = "FRA", level = 0)

获取质心。在这种情况下，你可以做

centr <- coordinates(s)

如果有多个多边形，你可以这样做

centr <- apply(coordinates(s), 2, mean)

计算 UTM 区域。（请注意，法国有 32 个，这不好）

zone <- ceiling((centr[1] + 180)/6)
zone
#[1] 31

然后像这样使用它

crs <- paste0("+proj=utm +datum=WGS84 +unit=m +zone=", zone)
st <- spTransform(s, crs)

对于泰国，你会得到

s <- getData(country = "THA", level = 0)
centr <- apply(coordinates(s), 2, mean)
zone <- ceiling((centr[1] + 180)/6)
zone
#[1] 47

然而，这并不是一种适用于所有国家的方法。UTM 区域宽 6 度，许多国家/地区跨越多个区域（俄罗斯以 28 个区域占据了蛋糕）。因此，根据您的目标，您可能需要使用另一个坐标参考系统 (crs)。

之后，获得正方形多边形的另一种方法是创建一个范围为 s 和分辨率选择的 RasterLayer。但我怀疑这是从 GEE 中获取数据的最佳方式。我建议改为上传国家大纲。

r <- raster(st, res=10000)
r <- rasterize(st, r, 1)
x <- as(r, "SpatialPolygons")

# write to file
shapefile(x, "test.shp")

# view
plot(x)