Libraries used:
library(sp)
library(sf)
library(ggplot2)
library(ggmap)
Created a dataframe called "coordinate.data" with longitude & latitude as column names, & weather station locations as row names.
longitude <- c(-73.964482,-73.953678,-73.893522,-73.815856,-74.148499)
latitude <- c(40.767544,40.631762,40.872481,40.734335,40.604014)
coordinate.data <- data.frame(longitude,latitude)
rownames(coordinate.data) <- c("MANH","BKLN","BRON","QUEE","STAT")
I then retreived shapefile data of NJ counties and NYC boroughs, & deleted all unnecessary columns so only the geometry field was left in both shapefiles. The NYC Boroughs shapefile data was downloaded from NYC Open Data, & the NJ county boundaries was downloaded from NJGIN Open Data.
nj.shp <- st_read("~/Downloads/NJ/NJ_Counties.shp")
nj <- nj.shp[,-(1:21)]
nyc.shp <- st_read("~/Downloads/NY/NYC_Boroughs.shp")
nyc <- nyc.shp[,-(1:4)]
I formatted both shapefiles to have the same projection (ESPG code 3857) and combined them into a shapefile dataframe with 26 observations (counties/boroughs) in one variable (geometry).
same.projection <- CRS("+init=EPSG:3857")
nj.data <- st_transform(nj,same.projection)
new.projection <- CRS("+init=EPSG:3857")
nyc.data <- st_transform(nyc,new.projection)
combined.data <- rbind(nj.data,nyc.data)
I am now attempting to plot the combined shapefile ("combined.data") on a map, in addition to the weather station locations ("coordinate.data"). When I attempt this, it runs inevitably & R shuts down. If I remove geom_sf(...), it plots the stations & formats everything correctly, so I assume the issue is with this line of code.
mesonet.map <-ggplot() +
ggtitle("NY Mesonet Site Locations") +
xlab("Longitude") +
ylab("Latitude") +
geom_point(data=coordinate.data,aes(x=longitude,y=latitude))+
geom_text(aes(x=longitude,y=latitude,label=rownames(coordinate.data)),size=3.25,nudge_y=0.02)+
geom_sf(data=combined.data,fill='darkgreen') +
mesonet.map + theme(
panel.background=element_rect(fill="lightblue",color="lightblue",size=0.5,linetype="solid"),
panel.grid.major=element_line(size=0.5,linetype='solid',color="white"),
panel.grid.minor=element_line(size=0.25,linetype='solid',color="white")
)
I'm not too sure what you were deleting from those shape files. I didn't delete anything. Nor did I combine anything. There are four separate layers projected in my output: the base map, NY shape file, NJ shape file, and the shape file called sites that was created and described below.
This map was created using 3 existing shape files, and the site shape file that was created in the following steps. The variables that were created from each step are printed and shown with the descriptions of each step.
First, created the spatial geometry variable:
MULTIPOINT ((-73.96448 40.76754), (-73.95368 40.63176), (-73.89352
40.87248), (-73.81586 40.73434), (-74.1485 40.60401))
Then, created the geometry column (combined the multipoint variable with the crs)
MULTIPOINT ((-73.96448 40.76754), (-73.95368 40...
Geometry set for 1 feature
geometry type: MULTIPOINT
dimension: XY
bbox: xmin: -74.1485 ymin: 40.60401 xmax: -73.81586 ymax: 40.87248
projected CRS: WGS 84 / Pseudo-Mercator
Then created a regular df consisting or one column, the site names.
Then created the simple feature object - (combined the geometry column with the df). Output shown:
A sf: 5 × 2
site pts.sfc
<fct> <MULTIPOINT [m]>
MANH MULTIPOINT ((-73.96448 40.7...
BKLN MULTIPOINT ((-73.96448 40.7...
BRON MULTIPOINT ((-73.96448 40.7...
QUEE MULTIPOINT ((-73.96448 40.7...
STAT MULTIPOINT ((-73.96448 40.7...
Then this sf object was written to the shape folder:
Writing layer `Weather_sites' to data source `C:/pathR' using driver `ESRI Shapefile'
Writing 5 features with 1 fields and geometry type Multi Point.
The output shows the 5 sites in NY and the states of NJ and NY. The NY shape file seems less inclusive and only includes some limited NY data.
Getting the entire mapping thing right involves doing the four steps in order. Geometry, geometry column, sf object, write to shape folder.
