rgisgeospatialspatialr-sp

Converting the Coordinate Reference System (CRS) to Create a Spatial Dataframe using SpatialPoints() in the SP Package


Issue:

I have a shapefile that I've imported into R and I selected the variables of interest for the analysis that I'm conducting. My ultimate aim is to interpolate the point data (dolphin IDs) to obtain sea surface temperature (SST) values from each individual raster file within a stack of 70+ rasters from an object called ncin_SST. This object was created using the function stack::raster() from multiple Aqua Modis netCDF files contained in one folder and I downloaded the data from the 'Ocean Color Project' from NASA. My goal is to extract the average SST per ID over the time period across all 70+ raster files from 2016-to 2021.

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

I imported my shapefile using the function shapefile() in the rgdal package and I want to extract three variables of interest involving variable 1 = ID, variable 2 = LONGITUDE, and variable 3 = LATITUDE. Before I can interpolate the IDs with the object ncin_SST, I need to produce a spatial data frame from coordinates using the SpatialPoints() function from the sp package.

When I am attempting to transform the CRS to WG85/UTM 34N, I am getting this error message below.

If anyone is able to help, then many thanks.

R-code

#Read our shapefile using the function shapefile() from the raster package

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

Check variable headers

head(Points_shp)

Select the variables of interest

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]]] 

Finalising the spatial coordinate reference data frame

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

Error Message

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

Solution

  • #Read our shapefile using the function shapefile() from the raster package

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

    Check variable headers

    head(Points_shp)
    

    Select the variables of interest from the 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) <- ~Longitude + Latitude 
    
    #Finalise the coordinate reference system dataframe
    points_spdf = SpatialPoints(coords_N, proj4string = crs("+proj=utm + zone=34 + datum=WGS84"))