-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathSpatial_Join
47 lines (35 loc) · 2.15 KB
/
Spatial_Join
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#This section of the script joins the names of states from a states polygon
#to a hurricane points shapefile.
#Import rgdal
library(rgdal)
#Import point Hurricane data. This dataset was clipped to only hurricane points that
#fall over the United States. Hurricane points that fell over the Atlantic Ocean were
#omitted.
HurricanesUS <- readOGR(dsn = "data/HurricanesUS.shp", layer = "HurricanesUS")
#Import states polygon shapefile
states <- readOGR(dsn = "data\\UnitedStates.shp", layer = "UnitedStates")
#Use the over() function to create a new. "STATE_NAME" is a field name in the states shapefile.
Polys2Pts <- over(HurricanesUS, states[,"STATE_NAME"])
#Write the new data to the polygon (or shapefile) that you want to have the join.
#First write the object shapefile that you want to write the new data to and create
#a new field by writing $new field. Then write the new object created from the over()
#function with the field name. This will create a new data frame with no spatial components.
#Do not filter or sort either tables before writing the next command.
HurricanesUS$states <- Polys2Pts$STATE_NAME
#Verify the spatial join worked by creating a new shapefile and open in GIS program.
writeOGR(obj = HurricanesUS, dsn = "tempdir", layer = "Hurricane_states", driver = "ESRI Shapefile")
############################################################################################
#This section of the script calculates the percentage and sum of hurricane points
#that fell over the United States.
#Import sp library
library(sp)
#Import large point Hurricane data. Data is within the US and Atlantic Ocean.
Hurricane <- readOGR(dsn = "data/Hurricane.shp", layer = "Hurricane")
#Get the percentage of the hurricanes that land in the United States. Use is.na()
#and over() function to calculate the total. Use the "SpatialPolygons" expression to
#write the output as a SpatialPolygons object instead of a SpatialPolygonsDataFrame.
inside.states <- !is.na(over(Hurricane, as(states, "SpatialPolygons")))
#View the mean of the number of hurricanes that fall in the US
mean(inside.states)
#View the total number of hurricanes that fall in the US
sum(inside.states)