如何将以网格间隔收集的点数据转换为 r 中的地理参考数据集？

Posted 2023-02-23

【中文标题】如何将以网格间隔收集的点数据转换为 r 中的地理参考数据集？

【英文标题】：How to convert point data collected at grid interval to a georeferenced dataset in r?

【发布时间】：2015-05-26 23:50:21

【问题描述】：

我有这个数据集：https://www.dropbox.com/s/k06n9l05t25r6x2/newdata.csv?dl=0

（示例）

"","row","col","flagrv"
"1",2361,530,2
"2",2378,531,2
"3",2360,531,2
"4",2355,531,2
"5",2363,532,2
"6",2359,532,2
"7",2368,533,2
"8",2367,533,2
"10",2359,533,2

如果我使用此代码进行绘图：

gs.pal <- colorRampPalette(c("blue", "green","yellow","orange","red"),bias=1,space="rgb")
ggplot(data=ndata,aes(x=col,y=row,color=flagrv)) + 
  geom_point(size = 0.01)+
  scale_colour_gradientn(name = "Scale",colours = gs.pal(5))+
  xlab('Longitude')+
  ylab('Latitude')+
  theme_bw()+
  theme(line = element_blank())+
  theme(legend.position = c(.93,.20),panel.grid.major = element_line(colour = "#854440"))+
  ggsave("test.png",width=10, height=8,dpi=300)

我们得到这个数字：

现在，问题是我没有 Lat-Long 值。我想覆盖州边界但不能使用 Maps 包。有人建议我使用 gdal，但我不知道如何。您能否告诉我如何将其映射到 Lat-Long 域，以便我可以轻松地对其进行操作。

编辑：

我从别人那里了解到我可以使用这个：

gdal_translate -a_srs EPSG:4269 FILE.asc FILE.tif

#

答案 1 的错误

Error: unexpected ']' in "spdf = SpatialPointsDataFrame(coords, all_data[, c("flagrv"]"

然后我把代码改成：

spdf = SpatialPointsDataFrame(coords, all_data[, c("flagrv")]) 

但现在我有这个错误：

Error in validObject(.Object) : invalid class “SpatialPointsDataFrame” object: invalid object for slot "data" in class "SpatialPointsDataFrame": got class "integer", should be or extend class "data.frame"

【参考方案1】：

在至少不知道数据集的投影和基准面的情况下（但希望有更多信息，例如分辨率和范围），没有简单的方法可以做到这一点。如果这是一个派生地图，请尝试查找用于生成它的内容。 有了这些信息，您就可以使用 raster 包中的投影函数来定义数据集的投影。

编辑（根据提供的其他信息，有一个可行的解决方案）： 这是一个可行的解决方案，因为数据集的左下角坐标为 24.55，-130，行/列之间的间距为 0.01 度，投影为 nad83。请注意，提供的元数据信息是错误的，因为最小纬度值不是 20 度，但可以从最南端（基韦斯特）估计为 24.55。

#load dataset 
all_data=(read.csv('new_data.csv',header=T, stringsAsFactors=F))
res=0.01 #spacing of row and col coords pre-specified
#origin_col_row=c(0, 0) 
origin_lat_lon=c(24.55, -130) 
all_data$row=(all_data$row)*res+origin_lat_lon[1] 
all_data$col=(all_data$col)*res+origin_lat_lon[2]

#now that we have real lat/lon, we can just create a spatial dataframe
library(rgdal)
library(sp)
coords = cbind(all_data$col, all_data$row)
spdf = SpatialPointsDataFrame(coords, data=all_data) #sp = SpatialPoints(coords)
proj4string(spdf) <- CRS("+init=epsg:4269") 

r 似乎在尝试绘制那么多点时窒息，所以为了检查答案是否有意义，我将数据集保存为 shapefile 并将其绘制在 arcgis 上：

writeOGR(spdf,"D:/tmp_shapefile4.shp", "flagrv", driver="ESRI Shapefile")

我设法使用 ggplot2 和下面的代码来绘制它，请耐心等待，因为绘制它需要一段时间：

df=as.data.frame(spdf)
library(ggplot2)
ggplot(data=df,aes(x=col,y=row,color=flagrv))+ 
   geom_point(size = 0.01)+
  xlab('Longitude')+
  ylab('Latitude')

【讨论】：

我收到了这个信息：如果有意义的话 - gdal_translate -a_srs EPSG:4269 FILE.asc FILE.tif

很好：这意味着，基于spatialreference.org，图像具有 NAD83 纬度/经度投影。问题是您的 X 和 Y 显然不是纬度/经度。如果不将 csv 文件中的 row 和 col 列转换为等效的 lat lon 值，我看不到这样做的明确方法。如果你这样做了，你就可以把你的数据变成一个空间数据框：检查maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/cheatsheet.html中的点和坐标部分@

如何将 X/Y 转换为等效的经纬度？你能帮忙写一些代码吗？我正在阅读该链接。如果我能做到这一点，我就可以做剩下的事情，比如覆盖状态边界。

这正是问题所在（不是代码问题）——您需要知道原始数据集的纬度/经度范围。如果丢失了，那么执行此操作的一种方法是通过查看 NAD83 地图中极端 N S E W 点的坐标来找到这些值。根据行数（对于 E/W）或列数（对于 N/S），您将获得将当前值转换为纬度/经度值所需的信息。同样，如果不从源文件中恢复此范围信息，我认为没有简单的方法解决它（但我可能错了）。

嗨，我也有这个信息：ncols 7000 nrows 3500 xllcorner -130 yllcorner 20 cellsize 0.01 NODATA_value -9999