Import Packages
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.1.1, released 2020/06/22
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
Import shapefile of plot south west corner
SWC<-shapefile('/Users/lavran_pagano/Downloads/Creekshead_Plots (1)/Creekshead_Plots.shp')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum WGS_1984 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs
## Warning in showSRID(wkt2, "PROJ"): Discarded ellps WGS 84 in CRS definition:
## +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m
## +nadgrids=@null +wktext +no_defs
## Warning in showSRID(wkt2, "PROJ"): Discarded datum WGS_1984 in CRS definition
Get Lat-Long Cordinates from Plots
long_lat<-data.frame(Long = SWC$ESRIGNSS_1 ,Lat = SWC$ESRIGNSS_L )#get lat and long cordinates
long_lat<-na.omit(long_lat) # remove NAs
Create a function to covert to UTM projection
# convert to UTM FUnction
LongLatToUTM<-function(x,y,zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
Apply UTM Function
utm<-LongLatToUTM(long_lat$Long,long_lat$Lat,17)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown based on WGS84 ellipsoid in CRS definition
Create a dataframe file with plot the UTM plot center cordinates
utmEasting<-utm$X +5
utmNorthing<-utm$Y +5
utm_center <-data.frame(ID=utm$ID, Easting=utmEasting, Northing = utmNorthing)
Renane if not already in the appropriate format
centroids <- utm_center
colnames(centroids)[1]<- "plot"
colnames(centroids)[2]<-"Easting"
colnames(centroids)[3]<-"Northing"
Get cordinates of plot corners
radius <- 5
yPlus <- centroids$Northing+radius
xPlus <- centroids$Easting+radius
yMinus <- centroids$Northing-radius
xMinus <- centroids$Easting-radius
Calculate polygon coordinates for each plot centroid
square=cbind(xMinus,yPlus, # NW corner
xPlus, yPlus, # NE corner
xPlus,yMinus, # SE corner
xMinus,yMinus, # SW corner
xMinus,yPlus) # NW corner again - close ploygon
Extract the plot ID information
# Extract the plot ID information
ID=centroids$plot
Create spatial polygon dataframe
polys <- SpatialPolygons(mapply(function(poly, id)
{
xy <- matrix(poly, ncol=2, byrow=TRUE)
Polygons(list(Polygon(xy)), ID=id)
},
split(square, row(square)), ID),
proj4string=CRS(as.character('+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
## but +towgs84= values preserved
Create spatial polygon DataFrame and export as a shapefile
## create spatial polygon DataFrame and export as a shapefile
polys.df <-SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))
# write to a polygon shapefile
writeOGR(polys.df, '.', 'CreeksheadPlotPolygon2', 'ESRI Shapefile',overwrite=TRUE)