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

Issue:

我有一个形状文件我已将其导入到 R 中，并为正在进行的分析选择了感兴趣的变量。我的最终目标是插值点数据（海豚 ID）获取海面温度（SST）堆栈中每个单独的光栅文件的值70+ 栅格来自名为 ncin_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")

我使用 shapefile() 中的函数导入了 shapefilergdal包我想提取三个变量感兴趣的涉及变量 1 = ID，变量 2 = 经度，变量 3 = 纬度。在我可以之前插入 ID与对象 ncin_SST，我需要使用以下坐标从坐标生成空间数据框空间点()函数从sp包.

当我尝试将 CRS 转换为 WG85/UTM 34N 时，收到以下错误消息。

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

R-code

#使用栅格包中的函数 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")

Results

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

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

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

检查变量头

head(Points_shp)

从 shapefile 中选择感兴趣的变量

#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) <- ~Latitude + Longitude

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