-
Notifications
You must be signed in to change notification settings - Fork 0
/
greatlakes.R
166 lines (130 loc) · 5.78 KB
/
greatlakes.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
# Maps connecting the Great Lakes to the ocean
# www.overfitting.net
# https://www.overfitting.net/2024/01/conectando-los-grandes-lagos-al-mar-con.html
library(terra) # read GeoTIFF, reprojection, crop and resample
library(tiff) # save 16-bit TIFF's
library(png) # save 8-bit PNG's
hillshademap=function(DEM, dx=25, dlight=c(0, 2, 3), gamma=1) {
# hillshademap() inputs DEM data and outputs a hillshade matrix
#
# DEM: digital elevation map matrix
# dx: DEM resolution (cell size in the same units as elevation values)
# dlight: lighting direction (3D vector defined from observer to light source):
# dlight=c(0, 2, 3) # sunrise
# dlight=c(0, 0, 1) # midday
# dlight=c(0,-2, 3) # sunset
# gamma: optional output gamma lift
DIMY=nrow(DEM)
DIMX=ncol(DEM)
# If array turn its first dimension into a matrix
if (!is.matrix(DEM)) DEM=matrix(DEM[,,1], nrow=DIMY, ncol=DIMX)
dlightM=sum(dlight^2)^0.5
# Vectorial product to calculate n (orthogonal vector)
nx = 2*dx*(DEM[1:(DIMY-2), 2:(DIMX-1)] - DEM[3:DIMY, 2:(DIMX-1)])
ny = 2*dx*(DEM[2:(DIMY-1), 1:(DIMX-2)] - DEM[2:(DIMY-1), 3:DIMX])
nz = 4*dx^2
nM = (nx^2 + ny^2 + nz^2)^0.5
# Dot product to calculate cos(theta)
dn = dlight[1]*nx + dlight[2]*ny + dlight[3]*nz # (DIMY-2)x(DIMX-2) matrix
# Reflectance (=cos(theta))
hillshadepre=dn/(dlightM*nM)
hillshadepre[hillshadepre<0]=0 # clip negative values
# Add 1-pix 'lost' borders
hillshademap=matrix(0, nrow=DIMY, ncol=DIMX)
hillshademap[2:(DIMY-1), 2:(DIMX-1)]=hillshadepre
rm(hillshadepre)
hillshademap[c(1,DIMY),]=hillshademap[c(2,DIMY-1),]
hillshademap[,c(1,DIMX)]=hillshademap[,c(2,DIMX-1)]
return(hillshademap^(1/gamma))
}
###########################################################
# 1. PROCESS GEOTIFF DATA TO GET THE DEM AS A MATRIX
# https://download.gebco.net/
# The GEBCO_2023 Grid is a global terrain model for ocean and land,
# providing elevation data, in meters, on a 15 arc-second interval grid
# of 43200 rows x 86400 columns, giving 3,732,480,000 data points.
# The data values are pixel-centre registered i.e. they refer to elevations,
# in meters, at the centre of grid cells.
lakes=rast("geotiff_northamerica.tif")
lakes
plot(lakes)
# CROP Great Lakes (in long/lant degrees)
cropdef=ext(-100, -70, 38, 52)
lakes=crop(x=lakes, y=cropdef, threads=TRUE)
lakes
plot(lakes)
# REPROJECT raster from Longitude Latitude (+proj=longlat)/WGS84
# to Lambert Conic Conformal (+proj=lcc)/WGS84
# https://pygis.io/docs/d_understand_crs_codes.html
# https://stackoverflow.com/questions/36868506/how-to-change-a-lambert-conic-conformal-raster-projection-to-latlon-degree-r
CRS="+proj=lcc +ellps=WGS84 +lat_1=33 +lat_2=45 +lon_0=-84 +units=km"
lakes=project(x=lakes, y=CRS, method='bilinear', threads=TRUE)
lakes
plot(lakes)
# CROP raster to area of interest (in km)
cropdef=ext(-675, 650, 4800, 5800)
lakes=crop(x=lakes, y=cropdef, threads=TRUE)
lakes
plot(lakes)
# resolution : 0.3538445, 0.3538445 (x, y) in kms
# 0.3538445km grid resolution
resolution=res(lakes)[1]*1000 # 353.8445m grid resolution
# Convert to matrix and save as TIFF
DEM=matrix(as.array(lakes), nrow=nrow(lakes))
hist(DEM, breaks=1000)
abline(v=0, col='red')
DEM=DEM-min(DEM)
writeTIFF(DEM/max(DEM), "lakes.tif", compression='LZW', bits.per.sample=16)
DEM=matrix(as.array(lakes), nrow=nrow(lakes))
DIMY=nrow(DEM)
DIMX=ncol(DEM)
###########################################################
# 2. PROCESS MATRIX TO OBTAIN MAP CONTOURS AND HILLSHADE
# Superior, Michigan, Huron, Erie contours: ~179m above sea level
# Ontario contours: 74m above sea level
# Sea level contours: 0m above sea level
lakes=c('currentmainlakes', 'currentontariolake', 'simulated')
depths=c(179, 74, 0)
for (i in 1:length(lakes)) {
# Calculate solid map contour
solid=DEM
if (depths[i]<=0) {
solid[solid>=depths[i]]=1 # set areas to 1 (land)
solid[solid<depths[i]]=0 # set areas to 0 (water)
} else {
solid[solid<depths[i]]=0 # set areas to 0 (water)
solid[solid>=depths[i]]=1 # set areas to 1 (land)
}
writePNG(solid, paste0("mapsolid_", lakes[i], ".png"))
# Calculate outline map from solid map
outline=solid*0
# 1 pixel thickness outline
outline[2:(DIMY-1), 2:(DIMX-1)]=
abs(solid[1:(DIMY-2), 2:(DIMX-1)] -
solid[2:(DIMY-1), 2:(DIMX-1)]) +
abs(solid[2:(DIMY-1), 1:(DIMX-2)] -
solid[2:(DIMY-1), 2:(DIMX-1)])
# increase to 2 pixel thickness outline
outline[2:(DIMY-1), 2:(DIMX-1)]=outline[2:(DIMY-1), 2:(DIMX-1)]+
outline[1:(DIMY-2), 2:(DIMX-1)]+outline[2:(DIMY-1), 3:(DIMX-0)]
# increase to 3 pixel thickness outline
outline[2:(DIMY-1), 2:(DIMX-1)]=outline[2:(DIMY-1), 2:(DIMX-1)]+
outline[1:(DIMY-2), 2:(DIMX-1)]+outline[3:(DIMY-0), 2:(DIMX-1)]+
outline[2:(DIMY-1), 1:(DIMX-2)]+outline[2:(DIMY-1), 3:(DIMX-0)]
outline[outline!=0]=1
writePNG(outline, paste0("mapoutline_", lakes[i], ".png"))
}
# Calculate grayscale hillshade
MIX=0.85 # two light sources are mixed to fill dark areas a bit
hill=hillshademap(DEM, dx=resolution, dlight=c(0, 2, 3))
hillfill=hillshademap(DEM, dx=resolution, dlight=c(0, 3, 2))
hill=hill*MIX+hillfill*(1-MIX)
gamma=1/2.0
hill=(hill/max(hill))^(1/gamma) # darken hillshade a bit
# Save hillshade
writeTIFF(hill, "hillshade.tif",
bits.per.sample=16, compression="LZW")
# Display hillshade
image(t(hill[nrow(hill):1,]), useRaster=TRUE,
col=c(gray.colors(256, start=0, end=1, gamma=0.5)),
asp=nrow(hill)/ncol(hill), axes=FALSE)