使用 SP 包中的 SpatialPoints() 转换坐标参考系 (CRS) 以创建空间数据框

Question

问题：

我有一个 shapefile，已将其导入到 R 中，并为我正在进行的分析选择了感兴趣的变量。我的最终目标是插值点数据（海豚 ID），以从名为 ncin_SST 的对象的 70 多个栅格堆栈中的每个单独的栅格文件中获取 海面温度 (SST) 值。该对象是使用函数 stack::raster() 从一个文件夹中包含的多个Aqua Modis netCDF 文件 创建的，我从 NASA 的“海洋颜色项目”下载了数据。我的目标是提取 2016 年至 2021 年期间所有 70 多个栅格文件中每个 ID 的平均 SST。

##Stack all the netCDF files and select the variable "sst" from the raster layers
ncin_SST <- raster::stack(filenames, varname = "sst")

我使用

rgdal 包中的函数 shapefile() 导入了 shapefile，我想提取感兴趣的三个变量，涉及变量 1 = ID、变量 2 = 经度和变量 3 = 纬度。在使用 object ncin_SST 插入 ID 之前，我需要使用 sp 包中的 SpatialPoints() 函数从坐标生成空间数据框。 当我尝试将 CRS 转换为 WG85/UTM 34N 时，我收到以下错误消息。

如果有人能够提供帮助，那么非常感谢。

R 代码

#使用光栅包中的函数 shapefile() 读取我们的 shapefile

Points_shp <- shapefile(".", point_ID.shp")

检查变量头

head(Points_shp)

选择感兴趣的变量

coords = cbind(Points_shp$ID, Points_shp$LATITUDE, Points_shp$LONGITUDE)

#Making a spatial data frame from coordinates
#The IDs were documented in WG84/UTM 34N
#Extract the project code for the CRS


CRS("+init=epsg:32634")

结果

Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["WGS 84 / UTM zone 34N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 34N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",21,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]],
        ID["EPSG",16034]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
    USAGE[
        SCOPE["unknown"],
        AREA["Between 18°E and 24°E, northern hemisphere between equator and 84°N, onshore and offshore. Albania. Belarus. Bosnia and Herzegovina. Bulgaria. Central African Republic. Chad. Croatia. Democratic Republic of the Congo (Zaire). Estonia. Finland. Greece. Hungary. Italy. Kosovo. Latvia. Libya. Lithuania. Montenegro. North Macedonia. Norway, including Svalbard and Bjornoys. Poland. Romania. Russian Federation. Serbia. Slovakia. Sudan. Sweden. Ukraine."],
        BBOX[0,18,84,24]]] 

最终确定空间坐标参考数据框

points_spdf = SpatialPoints(coords, proj4string = crs("+proj=utm + zone=34 + datum=WGS84 + units=m + no_defs"))

错误信息

Error in .local(obj, ...) : 
  cannot derive coordinates from non-numeric matrix

Answer 1

Points_shp <- shapefile(".", point_ID.shp")

检查变量头

head(Points_shp)

从形状文件中选择感兴趣的变量

#Check the number of layers in Point_shp
dim(Points_shp)
#[1] 635  19

#Check the variable names in the Point_shp object
names(Points_shp)

#making a spatial data frame from coordinates
#EPSG = 32634
CRS("+init=epsg:32634")
coords = cbind(Points_shp$ID, Points_shp$LATITUDE, Points_shp$LONGITUDE)

#Check the variables in the coords object
print(coords)

#check header
head(coords)

#Check the structure of the matrix coords
str(coords)

#Change the matrix coords into a dataframe
coords_N<-as.data.frame(coords)

#Rename Column Names of coords_N dataframe
#Cbind transformed the column headings to V1, V2, and V3
#Change the column names
colnames(coords_N)[colnames(coords_N) == "V1"] <- "ID"
colnames(coords_N)[colnames(coords_N) == "V2"] <- "Longitude"
colnames(coords_N)[colnames(coords_N) == "V3"] <- "Latitude"

#Check the structure of the data frame coords
str(coords_N)

#Change the format of the variables from characters to numerica format
coords_N$ID <- as.numeric(as.character(coords_N$ID))
coords_N$Longitude <- as.numeric(as.character(coords_N$Longitude))
coords_N$Latitude <- as.numeric(as.character(coords_N$Latitude))

#Change latitude and longitude into coordinates using the sp package
coordinates(coords_N) <- ~Longitude + Latitude 

#Finalise the coordinate reference system dataframe
points_spdf = SpatialPoints(coords_N, proj4string = crs("+proj=utm + zone=34 + datum=WGS84"))