-
Notifications
You must be signed in to change notification settings - Fork 0
/
select_offshore_points.R
124 lines (87 loc) · 4.52 KB
/
select_offshore_points.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
## ---------------------------------------------------------------------------------------------
data("sample_raster", package="basics")
df <- sample_raster$df
ras <- sample_raster$raster
lons <- sample_raster$lons
lats <- sample_raster$lats
## ---------------------------------------------------------------------------------------------
require(raster)
require(ggplot2)
## ---------------------------------------------------------------------------------------------
world <- rnaturalearth::ne_countries(scale = "small", returnclass = "sp")
world <- rgeos::gUnaryUnion(world)
## ---------------------------------------------------------------------------------------------
newcrs <- "+proj=wintri +lon_0=0 +lat_1=0 +x_0=0 +y_0=0 +datum=WGS84 +units=km +no_defs"
mworld <- sp::spTransform(world, newcrs)
mras <- projectRaster(ras, crs=newcrs, over=TRUE)
mpts <- sp::SpatialPoints(cbind(lons, lats), proj4string=sp::CRS("+proj=longlat"))
mpts <- sp::spTransform(mpts, newcrs)
## ---------------------------------------------------------------------------------------------
buff1 <- rgeos::gBuffer(mworld, width = 300, byid=TRUE)
e<-raster::erase(buff1, mworld)
plot(e)
## ---------------------------------------------------------------------------------------------
e300 <- spatialEco::remove.holes(spatialEco::remove.holes(e))
plot(e300)
## ---------------------------------------------------------------------------------------------
plot(e300)
plot(mworld, add=TRUE, col="grey")
## ---------------------------------------------------------------------------------------------
el <- as(mworld, "SpatialLines")
el <- raster::crop(el, mpts@bbox)
set.seed(123)
pts <- sp::spsample(el, n = 1, type = "regular")
## ---------------------------------------------------------------------------------------------
plot(e300, border="red", axes=TRUE, xlim=mpts@bbox[1,], ylim=mpts@bbox[2,])
plot(mworld, add=TRUE, col="grey")
plot(pts, add=TRUE, pch=19)
## ---------------------------------------------------------------------------------------------
el300 <- as(e300, "SpatialLines")
close_pt <- maptools::snapPointsToLines(pts, el300, maxDist=500)
plot(mworld)
plot(close_pt, add=TRUE, pch=19, col="red")
plot(pts, add=TRUE, pch=19)
## ---------------------------------------------------------------------------------------------
el300 <- as(e300, "SpatialLines")
df <- c()
n <- length(el300@lines[[1]]@Lines)
for (i in 1:n){
df <- rbind(df, cbind(el300@lines[[1]]@Lines[[i]]@coords, ID=i))
}
close_pt <- maptools::nearestPointOnLine(df, pts@coords)
close_pt <- sp::SpatialPoints(matrix(close_pt, ncol=2), proj4string=CRS(newcrs))
## ---------------------------------------------------------------------------------------------
plot(e300, border="red", axes=TRUE, xlim=mpts@bbox[1,], ylim=mpts@bbox[2,])
plot(mworld, add=TRUE, col="grey")
plot(pts, add=TRUE, pch=19)
plot(close_pt, add=TRUE, pch=1)
## ---------------------------------------------------------------------------------------------
plot(mras, axes=TRUE, xlim=bbox(mras)[1,], ylim=bbox(mras)[2,])
plot(e300, border="red", add=TRUE)
plot(mworld, add=TRUE, col="grey")
plot(pts, add=TRUE, pch=19)
plot(close_pt, add=TRUE, pch=1)
## ---------------------------------------------------------------------------------------------
raster::extract(mras, close_pt)
## ---------------------------------------------------------------------------------------------
circle_pt <- raster::buffer(close_pt, width = 100)
## ----echo=FALSE-------------------------------------------------------------------------------
plot(mras, axes=TRUE, xlim=bbox(mras)[1,], ylim=bbox(mras)[2,])
plot(e300, border="red", add=TRUE)
plot(mworld, add=TRUE, col="grey")
plot(pts, add=TRUE, pch=19)
plot(close_pt, add=TRUE, pch=19, col="red")
plot(circle_pt, add=TRUE)
## ---------------------------------------------------------------------------------------------
vals <- raster::extract(mras, circle_pt)[[1]]
vals
## ---------------------------------------------------------------------------------------------
mean(vals, na.rm=TRUE)
## ----eval=FALSE-------------------------------------------------------------------------------
## crs.wintri <- newcrs
## save(crs.wintri, file=file.path(here::here(), "data/crs_wintri.rda"))
## world.wintri <- mworld
## save(world.wintri, file=file.path(here::here(),"data/world_wintri.rda"))
## save(world, file=file.path(here::here(), "data/world.rda"))
## buffer300.wintri <- list(buf_line = el300, buf_poly = e300, buf_df = df, crs = newcrs)
## save(buffer300.wintri, file=file.path(here::here(), "data/buffer300_wintri.rda"))