问题概述:已有某地区站点数据,只有站点名称以及相应变量;另有较多站点数据,有站点名称,站号,经纬度信息,需要对该地区站点数据匹配站号以及经纬度信息,然后存为shp,最后叠加在行政区上:
实操借助chatgpt,个人感觉新必应在写代码方面不如chatgpt全面;;;但在搜寻数据集方面新必应更为精准;
代码完全借助chatgpt查找,出错的地方再次返回聊天框解决,大部分可解决:
导出的点shp文件用ARCGIS打开没问题,但在R中和行政矢量叠加绘图时有问题,点没出现,目前没解决,后续解决再更新~~~~
###站点信息整合匹配
library(openxlsx)
phesmaize<-read.xlsx("E:/seasonyield_predict/datapre/yield/phe_sprmaize.xlsx",sheet=1)
###将phesmaize和sprm相同属性值列转化统一名称,便于后续操作
colnames(phesmaize)[colnames(phesmaize)=="Sprm$NAME"] <- "NAME"
sprm<-read.xlsx("E:/seasonyield_predict/datapre/yield/Sprm.xlsx",sheet=1)
#以左侧为基础,按NAME列合并
merged_data <- merge(phesmaize, sprm, by = "NAME")
matched_data <- subset(merged_data, select = c(NAME,ID,LAT,LON))
#去除重复数据
sta_m<-unique(matched_data)
write.xlsx(sta_m,file=paste0("E:/seasonyield_predict/datapre/yield/stations/crop_station/",'sta_m','.xlsx'),colNames = TRUE)###存shp并叠加绘图
library(rgdal)
library(sp)
library(ggplot2)
points <- data.frame(longitude = sta_m$LON, latitude = sta_m$LAT)###转对象时的出错语句###
【#coordinates(points) <- c("longitude", "latitude")#创建一个包含点经纬度信息的数据框,
#假设经度信息保存在longitude列中,纬度信息保存在latitude列中:
#proj4string(points) <- CRS("+proj=longlat +datum=WGS84")#将点经纬度信息转换为SpatialPoints对象,并指定坐标系
#将SpatialPoints对象转换为SpatialPointsDataFrame对象,
#points_spdf <- SpatialPointsDataFrame(points, data.frame(ID = 1:nrow(points)), proj4string = CRS("+proj=longlat +datum=WGS84"))】
#Error in 1:nrow(points) : argument of length 0
In addition: Warning message:
In SpatialPointsDataFrame(points, data.frame(ID = 1:nrow(points)), :If the coords argument is a SpatialPoints object, set its CRS first;the proj4string argument to this function is ignored.# 报错信息重新输入对话框后,解决:::假设经度列名为"lon",纬度列名为"lat"
points_spdf <- SpatialPointsDataFrame(coords = points[, c("longitude", "latitude")], data = data.frame(ID = 1:nrow(points)), proj4string = CRS("+proj=longlat +datum=WGS84"))
#将其保存为shp文件:
writeOGR(points_spdf, dsn = "E:/seasonyield_predict/datapre/yield/stations/crop_station/", layer = "sta_m", driver = "ESRI Shapefile")# 读取已有的行政区矢量图
db <- readOGR(dsn = "E:/seasonyield_predict/datapre/region/roi/", layer = "roi_db")
# 绘制行政区矢量图
ggplot() + geom_polygon(data = db, aes(x = long, y = lat, group = group), fill = "gray80", color = "gray60", size = 0.2) +# 绘制点图层geom_point(data = points, aes(x = longitude, y = latitude), color = "red", size = 3)