Make a Survey Overview¶
Goal: Make a nice of all the data on a nice DEM
Approach:
Create an overview DEM with hillShade
Grab locations of single location type data
Grab centroids of all the raster tiles
Plot it all!
Process:¶
Step 1. Create an Overview Raster with HillShade¶
[5]:
from snowexsql.db import get_db
from snowexsql.data import ImageData, LayerData, PointData, SiteData
from snowexsql.conversions import raster_to_rasterio
from rasterio.plot import show
from sqlalchemy.sql import func
from geoalchemy2.types import Raster
import geoalchemy2.functions as gfunc
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import Polygon
from geoalchemy2.shape import from_shape, to_shape
# Connect to the database we made.
db_name = 'db.snowexdata.org/snowex'
engine, session = get_db(db_name, credentials='./credentials.json')
# DEM data name and surveyor
data_name = 'DEM'
observers = 'USGS'
# Resolution to make our DEM/Hillshade
res = 20
# Form a polygon to clip at the max extent (figured out in advance)
x1 = 735200.0
x2 = 760000.0
y1 = 4319989.0
y2 = 4329803.0
extent = Polygon([[x1, y1], [x1, y2], [x2, y2], [x2, y1]]) # Construct it using shapely
# Make the polygon usable to the db
extent_ewkt = from_shape(extent, srid=26912).ST_AsEWKT()
# Make polygon plottable for later
extent_df = gpd.GeoDataFrame({'geometry':[extent]})
[6]:
# Define a function to reduce the same code used for getting the dem and hillshade
def filter_and_return(session, base, data_name, observers, extent):
'''
Small function to apply redundent filters and raster making
'''
# Save our query as a Tiff and clip it along the extents polygon
q = session.query(func.ST_AsTiff(func.ST_Clip(base, extent)))
# Filter by data name and surveyor
q = q.filter(ImageData.type == data_name)
q = q.filter(ImageData.observers == observers)
# Execute the query
rasters = q.all()
# Convert the dataset from the DB to rasterio
dataset = raster_to_rasterio(session, rasters)[0]
return dataset
[7]:
# Create the base query to reduce code: Collect Rasters and rescale them to our resolution and use bilinear interpolation
bq = func.ST_Rescale(ImageData.raster, res, -1 * res, 'blinear')
# Retrieve the dem, join all the tiles retrieved
base = gfunc.ST_Union(bq, type_=Raster) # Is NOT EXECUTED until query is executed.
dem = filter_and_return(session, base, data_name, observers, extent_ewkt)
# Merge all the tiles retrieved, then make the hillshade.
base = func.ST_Hillshade(gfunc.ST_Union(bq, type_=Raster)) # Is NOT EXECUTED until query is executed.
hillshade = filter_and_return(session, base, data_name, observers, extent_ewkt)
---------------------------------------------------------------------------
UndefinedColumn Traceback (most recent call last)
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/base.py in _execute_context(self, dialect, constructor, statement, parameters, execution_options, *args, **kw)
1705 if not evt_handled:
-> 1706 self.dialect.do_execute(
1707 cursor, statement, parameters, context
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/default.py in do_execute(self, cursor, statement, parameters, context)
716 def do_execute(self, cursor, statement, parameters, context=None):
--> 717 cursor.execute(statement, parameters)
718
UndefinedColumn: column images.observers does not exist
LINE 3: WHERE public.images.type = 'DEM' AND public.images.observers...
^
The above exception was the direct cause of the following exception:
ProgrammingError Traceback (most recent call last)
<ipython-input-7-8ae6983f1c2d> in <module>
4 # Retrieve the dem, join all the tiles retrieved
5 base = gfunc.ST_Union(bq, type_=Raster) # Is NOT EXECUTED until query is executed.
----> 6 dem = filter_and_return(session, base, data_name, observers, extent_ewkt)
7
8 # Merge all the tiles retrieved, then make the hillshade.
<ipython-input-6-24f4894732c7> in filter_and_return(session, base, data_name, observers, extent)
12
13 # Execute the query
---> 14 rasters = q.all()
15
16 # Convert the dataset from the DB to rasterio
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/orm/query.py in all(self)
2697 :ref:`faq_query_deduplicating`
2698 """
-> 2699 return self._iter().all()
2700
2701 @_generative
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/orm/query.py in _iter(self)
2832
2833 statement = self._statement_20()
-> 2834 result = self.session.execute(
2835 statement,
2836 params,
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/orm/session.py in execute(self, statement, params, execution_options, bind_arguments, _parent_execute_state, _add_event, **kw)
1673
1674 conn = self._connection_for_bind(bind, close_with_result=True)
-> 1675 result = conn._execute_20(statement, params or {}, execution_options)
1676
1677 if compile_state_cls:
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/base.py in _execute_20(self, statement, parameters, execution_options)
1519 )
1520 else:
-> 1521 return meth(self, args_10style, kwargs_10style, execution_options)
1522
1523 def exec_driver_sql(
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/sql/elements.py in _execute_on_connection(self, connection, multiparams, params, execution_options, _force)
311 ):
312 if _force or self.supports_execution:
--> 313 return connection._execute_clauseelement(
314 self, multiparams, params, execution_options
315 )
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/base.py in _execute_clauseelement(self, elem, multiparams, params, execution_options)
1388 linting=self.dialect.compiler_linting | compiler.WARN_LINTING,
1389 )
-> 1390 ret = self._execute_context(
1391 dialect,
1392 dialect.execution_ctx_cls._init_compiled,
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/base.py in _execute_context(self, dialect, constructor, statement, parameters, execution_options, *args, **kw)
1747
1748 except BaseException as e:
-> 1749 self._handle_dbapi_exception(
1750 e, statement, parameters, cursor, context
1751 )
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/base.py in _handle_dbapi_exception(self, e, statement, parameters, cursor, context)
1928 util.raise_(newraise, with_traceback=exc_info[2], from_=e)
1929 elif should_wrap:
-> 1930 util.raise_(
1931 sqlalchemy_exception, with_traceback=exc_info[2], from_=e
1932 )
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/util/compat.py in raise_(***failed resolving arguments***)
209
210 try:
--> 211 raise exception
212 finally:
213 # credit to
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/base.py in _execute_context(self, dialect, constructor, statement, parameters, execution_options, *args, **kw)
1704 break
1705 if not evt_handled:
-> 1706 self.dialect.do_execute(
1707 cursor, statement, parameters, context
1708 )
~/projects/venv/snowexenv/lib/python3.8/site-packages/sqlalchemy/engine/default.py in do_execute(self, cursor, statement, parameters, context)
715
716 def do_execute(self, cursor, statement, parameters, context=None):
--> 717 cursor.execute(statement, parameters)
718
719 def do_execute_no_params(self, cursor, statement, context=None):
ProgrammingError: (psycopg2.errors.UndefinedColumn) column images.observers does not exist
LINE 3: WHERE public.images.type = 'DEM' AND public.images.observers...
^
[SQL: SELECT ST_AsTIFF(ST_Clip(ST_Union(ST_Rescale(public.images.raster, %(ST_Rescale_1)s, %(ST_Rescale_2)s, %(ST_Rescale_3)s)), ST_AsEWKT(ST_GeomFromWKB(%(ST_GeomFromWKB_1)s, %(ST_GeomFromWKB_2)s)))) AS "ST_AsTIFF_1"
FROM public.images
WHERE public.images.type = %(type_1)s AND public.images.observers = %(observers_1)s]
[parameters: {'ST_Rescale_1': 20, 'ST_Rescale_2': -20, 'ST_Rescale_3': 'blinear', 'ST_GeomFromWKB_1': <memory at 0x7ff354c411c0>, 'ST_GeomFromWKB_2': 26912, 'type_1': 'DEM', 'observers_1': 'USGS'}]
(Background on this error at: http://sqlalche.me/e/14/f405)
Step 2. Grab locations of single location type data¶
[4]:
# Write a function to retrieve positional points of data
def grab_points(session, TableClass, data_name, distinct=False, instrument=None, downsample=None):
"""
Returns a dataframe of geometric points of requested data. Use distinct to avoid collecting
identical locations from layer data. Use instrument to isolate snow depths
Returns: df: Geopandas dataframe with one column geometry for plotting
"""
q = session.query(TableClass.geom).filter(TableClass.type == data_name)
if instrument != None:
# Filter by what instruments are in a list provided
if type(instrument) == list:
q = q.filter(TableClass.instrument.in_(instrument))
else:
q = q.filter(TableClass.instrument == instrument)
if downsample is not None:
# Only sample some of the data
q = q.filter(TableClass.id % downsample == 0)
if distinct:
q = q.distinct()
points = q.all()
# Make the points into a geopandas dataframe so we can easily plot it
df = gpd.GeoDataFrame({'geometry':[to_shape(p[0]) for p in points]})
print('Found {} points for {}...'.format(len(df.index), data_name))
return df
[5]:
# Show all instruments used to gather snow depths
tools = session.query(PointData.instrument).filter(PointData.type == 'depth').distinct().all()
# Grab our pits by selecting hand hardness profiles
data = {}
# Use distinct locations from hand hardness profiles to get pit locations
pits = grab_points(session, LayerData, 'hand_hardness', distinct=True)
# Get distinct locations of smp profiles
smp = grab_points(session, LayerData, 'force', distinct=True)
# Grab all snow depths collected with magnaprobe, mesa, or a pit rule
depths = grab_points(session, PointData, 'depth', instrument=['magnaprobe', 'mesa'])
# Grab all the GPR point data
gpr = grab_points(session, PointData, 'two_way_travel', downsample=100)
# Grab all the camera locations point data
cameras = grab_points(session, PointData, 'depth', instrument=['camera'], distinct=True)
Found 155 points for hand_hardness...
Found 931 points for force...
Found 37754 points for depth...
Found 12649 points for two_way_travel...
Found 29 points for depth...
Step 3. Grab centroids of all the raster tiles¶
[6]:
# Define a function to grab the center of each raster tile
def get_tile_centers(session, data_name, observers=None):
'''
Simple function to grab the center of each tile given a data type and optionally a surveyor name
'''
# Use database to grab the centroid of each tile outline (envelope) filtering on type. Also return the surveyor.
q = session.query(func.ST_Centroid(func.ST_Envelope(ImageData.raster))).filter(ImageData.type == data_name)
# If observers is provided, filter on that too
if observers != None:
q = q.filter(ImageData.observers == observers)
centers = q.all()
# Form the data into plottable geopandas dataframe with only a geometry column
df = gpd.GeoDataFrame({'geometry':[to_shape(p[0]) for p in centers]})
print('Found {} tiles for {}...'.format(len(df.index), data_name))
return df
[7]:
# Grab all the names of the raster data so we know what to reference
names = session.query(ImageData.type).distinct().all()
# Use the get_tile_centers function to grab the dataframe containing the tile centroids and assign a color and marker to it for plotting
tiles = {}
# Grab all the ASO DEM centers, assign color as blue and use squares for symbols
tiles['ASO Depths'] = (get_tile_centers(session, 'depth', observers='ASO Inc.'), 'steelblue', 's')
# tiles['ASO SWE'] = (get_tile_centers(session, 'swe', observers='ASO Inc.'), 'plum', 's')
# Grab all the USGS DEM centers, assign color as light red and use pentagons for symbols
tiles['USGS DEM'] = (get_tile_centers(session, 'DEM', observers='USGS'), 'plum','p')
# Grab all the insar data centers, assign color as some shade of orange and use diamonds for symbols
tiles['INSAR Amplitudes'] = (get_tile_centers(session, 'insar amplitude'), 'gold', 'D')
tiles['INSAR Derived'] = (get_tile_centers(session, 'insar interferogram real'), 'goldenrod', 'D') # Since imaginary and real are in the same location we only need one of them
# tiles['INSAR Correlation.'] = (get_tile_centers(session, 'insar correlation'), 'bisque', 'D')
Found 256 tiles for depth...
Found 1713 tiles for DEM...
Found 1344 tiles for insar amplitude...
Found 672 tiles for insar interferogram real...
[8]:
metadata = {'instrument': [], 'date':[],'doi':[]}
for k in metadata.keys():
for tbl in [PointData, LayerData, ImageData]:
print(f'Counting {k} in {tbl.__name__}...')
metadata[k] += [v[0] for v in session.query(getattr(tbl,k)).filter(getattr(tbl,k).isnot(None)).distinct()]
metadata[k] = list(set(metadata[k]))
# Find all observers (Temporarily not working because surveyor upload is incorrect)
metadata['observers'] = []
for tbl in [PointData, LayerData, ImageData]:
print(f'Counting observers in {tbl.__name__}...')
# Manage multiple names in the surveyor
observers = session.query(tbl.observers).filter(tbl.observers.isnot(None)).distinct().all()
for name in observers:
if ',' in name[0] and 'JPL' not in name[0]:
metadata['observers'] += name[0].split(',')
elif '&' in name[0]:
metadata['observers'] += name[0].split('&')
else:
metadata['observers'] += [name[0]]
metadata['observers'] = len(list(set(metadata['observers'])))
metadata['date'] = (min(metadata['date']).strftime('%Y-%m-%d'), max(metadata['date']).strftime('%Y-%m-%d'))
metadata['instrument'] = len(metadata['instrument'])
metadata['doi'] = len(metadata['doi'])
Counting instrument in PointData...
Counting instrument in LayerData...
Counting instrument in ImageData...
Counting date in PointData...
Counting date in LayerData...
Counting date in ImageData...
Counting doi in PointData...
Counting doi in LayerData...
Counting doi in ImageData...
Counting surveyors in PointData...
Counting surveyors in LayerData...
Counting surveyors in ImageData...
Step 4. Plot it all!¶
[9]:
# Create a figure with one subplot. Give it size
fig, ax = plt.subplots(1,1, figsize=(32, 16))
# Plot the hillshade raster in gray
show(hillshade, ax=ax, cmap='gray', transform=hillshade.transform)
# Plot the DEM with partial transparency so the hillshade can be seen
show(dem, ax=ax, alpha=0.5, cmap='terrain', transform=dem.transform)
# Plot raster centers
for n, d in tiles.items():
df = d[0]
color = d[1]
marker = d[2]
# Filter the data only to our map extent defined in the beginning
ind = df.geometry.within(extent)
# Plot with a name that shows the total tiles not just the tiles in the extent.
df.loc[ind].plot(ax=ax, marker=marker, markersize=30, edgecolor='black', alpha=0.5, color=color, label='{} tiles of {}'.format(len(df.index), n))
# Plot the GRP points as orange pixels
gpr.plot(ax=ax, marker=',' , color='orange', markersize=1, label='{} GPR Points'.format(len(gpr.index)))
# Plot the snow depths using aqua colored pixels
depths.plot(ax=ax, marker=',' , color='aqua', markersize=1, label='{} Manual Depths'.format(len(depths.index)))
# Plot the pits as magenta triangles
pits.plot(ax=ax, marker='^' , color='magenta', edgecolor='black', markersize=70, label=' {} pits'.format(len(pits.index)))
# Plot the SMP positions with a red plus
smp.plot(ax=ax, marker='+' , color='lightgreen', markersize=5, label='{} SMP Profiles'.format(len(smp.index)))
# Plot the pits as magenta triangles
cameras.plot(ax=ax, marker='s' , color='lime', edgecolor='black', markersize=70, label=' {} Camera Traps (not all shown)'.format(len(cameras.index)))
# Don't use scientific notation on the axis ticks
ax.ticklabel_format(style='plain', useOffset=False)
# Add some labeling
ax.legend(loc='lower right', fontsize='xx-large', framealpha=1.0)
ax.set_title('Overview of Measurements in the SnowEx Database for Grand Mesa', size=20)
ax.set_xlabel('Easting [m]', size=20)
ax.set_ylabel('Northing [m]', size=20)
ax.set_xlim([extent.bounds[0], extent.bounds[2]])
ax.set_ylim([extent.bounds[1], extent.bounds[3]])
# Add a block showing off details
textstr = '\n'.join([f"No. observers: {metadata['observers']}",
f"No. Instruments: {metadata['instrument']}",
f"Temporal Extent: {metadata['date'][0]} - {metadata['date'][1]}",
f"Published Datasets: {metadata['doi']}"])
props = dict(boxstyle='round', facecolor='white', alpha=0.9)
# place a text box in upper left in axes coords
s = ax.text(0.78, 0.95, textstr, transform=ax.transAxes, fontsize='xx-large',
verticalalignment='top', bbox=props)
# Save the figure if you want
#plt.savefig('/SOME/LOCATION/')
[10]:
session.close()
[ ]: