-
Notifications
You must be signed in to change notification settings - Fork 0
/
agu_raster_LAvisualization.py
172 lines (147 loc) · 8.82 KB
/
agu_raster_LAvisualization.py
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
import main
import pandas as pd
import numpy as np
import funcs
import rasterio
# from rasterio import show
# Everything I need for this should be within the file "D:\Etienne\fall2022\agu_data"
## Data from CIMS
data = main.load_data()
bysite = main.average_bysite(data)
## Data from CRMS
perc = pd.read_csv(r"D:\Etienne\fall2022\agu_data\percentflooded.csv",
encoding="unicode escape")
perc['Simple site'] = [i[:8] for i in perc['Station_ID']]
perc = perc.groupby('Simple site').median()
wl = pd.read_csv(r"D:\Etienne\fall2022\agu_data\waterlevelrange.csv",
encoding="unicode escape")
wl['Simple site'] = [i[:8] for i in wl['Station_ID']]
wl = wl.groupby('Simple site').median()
marshElev = pd.read_csv(r"D:\Etienne\fall2022\CRMS_data\bayes2year\12009_Survey_Marsh_Elevation\12009_Survey_Marsh_Elevation.csv",
encoding="unicode escape").groupby('SiteId').median().drop('Unnamed: 4', axis=1)
SEC = pd.read_csv(r"D:\Etienne\fall2022\agu_data\12017_SurfaceElevation_ChangeRate\12017.csv",
encoding="unicode escape")
SEC['Simple site'] = [i[:8] for i in SEC['Station_ID']]
SEC = SEC.groupby('Simple site').median().drop('Unnamed: 4', axis=1)
acc = pd.read_csv(r"D:\Etienne\fall2022\agu_data\12172_SEA\Accretion__rate.csv", encoding="unicode_escape")[
['Site_ID', 'Acc_rate_fullterm (cm/y)']
].groupby('Site_ID').median()
## Data from Gee and Arc
jrc = pd.read_csv(r"D:\Etienne\summer2022_CRMS\run_experiments\CRMS_GEE_JRCCOPY2.csv", encoding="unicode_escape")[
['Simple_sit', 'Land_Lost_m2']
].set_index('Simple_sit')
gee = pd.read_csv(r"D:\Etienne\fall2022\agu_data\CRMS_GEE60pfrom2007to2022.csv",
encoding="unicode escape")[['Simple_sit', 'NDVI', 'tss_med', 'windspeed']]\
.groupby('Simple_sit').median().fillna(0) # filling nans with zeros cuz all nans are in tss because some sites are not near water
distRiver = pd.read_csv(r"D:\Etienne\fall2022\CRMS_data\totalDataAndRivers.csv",
encoding="unicode escape")[['Field1', 'distance_to_river_m', 'width_mean']].groupby('Field1').median()
nearWater = pd.read_csv(r"D:\Etienne\fall2022\agu_data\ALLDATA2.csv", encoding="unicode_escape")[
['Simple site', 'Distance_to_Water_m']
].set_index('Simple site')
floodfreq = pd.read_csv(r"D:\\Etienne\\fall2022\\agu_data\\floodFrequencySitePerYear.csv", encoding="unicode_escape")[[
'Simple site', 'Flood Freq (Floods/yr)'
]].set_index('Simple site')
# Concatenate
df = pd.concat([bysite, distRiver, nearWater, gee, jrc, marshElev, wl, perc, SEC, acc, floodfreq], axis=1, join='outer')
# Now clean the columns
# First delete columns that are more than 1/2 nans
tdf = df.dropna(thresh=df.shape[0]*0.5, how='all', axis=1)
# Drop uninformative features
udf = tdf.drop([
'Year (yyyy)', 'Accretion Measurement 1 (mm)', 'Year',
'Accretion Measurement 2 (mm)', 'Accretion Measurement 3 (mm)',
'Accretion Measurement 4 (mm)','Basins',
'Month (mm)', 'Average Accretion (mm)', 'Delta time (days)', 'Wet Volume (cm3)',
'Delta Time (decimal_years)', 'Wet Soil pH (pH units)', 'Dry Soil pH (pH units)', 'Dry Volume (cm3)',
'Measurement Depth (ft)', 'Plot Size (m2)', '% Cover Shrub', '% Cover Carpet', 'Direction (Collar Number)',
'Direction (Compass Degrees)', 'Pin Number', 'Observed Pin Height (mm)', 'Verified Pin Height (mm)',
'calendar_year', 'percent_waterlevel_complete',
'Average Height Shrub (cm)', 'Average Height Carpet (cm)' # I remove these because most values are nan and these vars are unimportant really
], axis=1)
# Address the vertical measurement for mass calculation (multiple potential outcome problem)
vertical = 'Accretion Rate (mm/yr)'
if vertical == 'Accretion Rate (mm/yr)':
udf = udf.drop('Acc_rate_fullterm (cm/y)', axis=1)
# Make sure multiplier of mass acc is in the right units
# udf['Average_Ac_cm_yr'] = udf['Accretion Rate (mm/yr)'] / 10 # mm to cm conversion
# Make sure subsidence and RSLR are in correct units
udf['Shallow Subsidence Rate (mm/yr)'] = udf[vertical] - udf['Surface Elevation Change Rate (cm/y)'] * 10
udf['Shallow Subsidence Rate (mm/yr)'] = [0 if val < 0 else val for val in udf['Shallow Subsidence Rate (mm/yr)']]
udf['SEC Rate (mm/yr)'] = udf['Surface Elevation Change Rate (cm/y)'] * 10
# Now calcualte subsidence and RSLR
# Make the subsidence and rslr variables: using the
udf['SLR (mm/yr)'] = 2.0 # from jankowski
udf['Deep Subsidence Rate (mm/yr)'] = ((3.7147 * udf['Latitude']) - 114.26) * -1
udf['RSLR (mm/yr)'] = udf['Shallow Subsidence Rate (mm/yr)'] + udf['Deep Subsidence Rate (mm/yr)'] + udf[
'SLR (mm/yr)']
udf = udf.drop(['SLR (mm/yr)', 'Latitude'],
axis=1) # obviously drop because it is the same everywhere ; only used for calc
elif vertical == 'Acc_rate_fullterm (cm/y)':
udf = udf.drop('Accretion Rate (mm/yr)', axis=1)
# Make sure multiplier of mass acc is in the right units
# udf['Average_Ac_cm_yr'] = udf[vertical]
# Make sure subsidence and RSLR are in correct units
udf['Shallow Subsidence Rate (mm/yr)'] = (udf[vertical] - udf['Surface Elevation Change Rate (cm/y)'])*10
udf['SEC Rate (cm/yr)'] = udf['Surface Elevation Change Rate (cm/y)']
# Now calcualte subsidence and RSLR
# Make the subsidence and rslr variables: using the
udf['SLR (mm/yr)'] = 2.0 # from jankowski
udf['Deep Subsidence Rate (mm/yr)'] = ((3.7147 * udf['Latitude']) - 114.26) * -1
udf['RSLR (mm/yr)'] = udf['Shallow Subsidence Rate (mm/yr)'] + udf['Deep Subsidence Rate (mm/yr)'] + udf[
'SLR (mm/yr)']*0.1
udf = udf.drop(['SLR (mm/yr)'],
axis=1) # obviously drop because it is the same everywhere ; only used for calc
else:
print("NOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO")
####### Define outcome as vertical component
outcome = vertical
# Try to semi-standardize variables
des = udf.describe() # just to identify which variables are way of the scale
udf['distance_to_river_km'] = udf['distance_to_river_m']/1000 # convert to km
udf['river_width_mean_km'] = udf['width_mean']/1000
udf['distance_to_water_km'] = udf['Distance_to_Water_m']/1000
udf['land_lost_km2'] = udf['Land_Lost_m2']*0.000001 # convert to km2
# Drop remade variables
udf = udf.drop(['distance_to_river_m', 'width_mean', 'Distance_to_Water_m', 'Soil Specific Conductance (uS/cm)',
'Soil Porewater Specific Conductance (uS/cm)',
'Land_Lost_m2'], axis=1)
udf = udf.rename(columns={'tss_med': 'tss med mg/l'})
# conduct outlier removal which drops all nans
rdf = funcs.outlierrm(udf.drop('Community', axis=1), thres=2.6)
# transformations (basically log transforamtions) --> the log actually kinda regularizes too
rdf['log_distance_to_water_km'] = [np.log10(val) if val > 0 else 0 for val in rdf['distance_to_water_km']]
rdf['log_river_width_mean_km'] = [np.log10(val) if val > 0 else 0 for val in rdf['river_width_mean_km']]
rdf['log_distance_to_river_km'] = [np.log10(val) if val > 0 else 0 for val in rdf['distance_to_river_km']]
# drop the old features
rdf = rdf.drop(['distance_to_water_km', 'distance_to_river_km', 'river_width_mean_km'], axis=1)
# Now it is feature selection time
# drop any variables related to the outcome
rdf = rdf.drop([ # IM BEING RISKY AND KEEP SHALLOW SUBSIDENCE RATE
'Surface Elevation Change Rate (cm/y)', 'Deep Subsidence Rate (mm/yr)', 'RSLR (mm/yr)', 'SEC Rate (mm/yr)',
'Shallow Subsidence Rate (mm/yr)', # potentially encoding info about accretion
# taking out water level features because they are not super informative
# Putting Human in the loop
'90th%Upper_water_level (ft NAVD88)', '10%thLower_water_level (ft NAVD88)', 'avg_water_level (ft NAVD88)',
'std_deviation_water_level(ft NAVD88)', 'Staff Gauge (ft)', 'Soil Salinity (ppt)',
'log_river_width_mean_km', # i just dont like this variable because it has a sucky distribution
'Soil Porewater Temperature (°C)',
'Bulk Density (g/cm3)', 'Organic Density (g/cm3)',
'Soil Moisture Content (%)', 'Organic Matter (%)',
], axis=1)
# Rename some variables for better text wrapping
rdf = rdf.rename(columns={
'Tide_Amp (ft)': 'Tide Amp (ft)',
'avg_percentflooded (%)': ' avg percent flooded (%)',
'Average_Marsh_Elevation (ft. NAVD88)': 'Average Marsh Elevation (ft. NAVD88)',
'log_distance_to_water_km': 'log distance to water km',
'log_distance_to_river_km': 'log distance to river km',
'10%thLower_flooding (ft)': '10%thLower flooding (ft)',
'90%thUpper_flooding (ft)': '90%thUpper flooding (ft)',
'avg_flooding (ft)': 'avg flooding (ft)',
'std_deviation_avg_flooding (ft)': 'std dev avg flooding (ft)'
})
rdf.to_csv(r"D:\Etienne\fall2022\agu_data\results\AGU_rdf_forEDA_dataset.csv")
# ### Raster visulaization
# fp = r"D:\\Etienne\\fall2022\\agu_data\\LAraster_Output.shp"
# raster = rasterio.open(fp)
# rasterio.show(raster)