-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy path1.5-ExtractGeoData-Water-PointsToVector.R
173 lines (124 loc) · 5.66 KB
/
1.5-ExtractGeoData-Water-PointsToVector.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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
# This script connects the point sources to the river and lakes network.
# The paths to the original data and the location for saving the processed data
# should be adapted.
#
# Author: Delphine Kawecki-Wenger
# Date of last modification: 02.10.2019
library(sf)
# load vector water data
load("Data_Modified/SurfaceWaters.Rdata")
# load point water data
load("Data_Modified/WWTP_CSO_sf.Rdata")
points_sf <- Data.sf
rm(Data.sf)
st_crs(points_sf) <- "+init=epsg:2056"
### HOW TO HANDLE POINT EMISSIONS
# for each emission point, create a vector of indices giving the nearest feature
# with lakes
ind.lakes <- st_nearest_feature(points_sf, lakes_emi)
# with rivers
ind.rivers <- st_nearest_feature(points_sf, rivers_emi)
# calculate minimal distance with lakes
dist.lakes <- rep(NA, length(ind.lakes))
for(point in 1:length(points_sf$geometry)){
dist.lakes[point] <- st_distance(points_sf[point,], lakes_emi[ind.lakes[point],])
}
# calculate minimal distance with rivers
dist.rivers <- rep(NA, length(ind.rivers))
for(point in 1:length(points_sf$geometry)){
dist.rivers[point] <- st_distance(points_sf[point,], rivers_emi[ind.rivers[point],])
}
# compare (= create a logical vector, TRUE if river is closer, FALSE if lake is closer)
log.river.closer <- (dist.rivers < dist.lakes)
# create a layer in the river or lake sf object corresponding to the points in the right order
add.layernames <- names(points_sf)[1:3]
for(lay in add.layernames){
lakes_emi[lay] <- rep(0,length(lakes_emi$geometry))
rivers_emi[lay] <- rep(0,length(rivers_emi$geometry))
}
### First: rivers
# create a nested empty list (first emission ayers; then river segments)
addpoint <- sapply(add.layernames, function(x) NULL)
for(i in 1:length(addpoint)){ addpoint[[i]] <- rep(0,length(rivers_emi$geometry)) }
# store point emission to layer (whole loop takes about 50 minutes)
message(paste(format(Sys.time(), "%H:%M:%S"),"started rivers"))
for(seg in 1:length(rivers_emi$geometry)){
# find point that needs to be added to the river segment
ind.point <- which(seg == ind.rivers)
# skip if river segment is never the closest to an emission point
if(length(ind.point) == 0){
next
} else if(length(ind.point) == 1){ # if only one point to add, easy
# check if the distance between point and lake is smaller, if yes, skip (it will be added to a lake instead)
if(!log.river.closer[ind.point]){ next }
# else, add the point to the river segment (except if NA)
for(lay in add.layernames){
if(!is.na(points_sf[lay][[1]][ind.point])){
addpoint[[lay]][ind.rivers[ind.point]] <- addpoint[[lay]][ind.rivers[ind.point]] + points_sf[lay][[1]][ind.point]
}
}
} else { # if more than one point to add per segment, loop over points and do the same
for(ip in ind.point){
# check if the distance between point and lake is smaller, if yes, skip (it will be added to a lake instead)
if(!log.river.closer[ip]){ next }
# else, add the point to the river segment (except if NA)
for(lay in add.layernames){
if(!is.na(points_sf[lay][[1]][ip])){
addpoint[[lay]][ind.rivers[ip]] <- addpoint[[lay]][ind.rivers[ip]] + points_sf[lay][[1]][ip]
}
}
}
}
}
message(paste(format(Sys.time(), "%H:%M:%S"),"finished rivers"))
# add stored emissions to final sf object
for(lay in add.layernames){
rivers_emi[lay] <- addpoint[[lay]]
}
### Second: lakes
# create a nested empty list (first lemission ayers; then lakes)
addpoint <- sapply(add.layernames, function(x) NULL)
for(i in 1:length(addpoint)){ addpoint[[i]] <- rep(0,length(lakes_emi$geometry)) }
# store point emission to layer (whole loop takes about 3-50 minutes)
message(paste(format(Sys.time(), "%H:%M:%S"),"started lakes"))
for(seg in 1:length(lakes_emi$geometry)){
# find point that needs to be added to the lake
ind.point <- which(seg == ind.lakes)
# skip if lake is never the closest to an emission point
if(length(ind.point) == 0){
next
} else if(length(ind.point) == 1){ # if only one point to add, easy
# check if the distance between point and lake is smaller, if not, skip
if(log.river.closer[ind.point]){ next }
# else, add the point to the lake (except if NA)
for(lay in add.layernames){
if(!is.na(points_sf[lay][[1]][ind.point])){
addpoint[[lay]][ind.lakes[ind.point]] <- addpoint[[lay]][ind.lakes[ind.point]] + points_sf[lay][[1]][ind.point]
}
}
} else { # if more than one point to add per segment, loop over points and do the same
for(ip in ind.point){
# check if the distance between point and lake is smaller, if yes, skip (it will be added to a lake instead)
if(log.river.closer[ip]){ next }
# else, add the point to the lake (except if NA)
for(lay in add.layernames){
if(!is.na(points_sf[lay][[1]][ip])){
addpoint[[lay]][ind.lakes[ip]] <- addpoint[[lay]][ind.lakes[ip]] + points_sf[lay][[1]][ip]
}
}
}
}
}
message(paste(format(Sys.time(), "%H:%M:%S"),"finished lakes"))
# add stored emissions to final sf object
for(lay in add.layernames){
lakes_emi[lay] <- addpoint[[lay]]
}
### check addition is correct
for(lay in add.layernames){
sumpoints <- sum(points_sf[lay][[1]], na.rm = T)
sumpolys <- sum(rivers_emi[lay][[1]], na.rm = T) + sum(lakes_emi[lay][[1]], na.rm = T)
stopifnot( sumpoints == sumpolys )
}
### save
save(lakes_emi, lakes_notCH, rivers_emi, file = "Data_Modified/SurfaceWaters_wPoints.Rdata")