-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathswe_inversion_feb19-feb26.R
164 lines (124 loc) · 4.71 KB
/
swe_inversion_feb19-feb26.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
# SWE inversion for 2/19-2/26
# jack tarricone
# feb 23, 2021
# update march 7th using pit location!
# systematically extract values from 8 points around cell pit
library(terra)
library(ggplot2)
# set home folder
setwd("/Users/jacktarricone/ch1_jemez_data/gpr_rasters_ryan/")
list.files() #pwd
# import corrected unwrapped phase data
unw_raw <-rast("unw_feb19-feb26.tif")
plot(unw_raw)
# import i_angle raster and resample to unw grid bc of slight extent difference
lidar_inc_raw <-rast("lidar_inc_deg.tif")
lidar_inc_v1 <-resample(lidar_inc_raw, unw_raw)
# set crop extent and crop for better visualization
crop_ext <-ext(-106.57, -106.38, 35.81, 35.96)
lidar_inc <-crop(lidar_inc_v1, crop_ext)
# mask and crop unwrapped phase data down to extent on new inc angle raster
unw_crop <-crop(unw_raw, crop_ext)
unw <-mask(unw_crop, lidar_inc)
# plot results
plot(lidar_inc)
plot(unw)
####################################
###### bring in fsca layers ########
####################################
# fsca
fsca_raw <-rast("landsat_fsca_2-18.tif")
# crop to .inc angle extent and mask
fsca_crop <-crop(fsca_raw, ext(lidar_inc))
fsca <-mask(fsca_crop, lidar_inc)
plot(fsca)
# create snow mask
snow_mask <-fsca
values(snow_mask)[values(snow_mask) > 1] = 1
plot(snow_mask)
# writeRaster(snow_mask,"02_18_2020_snow_mask.tif")
# masked the unwrapped phase with snow mask
unw_snow_mask <- mask(unw, snow_mask, maskvalue = NA)
plot(unw)
plot(unw_snow_mask)
########################################################
######### converting phase change to SWE ##############
########################################################
#######################################################################
# don't pick denisty and di_elec value
# pick a density and LWC (from ryans equations and field measurements)
# vary density and LWC over range over measured values
########################################################################
# table ryan sent over
pit_info <-read.csv("/Users/jacktarricone/ch1_jemez_data/pit_data/perm_pits.csv")
head(pit_info)
# ## define static information from pits
# # calculate density
# mean_density_feb12 <- pit_info$mean_density[1]
# mean_density_feb20 <- pit_info$mean_density[2]
#
# # mean density between two flights
# mean_density_feb12_20 <-(mean_density_feb12 + mean_density_feb20)/2
# dielctric constant k
k_feb20 <- pit_info$mean_k[2]
k_feb26 <- pit_info$mean_k[3]
# mean k between two flights
mean_k_feb20_26 <-(k_feb20+k_feb26)/2
#######################
#### swe inversion ####
#######################
# first step, define function for insar constant
# inc = incidence angle raster [deg]
# wL = sensor save length [cm]
# k = dielectric permittivty
# radar wave length from uavsar annotation file
uavsar_wL <- 23.8403545
insar_constant <-function(inc, wL, k){
(-(4*pi)/wL)*(cos(inc) - sqrt(k - sin((inc)^2)))
}
# create the insar constant raster
# because of the variation of .inc angle, this value you will change on pixelwise basis
insar_constant_rast <-insar_constant(inc = lidar_inc, wL = uavsar_wL, k = mean_k_feb20_26)
hist(insar_constant_rast, breaks = 100)
plot(insar_constant_rast)
#do swe change calc with masked unwrapped phase data
dswe_raw <-insar_constant_rast*unw_snow_mask
plot(dswe_raw)
hist(dswe_raw, breaks = 100)
# writeRaster(dswe_raw,"./final_swe_change/raw_dswe_feb20_26.tif")
#######################################
### calculating absolute SWE change ###
#######################################
# using swe change from the pit as "known" change point
# extent around gpr transect
gpr <-ext(-106.5255, -106.521, 35.856, 35.8594)
dswe_crop <-crop(dswe_raw, gpr)
plot(dswe_crop)
# pull out location info into separate df
loc <-data.frame(lat = pit_info$lat[1],
lon = pit_info$lon[1])
# plot pit location using terra vector funcitonality
pit_point <-vect(loc, geom = c("lon","lat"), crs = crs(unw))
points(pit_point, cex = 1)
# extract cell number from pit lat/lon point
pit_cell_v1 <-cells(dswe_raw, pit_point)
cell_number <-pit_cell_v1[1,2]
cell_number
# define neighboring cells by number and create a vector
neighbor_cells <-c(adjacent(dswe_raw, cells = cell_number, directions ="8"))
# add orginal cell back to vector
cell_vector <-c(cell_number, neighbor_cells)
# extract using that vector
nine_cell_dswe <-terra::extract(dswe_raw, cell_vector, cells = TRUE, xy = TRUE)
nine_cell_dswe
# mean of 9 swe changes around
mean_pit_dswe <-mean(nine_cell_dswe[1:9,1])
mean_pit_dswe
# subtract average swe change value to great absolute swe change
# therefor pixels around the pit will show no swe change
# which is consistent with what was observed on the ground
dswe_abs <-dswe_raw - mean_pit_dswe
plot(dswe_abs)
hist(dswe_abs, breaks = 100)
# save
# writeRaster(dswe_abs,"./final_swe_change/dswe_feb20-26.tif")