R/GIS：用其他来源的数据覆盖 GADM 数据？

Posted 2023-02-25

【中文标题】R/GIS：用其他来源的数据覆盖 GADM 数据？

【英文标题】：R/GIS: Overlay GADM data with data from other sources?

【发布时间】：2012-02-09 12:29:27

【问题描述】：

我有来自Global Administrative Areas Database 的行政边界数据。使用这些数据绘制细分非常容易。例如，如果我想对佛罗里达州县的文件mydata 中包含的地理编码数据进行重叠绘制，我下载美国的“2 级”，然后使用以下代码：

load("/home/anindya/Desktop/DELETE/USA_adm2.RData")
temp = as.data.frame(gadm) *Code for eyeballing structure of data
florida = gadm[344:410,]
plot(florida); points(mydata$longitude, mydata$latitude)

但是你如何在上面叠加其他信息，比如地形和水道？我对来自Natural Earth 和Global Lakes and Wetlands Database 的数据特别感兴趣。

后续编辑：非常感谢 Paul 解决了这个问题。我在here at gis.stackexchange 上交叉发布了这个问题的一个变体，R.K.给出了很好的答案。也请阅读它。

【问题讨论】：

看看gis.stackexchange.com

抱歉，我找不到直接相关的内容。当然应该有一个简单的函数来覆盖同一地理区域上具有不同类型细节的两个地块？！

【参考方案1】：

如果我理解正确，您问的是如何叠加不同的地理图层，就像使用 GUI GIS 一样。

如果数据位于 Spatial* 对象中（例如 sp 包中的 SpatialPoints, SpatialLines, SpatialPolygons 等）或 raster 包中的 Raster* 对象，这很简单。如果给定参数add=TRUE，这两个包中的plot()方法都可以处理这种叠加。

要将基于矢量的地理数据导入Spatial* 对象，您可以使用maptools 包中的函数（例如readShapePoly() 读取ESRI 多边形形状文件）。可以使用raster() 加载来自各种文件规范的栅格，如果安装了rgdal 包，则可以使用其他格式。

这是一个使用伪造的地理数据叠加地理数据类型的示例。首先在 UTM 网格上创建一些示例数据：

library(sp)
library(raster)

## Create a RasterLayer object and fill with random values
baseRaster <- raster(nrow=100, ncol=100, 
    xmn=0, xmx=100, ymn=0, ymx=100,
    crs= "+proj=utm")
baseRaster[] <- runif(ncell(baseRaster))

## Create a second raster to superimpose
## It should contain NA values where it is transparent
supRaster <- baseRaster
supRaster[] <- NA
supRaster[cellFromCol(supRaster, 48:52)] <- 1


## Create SpatialPoints object to superimpose on these
loc <- SpatialPoints(cbind(seq(10, 90, by=10), seq(10, 90, by=10)))

现在做图：

## Plot base raster
plot(baseRaster)

## Superimpose second raster in a different colour
## Turn off legend
plot(supRaster, add=TRUE, col="blue", legend=FALSE)

## Superimpose points and make them big and colourful
plot(loc, add=TRUE, pch=20, cex=3, col="red")

结果如下：

【讨论】：

请注意，只有当数据集具有相同的坐标参考系时，这样的叠加才有效。 R 不会即时进行坐标转换，因此您必须使用 rgdal 中的 spTransform，例如，如果您在经纬度坐标中有一组数据，而在某个 UTM 度量网格中有一组数据。

@Spacedman - 转换为通用 CRS 通常是必不可少的，但在 R 中进行转换对我来说仍然有点神秘。我已经在***.com/questions/9230318

提出了这个后续问题