-
Notifications
You must be signed in to change notification settings - Fork 0
/
3.0_make_NatCap_files.r
153 lines (118 loc) · 7.78 KB
/
3.0_make_NatCap_files.r
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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
# Make a shapefile of our Areas of Interest, for use in the NatCap's InVest software
# Also make a csv of lat/lon for landing points and grid connection points
#######################
## Load libraries
#######################
#require(maptools)
require(RColorBrewer)
#require(rgdal)
#require(PBSmapping)
#require(rgeos) # for gBuffer
require(sf)
##############################################################
## Read in and process data to create Areas of Interest (AOI)
##############################################################
climgrid <- readRDS('temp/SPsf2.rds') # the projection grid
# label by area of interest regions
climgrid$AOI <- NA
climgrid$AOI[climgrid$lon < -100 | climgrid$lon > 0] <- 'west'
climgrid$AOI[climgrid$lon >= -100 & climgrid$lon < 0] <- 'east'
sum(is.na(climgrid$AOI)) # 0
# plot(climgrid['AOI'], lwd=0.001, axes=TRUE) # Aleutians are on the far right
# plot(climgrid['AOI'], lwd=0.5, xlim=c(-70, -69), ylim=c(41,42), axes=TRUE) # zoom in (Cape Cod)
# project to North American Albers, but in WGS84 rather than NAD83
climgridp <- st_transform(climgrid, crs='+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=WGS84 +units=m +no_defs') # North America Albers Equal Area Conic from https://epsg.io/102008, except with WGS84
# plot(climgridp['AOI'], lwd=0.001, axes=TRUE) # Aleutians are on the far right
# buffer slightly (1000 m) to aid merging the grid cells together
climgridbuf <- st_buffer(climgridp, 1000)
# merge all into multipart polygons based on AOI
climpoly <- aggregate(climgridbuf[,c('geometry')], by = list(AOI = climgridbuf$AOI), do_union = TRUE, FUN=function(x, ...) return(1))
# plot(climpoly['AOI'], lwd=0.01, axes=TRUE)
# plot(climpoly['AOI'], lwd=1, xlim=c(1.9e6, 2.2e6), ylim=c(4e5,6e5), axes=TRUE) # zoom in (Cape Cod)
# buffer whole analysis area
climpolyb <- st_buffer(climpoly, dist = 100000) # 100km
# plot(climpolyb, lwd=0.2, axes=TRUE)
# write out each region as a separate shapefile
# will throw an error if file already exists
for(i in 1:nrow(climpolyb)){
st_write(climpolyb[i,], paste0('temp/AOI_', climpolyb$AOI[i], '.shp')) # write the .prj file as part of the shapefile. Overwrite
}
#######################################
## Make land and grid connection points
#######################################
# set parameters
crs <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs' # North America Albers Equal Area Conic from https://epsg.io/102008
crslatlong = "+init=epsg:4326"
# read in files
coastlines <- st_read('data/natcap/NAmainland_lines.shp') # global coast layer from NatCap, trimmed to North America
# plot(coastlines['id'], axes=TRUE) # a bit slow
townsUS <- st_read('dataDL/usgs/citiesx020_nt00007/citiesx020.shp') # US National Map cites and towns
st_crs(townsUS) <- 4269 # set CRS to NAD83
# plot(townsUS['NAME'], pch=16, cex=0.2, axes=TRUE)
# plot(st_geometry(townsUS), add=TRUE, col='red') # to add to coastlines plot. Not working for some reason.
townsCan <- st_read('dataDL/statcan/lpc_000b16a_e/lpc_000b16a_e.shp') # Canadian census population centres
# project to North American Albers
coastlinesp <- st_transform(coastlines, crs = crs)
townsUSp <- st_transform(townsUS, crs = crs)
townsCanp <- st_transform(townsCan, crs = crs)
# find centroids for Canada
townsCanpc <- st_centroid(townsCanp)
plot(st_geometry(coastlinesp), lwd=0.2, axes=TRUE)
plot(st_geometry(townsUSp), pch=16, cex=0.1, add=TRUE, col='blue') # plots on top of NA well
plot(st_geometry(townsCanpc), pch=16, cex=0.1, add=TRUE, col='green') # plots on top of NA well
# trim US towns to those >1000 people (Canada is already trimmed)
townsUS1000 <- townsUSp[which(townsUSp$POP_2000 >= 1000),]
dim(townsUSp) # 35432
dim(townsUS1000) # 9992
# combine US and Canada
names(townsCanpc)[names(townsCanpc)=='PCNAME'] <- 'NAME'
towns1000 <- rbind(townsUS1000[,c('NAME', 'geometry')], townsCanpc[,c('NAME', 'geometry')])
# trim town to those <50km from the NA coast
nearcoast <- st_is_within_distance(towns1000, coastlinesp, dist = 50*1000, sparse = FALSE) # slow (a few min). units in meters (50km). returns a matrix with columns corresponding to each coastline (includes a few major islands)
nearcoastany <- rowSums(nearcoast) > 0 # sum across rows: we don't care which coast a town is close to
sum(nearcoast) # 2233
sum(nearcoastany) # 1873. shows that some towns were close to multiple coastlines
plot(st_geometry(towns1000[which(nearcoastany),]), col='red', pch=16, cex=1, add=TRUE) # adds to the plot before
towns1000nearcoast <- towns1000[which(nearcoastany),]
# find nearest point on the coastline to each town. So much faster than the old method of sampling!
#coastpoint.near <- st_as_sf(rgeos::gNearestPoints(as(towns1000nearcoast,"Spatial"), as(coastlinesp,"Spatial"))[2,]) # from https://gis.stackexchange.com/questions/288570/find-nearest-point-along-polyline-using-sf-package-in-r. But only returns one point. Would have to implement in a loop
towncoastlines <- st_nearest_points(towns1000nearcoast, coastlinesp) # returns LINESTRINGS from first to second geometry
towncoastlengths <- st_length(towncoastlines) # length of each line, so that we can find the closest coastline to each town
nearestptinds <- aggregate(list(ind = towncoastlengths), by = list(town = rep(1:nrow(towns1000nearcoast), each = nrow(coastlinesp))), FUN = function(x) which.min(x)) # find index of shortest line from town to a coast. Works because st_nearest_points returns a vector where y cycles fastest and x cycles slowest
nearestptinds2 <- nearestptinds$ind + seq(0, length.out = nrow(nearestptinds), by = nrow(coastlinesp)) # convert to an index into towncoastlines
pts <- st_cast(towncoastlines[nearestptinds2], "POINT") # gives all start (towns) & end (coastlines) points, alternating
coastpts <- pts[seq(2, length(pts), 2)] # just the end points (on coastlines)
length(coastpts)
#plot(st_geometry(towns1000nearcoast[1:10,]), axes = TRUE, col = 'blue') # to plot a few towns
plot(st_geometry(towns1000nearcoast), axes = TRUE, col = 'blue', xlim = c(-2.3e6, -1.8e6), ylim = c(-4e5, 0)) # to plot a few towns
#plot(st_geometry(towns1000nearcoast[1:100,]), axes = TRUE, col = 'blue') # to plot many towns
plot(st_geometry(coastlinesp), add = TRUE)
plot(st_geometry(coastpts), add = TRUE, col = 'red')
# project back to latlong in order to make a table for NatCap InVEST
coastpts.ll <- st_transform(coastpts, crs = crslatlong)
towns1000nearcoast.ll <- st_transform(towns1000nearcoast, crs = crslatlong)
# make output table of town and nearest landing point locations
# format as specified by NatCap InVEST Wave
towns.coords <- st_coordinates(towns1000nearcoast.ll)
coast.coords <- st_coordinates(coastpts.ll)
outgrid <- data.frame(ID=1:nrow(towns.coords), LAT=towns.coords[,2], LONG=towns.coords[,1], TYPE='GRID', LOCATION=towns1000nearcoast.ll$NAME)
outland <- data.frame(ID=(1+nrow(towns.coords)):(nrow(towns.coords)+nrow(coast.coords)), LAT=coast.coords[,2], LONG=coast.coords[,1], TYPE='LAND', LOCATION=towns1000nearcoast.ll$NAME)
# make sure it looks OK
plot(outland$LONG, outland$LAT, pch=16, cex=0.5)
points(outgrid$LONG, outgrid$LAT, col='red', pch=16, cex=0.5) # the towns
# combine town and landings points
out <- rbind(outgrid, outland)
head(out[out$TYPE == 'GRID', ])
head(out[out$TYPE == 'LAND', ])
tail(out[out$TYPE == 'GRID', ])
tail(out[out$TYPE == 'LAND', ])
# write out
write.csv(out, file='output/landgridpts_northamerica.csv', row.names=FALSE)
# make output table of town and nearest landing point locations
# format as specified by NatCap InVEST Wind
out2 <- out[, c('ID', 'TYPE', 'LAT', 'LONG')]
names(out2)[names(out2)=='LAT'] <- 'LATI'
head(out2)
tail(out2)
# write out
write.csv(out2, file='output/landgridpts_northamerica_wind.csv', row.names=FALSE)