"""
SnowEx API for querying measurement data.
This module has the query logic for the Python API and the LambdaClient.
See the lambda_handler.py for requirements on exposing an endpoint to the client.
"""
import logging
import os
import pandas as pd
from sqlalchemy import Numeric, cast, exists, literal
from sqlalchemy.dialects import postgresql
from sqlalchemy.sql import func
from snowexsql.db import db_session_with_credentials
from snowexsql.tables import (
DOI,
Campaign,
CampaignObservation,
ImageData,
ImageObservation,
Instrument,
LayerData,
MeasurementType,
Observer,
PointData,
PointObservation,
Site,
)
# Initialize logger first
LOG = logging.getLogger(__name__)
def query_to_geopandas(query, engine, **kwargs):
"""
Convert SQLAlchemy query to GeoDataFrame (if geopandas available) or DataFrame.
Execution context:
- Local API: Returns GeoDataFrame with proper geometry objects
- Lambda client: Returns pandas DataFrame
- DataFrame is serialized to JSON by lambda_handler
- lambda_client receives JSON and converts to GeoDataFrame client-side
Args:
query: SQLAlchemy Query object
engine: SQLAlchemy engine
**kwargs: Additional arguments passed to read_postgis or read_sql
Returns:
Local API: geopandas.GeoDataFrame
Lambda client: pandas.DataFrame
"""
sql = query.statement.compile(dialect=postgresql.dialect())
try:
import geopandas as gpd
return gpd.read_postgis(sql, engine.connect(), **kwargs)
except ImportError:
# Geopandas not available (e.g., Lambda environment)
# Returns pandas DataFrame with geometry as WKB/WKT
# lambda_client will convert to GeoDataFrame client-side
return pd.read_sql(sql, engine, **kwargs)
def raster_to_rasterio(rasters):
"""Raster functionality requires rasterio"""
raise ImportError(
"Raster functionality not available in Lambda environment. "
"Use local API for raster operations."
)
[docs]
class LargeQueryCheckException(RuntimeError):
pass
def get_points():
"""
Get a single row from the points table
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(PointData).limit(1)
# Execute that query!
result = qry.all()
# TODO:
# * Possible enums
# * implement 'like' or 'contains' method
[docs]
class BaseDataset:
"""
Base API access class for all SnowEx measurement dataset providing common
query methods and shared attributes.
Provides filtering, querying, and property-based access to measurements
stored in the SnowEx database. Subclasses must set ``MODEL`` to the
appropriate SQLAlchemy table class.
"""
MODEL = None
ALLOWED_QRY_KWARGS = [
"campaign",
"date",
"instrument",
"type",
"utm_zone",
"date_greater_equal",
"date_less_equal",
"value_greater_equal",
"value_less_equal",
"doi",
"observer",
]
SPECIAL_KWARGS = ["limit"]
# Default max record count
MAX_RECORD_COUNT = 1000
[docs]
@staticmethod
def retrieve_single_value_result(result):
"""
When we only request a single thing we still get a list of
lists this function filters it out. This usually looks like a
list of tuples.
"""
final = []
if len(result) != 0:
final = [r[0] for r in result]
return final
@classmethod
def _build_select_clause(cls, verbose=False):
"""
Build the SELECT clause for queries.
Override in subclasses to customize columns returned based on verbose
mode.
Args:
verbose: If True, return denormalized data with related table
columns
Returns:
List of SQLAlchemy column expressions for session.query()
"""
# Default: just return the model (all its direct columns)
return [cls.MODEL]
@classmethod
def _check_size(cls, qry, kwargs):
"""
Safeguard against accidental giant requests
"""
count = qry.count()
if count > cls.MAX_RECORD_COUNT and "limit" not in kwargs:
raise LargeQueryCheckException(
f"Query will return {count} number of records,"
f" but we have a default max of {cls.MAX_RECORD_COUNT}."
f" If you want to proceed, set the 'limit' filter"
f" to the desired number of records."
)
@classmethod
def _filter_campaign(cls, qry, v):
qry = qry.filter(Site.campaign.has(Campaign.name == v))
return qry
@classmethod
def _filter_observers(cls, qry, v):
qry = qry.join(cls.MODEL.observers).filter(Observer.name == v)
return qry
@classmethod
def _filter_instrument(cls, qry, value):
return qry.filter(cls.MODEL.instrument.has(name=value))
@classmethod
def _filter_measurement_type(cls, qry, value):
return qry.join(cls.MODEL.measurement_type).filter(
MeasurementType.name == value
)
@classmethod
def _filter_doi(cls, qry, value):
return qry.join(cls.MODEL.doi).filter(DOI.doi == value)
[docs]
@classmethod
def extend_qry(cls, qry, check_size=True, **kwargs):
if cls.MODEL is None:
raise ValueError("You must use a class with a MODEL.")
# use the default kwargs
for k, v in kwargs.items():
# Handle special operations
if k in cls.ALLOWED_QRY_KWARGS:
qry_model = cls.MODEL
# Logic for filtering on date with LayerData
if "date" in k and cls.MODEL == LayerData:
qry = qry.join(LayerData.site)
qry_model = Site
elif cls.MODEL == PointData:
qry = qry.join(PointData.observation)
# standard filtering using qry.filter
if isinstance(v, list):
filter_col = getattr(qry_model, k)
if k == "date":
raise ValueError("We cannot search for a list of dates")
elif "_equal" in k:
raise ValueError(
"We cannot compare greater_equal or less_equal with a list"
)
elif k == "site":
# Skip list handling here, will be handled below
pass
else:
qry = qry.filter(filter_col.in_(v))
LOG.debug(f"Filtering {k} to value {v}")
else:
# Filter boundary
if "_greater_equal" in k:
key = k.split("_greater_equal")[0]
if key == "value":
qry = qry.filter(
cast(getattr(qry_model, key), Numeric) >= v
)
else:
qry = qry.filter(getattr(qry_model, key) >= v)
elif "_less_equal" in k:
key = k.split("_less_equal")[0]
if key == "value":
qry = qry.filter(
cast(getattr(qry_model, key), Numeric) <= v
)
else:
qry = qry.filter(getattr(qry_model, key) <= v)
# Filter linked columns
elif k == "instrument":
qry = cls._filter_instrument(qry, v)
elif k == "campaign":
qry = cls._filter_campaign(qry, v)
elif k == "site":
# Handle list of site names
if isinstance(v, list):
qry = qry.filter(
exists()
.where(qry_model.site_id == Site.id)
.where(Site.name.in_(v))
)
else:
qry = qry.filter(qry_model.site.has(name=v))
elif k == "observer":
qry = cls._filter_observers(qry, v)
elif k == "doi":
qry = cls._filter_doi(qry, v)
elif k == "type":
qry = cls._filter_measurement_type(qry, v)
# Filter to exact value
else:
qry = qry.filter(getattr(qry_model, k) == v)
LOG.debug(f"Filtering {k} to list {v}")
# to avoid limit before filter
elif k in cls.SPECIAL_KWARGS:
if k == "limit":
qry = qry.limit(v)
else:
# Error out for not-allowed kwargs
raise ValueError(f"{k} is not an allowed filter")
if check_size:
cls._check_size(qry, kwargs)
return qry
[docs]
@classmethod
def from_unique_entries(cls, columns_to_search, **kwargs):
"""
Returns unique values from a column to help with filtering
"""
columns = [getattr(cls.MODEL, column) for column in columns_to_search]
with db_session_with_credentials() as (_engine, session):
try:
qry = session.query(*columns)
# Hardcode the limit to
qry = cls.extend_qry(qry, check_size=False, **kwargs)
results = qry.distinct().all()
except Exception as e:
session.close()
LOG.error("Failed query finding options for filtering")
raise e
if len(columns_to_search) == 1:
results = cls.retrieve_single_value_result(results)
return results
[docs]
@classmethod
def from_filter(cls, verbose=False, **kwargs):
"""
Get data for the class by filtering by allowed arguments. The allowed
filters are cls.ALLOWED_QRY_KWARGS.
Args:
verbose: If True, return denormalized data with related table columns
kwargs: Filter arguments from ALLOWED_QRY_KWARGS
"""
with db_session_with_credentials() as (engine, session):
try:
select_clause = cls._build_select_clause(verbose)
qry = session.query(*select_clause)
# Add explicit joins for verbose mode to avoid cartesian
# products
if verbose and hasattr(cls, "_add_verbose_joins"):
qry = cls._add_verbose_joins(qry)
elif hasattr(cls, "_add_base_joins"):
# For verbose=False, still need basic joins (e.g.,
# Site for geom)
qry = cls._add_base_joins(qry)
qry = cls.extend_qry(qry, **kwargs)
# For debugging in the test suite and not
# recommended in production
# https://docs.sqlalchemy.org/en/20/faq/
# sqlexpressions.html#rendering-postcompile-
# parameters-as-bound-parameters
if "DEBUG_QUERY" in os.environ:
full_sql_query = qry.statement.compile(
compile_kwargs={"literal_binds": True}
)
print("\n ** SQL query **")
print(full_sql_query)
df = query_to_geopandas(qry, engine)
except Exception as e:
session.close()
LOG.error("Failed query for PointData")
raise e
return df
[docs]
@classmethod
def from_area(
cls, verbose=False, shp=None, pt=None, buffer=None, crs=26912, **kwargs
):
"""
Get data for the class within a specific shapefile or
within a point and a known buffer. Uses PostGIS functions via ORM
for spatial operations, eliminating dependency on geoalchemy2/shapely.
Args:
verbose: If True, return denormalized data with related table columns
shp: shapely geometry in which to filter, or WKT string
pt: shapely point that will have a buffer applied, or WKT string
buffer: buffer distance in same units as point (meters if using geography)
crs: integer SRID/EPSG code (default 26912 = UTM Zone 12N)
kwargs: for more filtering or limiting (cls.ALLOWED_QRY_KWARGS)
Returns:
pandas DataFrame with results (includes geom column with WKT)
"""
if shp is None and pt is None:
raise ValueError("Inputs must be a shape description or a point and buffer")
if (pt is not None and buffer is None) or (buffer is not None and pt is None):
raise ValueError("pt and buffer must be given together")
# Convert shapely objects to WKT if needed
if shp is not None and hasattr(shp, "wkt"):
shp_wkt = shp.wkt
elif isinstance(shp, str):
shp_wkt = shp
else:
shp_wkt = None
if pt is not None:
if hasattr(pt, "wkt"):
pt_wkt = pt.wkt
elif isinstance(pt, str):
pt_wkt = pt
elif isinstance(pt, (tuple, list)) and len(pt) == 2:
# Handle (x, y) tuple format
pt_wkt = f"POINT ({pt[0]} {pt[1]})"
else:
pt_wkt = None
else:
pt_wkt = None
# Determine table structure
table_name = cls.MODEL.__tablename__
needs_site_join = table_name == "layers"
with db_session_with_credentials() as (engine, session):
try:
# Detect database SRID to avoid transforming indexed column
# Query first non-null geometry to determine database SRID
if needs_site_join:
srid_qry = (
session.query(func.ST_SRID(Site.geom))
.filter(Site.geom.isnot(None))
.limit(1)
)
else:
srid_qry = (
session.query(func.ST_SRID(cls.MODEL.geom))
.filter(cls.MODEL.geom.isnot(None))
.limit(1)
)
try:
db_srid_result = session.execute(srid_qry).first()
if not db_srid_result or db_srid_result[0] is None:
# No data in table yet - use input CRS as default
# This allows empty table queries to work
# (will return empty)
LOG.warning(
f"No geometries found in {table_name}, "
f"using input CRS {crs} as default"
)
db_srid = crs
else:
db_srid = db_srid_result[0]
LOG.debug(
f"Detected database SRID: \
{db_srid} for table {table_name}"
)
except Exception as srid_error:
# If SRID detection fails, fall back to input CRS
LOG.warning(
f"SRID detection failed for {table_name}: {srid_error}."
f"Using input CRS {crs} as default"
)
db_srid = crs
# Build PostGIS search geometry
# Transform search geometry to match database SRID for index
# usage
if pt_wkt:
# Create point in input CRS, buffer it, then transform to
# DB SRID
# Buffer before transform to ensure correct distance units
search_geom = func.ST_Transform(
func.ST_Buffer(
func.ST_GeomFromText(literal(pt_wkt), literal(crs)),
literal(buffer),
),
literal(db_srid),
)
elif shp_wkt:
# Transform shape from input CRS to database SRID
search_geom = func.ST_Transform(
func.ST_GeomFromText(literal(shp_wkt), literal(crs)),
literal(db_srid),
)
else:
raise ValueError("Unable to parse geometry input")
# Build select clause (handles verbose mode)
select_clause = (
cls._build_select_clause(verbose)
if hasattr(cls, "_build_select_clause")
else [cls.MODEL]
)
qry = session.query(*select_clause)
# Add explicit joins for verbose mode to avoid
# cartesian products
if verbose and hasattr(cls, "_add_verbose_joins"):
qry = cls._add_verbose_joins(qry)
elif hasattr(cls, "_add_base_joins"):
# For verbose=False, still need basic joins
# (e.g., Site for geom)
qry = cls._add_base_joins(qry)
# Add spatial filter
if needs_site_join:
# For LayerData, join to Site for geometry
# SQLAlchemy handles duplicate joins from _add_*_joins above
qry = qry.join(cls.MODEL.site)
qry = qry.filter(func.ST_Intersects(Site.geom, search_geom))
else:
# For PointData, use direct geometry column
qry = qry.filter(func.ST_Intersects(cls.MODEL.geom, search_geom))
# Add standard filters using existing extend_qry
# This handles type, instrument, campaign, date ranges, etc.
qry = cls.extend_qry(qry, **kwargs)
# Execute and convert to GeoDataFrame
df = query_to_geopandas(qry, engine)
except Exception as e:
session.close()
LOG.error(f"Failed query for {cls.__name__}")
raise e
return df
@property
def all_campaigns(self):
"""
Return all campaign names
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(Campaign.name).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_types(self):
"""
Return all types of the data
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(MeasurementType.name).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_dates(self):
"""
Return all distinct dates in the data
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(self.MODEL.date).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_observers(self):
"""
Return all distinct observers in the data
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(Observer.name).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_dois(self):
"""
Return all distinct DOIs in the data
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(DOI.doi).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_units(self):
"""
Return all distinct units in the data
"""
with db_session_with_credentials() as (_engine, session):
qry = session.query(self.MODEL.units).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_instruments(self):
"""
Return all distinct instruments in the data
"""
with db_session_with_credentials() as (_engine, session):
# Use EXISTS for better performance on large datasets
# (29GB+ tables)
qry = session.query(Instrument.name).filter(
exists().where(self.MODEL.instrument_id == Instrument.id)
)
result = qry.all()
return self.retrieve_single_value_result(result)
[docs]
class PointMeasurements(BaseDataset):
"""
API class for access to PointData from a single location.
Examples: Snow Depth, SWE
"""
MODEL = PointData
@classmethod
def _build_select_clause(cls, verbose=False):
"""
Build SELECT clause for PointMeasurements queries.
Args:
verbose: If False, return only core point columns.
If True, return denormalized data with observation,
instrument, and measurement type info.
"""
if verbose:
# Return denormalized data with meaningful column names
# Note: Must also join these tables explicitly in the query
return [
cls.MODEL.value,
cls.MODEL.datetime.label("date"),
cls.MODEL.elevation,
cls.MODEL.geom,
CampaignObservation.name.label("observation_name"),
CampaignObservation.description.label("obs_description"),
MeasurementType.name.label("type"),
MeasurementType.units.label("units"),
MeasurementType.derived.label("derived"),
Instrument.name.label("instrument_name"),
Instrument.model.label("instrument_model"),
Instrument.specifications.label("instrument_specifications"),
Campaign.name.label("campaign_name"),
Observer.name.label("observer_name"),
]
else:
# Return only columns from points table
return [
cls.MODEL.id,
cls.MODEL.value,
cls.MODEL.datetime,
cls.MODEL.elevation,
cls.MODEL.geom,
]
@classmethod
def _add_verbose_joins(cls, qry):
"""
Add explicit joins needed for verbose mode to avoid cartesian products.
"""
qry = qry.join(cls.MODEL.observation)
qry = qry.join(cls.MODEL.measurement_type)
qry = qry.join(PointObservation.instrument)
qry = qry.join(PointObservation.campaign)
qry = qry.join(PointObservation.observer)
return qry
@classmethod
def _filter_campaign(cls, qry, value):
return (
qry.join(cls.MODEL.observation)
.join(PointObservation.campaign)
.filter(Campaign.name == value)
)
@classmethod
def _filter_instrument(cls, qry, value):
return (
qry.join(cls.MODEL.observation)
.join(PointObservation.instrument)
.filter(Instrument.name == value)
)
@classmethod
def _filter_doi(cls, qry, value):
return (
qry.join(cls.MODEL.observation)
.join(PointObservation.doi)
.filter(DOI.doi == value)
)
@classmethod
def _filter_observers(cls, qry, value):
return (
qry.join(cls.MODEL.observation)
.join(PointObservation.observer)
.filter(Observer.name == value)
)
@property
def all_types(self):
"""
Return all measurement types that have data in the points table
"""
with db_session_with_credentials() as (_engine, session):
# Use EXISTS for better performance on large points table
qry = session.query(MeasurementType.name).filter(
exists().where(PointData.measurement_type_id == MeasurementType.id)
)
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_instruments(self):
"""
Return all distinct instruments in the data
"""
with db_session_with_credentials() as (_engine, session):
result = (
session.query(Instrument.name)
.filter(
Instrument.id.in_(
session.query(PointObservation.instrument_id).distinct()
)
)
.distinct()
.all()
)
return self.retrieve_single_value_result(result)
class TooManyRastersException(Exception):
"""
Exception to report to users that their query will produce too
many rasters
"""
[docs]
class LayerMeasurements(BaseDataset):
"""
API class for access to LayerData.
Examples: Snow pit observations (depth, density, temperature), SMP
"""
MODEL = LayerData
ALLOWED_QRY_KWARGS = [
"campaign",
"site",
"date",
"instrument",
"observer",
"type",
"utm_zone",
"date_greater_equal",
"date_less_equal",
"doi",
"value_greater_equal",
"value_less_equal",
]
@classmethod
def _filter_campaign(cls, qry, v):
return qry.join(cls.MODEL.site).join(Site.campaign).filter(Campaign.name == v)
@classmethod
def _filter_observers(cls, qry, v):
return qry.join(cls.MODEL.site).join(Site.observers).filter(Observer.name == v)
@classmethod
def _filter_doi(cls, qry, value):
return qry.join(cls.MODEL.site).join(Site.doi).filter(DOI.doi == value)
@classmethod
def _build_select_clause(cls, verbose=False):
"""
Build SELECT clause for LayerMeasurements queries.
Args:
verbose: If False, return only core layer columns (no joins).
If True, return denormalized data with site, instrument,
and measurement type info.
"""
if verbose:
# Return denormalized data with meaningful column names
# Note: Must also join these tables explicitly in the query
return [
cls.MODEL.depth,
cls.MODEL.bottom_depth,
cls.MODEL.value,
Site.name.label("site_name"),
Site.description.label("site_description"),
Site.slope_angle.label("slope_angle"),
Site.aspect.label("aspect"),
Site.air_temp.label("air_temp"),
Site.total_depth.label("total_depth"),
Site.weather_description.label("weather_description"),
Site.precip.label("precip"),
Site.sky_cover.label("sky_cover"),
Site.wind.label("wind"),
Site.ground_condition.label("ground_condition"),
Site.ground_roughness.label("ground_roughness"),
Site.ground_vegetation.label("ground_vegetation"),
Site.vegetation_height.label("vegetation_height"),
Site.tree_canopy.label("tree_canopy"),
Site.comments.label("site_comments"),
Site.datetime.label("date"),
Site.geom,
MeasurementType.name.label("type"),
MeasurementType.units.label("units"),
MeasurementType.derived.label("type_derived"),
Instrument.name.label("instrument_name"),
Instrument.model.label("instrument_model"),
Instrument.specifications.label("instrument_specifications"),
func.ST_AsText(Site.geom).label("geom_wkt"),
]
else:
# Return only core columns from layers table plus geom for geopandas
return [
cls.MODEL.id,
cls.MODEL.depth,
cls.MODEL.bottom_depth,
cls.MODEL.value,
Site.geom, # Required for GeoDataFrame
]
@classmethod
def _add_base_joins(cls, qry):
"""
Add minimal joins needed for non-verbose mode (just Site for geom).
"""
qry = qry.join(cls.MODEL.site)
return qry
@classmethod
def _add_verbose_joins(cls, qry):
"""
Add explicit joins needed for verbose mode to avoid cartesian products.
"""
qry = qry.join(cls.MODEL.site)
qry = qry.join(cls.MODEL.measurement_type)
qry = qry.join(cls.MODEL.instrument)
return qry
@property
def all_types(self):
"""
Return all measurement types that have data in the layers table
"""
with db_session_with_credentials() as (_engine, session):
# Use EXISTS for better performance on 208M row table
qry = session.query(MeasurementType.name).filter(
exists().where(LayerData.measurement_type_id == MeasurementType.id)
)
result = qry.all()
return self.retrieve_single_value_result(result)
@property
def all_sites(self):
"""
Return all specific site names
"""
with db_session_with_credentials() as (_engine, session):
result = session.query(Site.name).distinct().all()
return self.retrieve_single_value_result(result)
@property
def all_dates(self):
"""
Return all distinct dates in the data
"""
with db_session_with_credentials() as (_engine, session):
result = session.query(Site.date).distinct().all()
return self.retrieve_single_value_result(result)
@property
def all_units(self):
"""
Return all distinct units in the data
"""
with db_session_with_credentials() as (_engine, session):
result = session.query(MeasurementType.units).distinct().all()
return self.retrieve_single_value_result(result)
[docs]
@classmethod
def get_sites(cls, site_names=None, **kwargs):
"""
Get site information including geometries.
Args:
site_names: List of site names or single site name
kwargs: Additional filters (campaign, date, etc.)
Returns:
GeoDataFrame with site information
"""
with db_session_with_credentials() as (engine, session):
qry = session.query(
Site.name, Site.geom, Site.description, Site.datetime
).distinct()
# others can be added
if site_names is not None:
if isinstance(site_names, list):
qry = qry.filter(Site.name.in_(site_names))
else:
qry = qry.filter(Site.name == site_names)
df = query_to_geopandas(qry, engine)
return df
class RasterMeasurements(BaseDataset):
"""
Early beta class for handling data stored in raster datasets.
"""
MODEL = ImageData
ALLOWED_QRY_KWARGS = BaseDataset.ALLOWED_QRY_KWARGS + ["description"]
@property
def all_types(self):
"""
Return all measurement types that have data in the images table
"""
with db_session_with_credentials() as (_engine, session):
result = (
session.query(MeasurementType.name)
.join(ImageData, ImageData.observation_id == ImageObservation.id)
.join(ImageObservation.measurement_type)
.distinct()
.all()
)
return self.retrieve_single_value_result(result)
@property
def all_descriptions(self):
with db_session_with_credentials() as (_engine, session):
qry = session.query(self.MODEL.description).distinct()
result = qry.all()
return self.retrieve_single_value_result(result)
@classmethod
def check_for_single_dataset(cls, **kwargs):
"""
At the moment there is not a clear path to how to deal with
multiple rasters so check that the user only requested one
dataset
"""
LOG.info("Checking raster query for single raster dataset...")
multi_raster_indicators = [
"instrument",
"date",
"observers",
"doi",
"type",
"description",
]
with db_session_with_credentials() as (_engine, session):
try:
# Form query and check if the query spans
# multiple rasters
for column in multi_raster_indicators:
values = cls.from_unique_entries([column], **kwargs)
if len(values) > 1:
options = [f"'{v}'" for v in values]
raise TooManyRastersException(
f"More than one `{column}` suggests "
f"there are multiple raster datasets. "
f"Try filter {column} to one of the "
f"following values {', '.join(options)}."
)
except Exception as e:
session.close()
LOG.error("Failed query for Raster Data")
raise e
@classmethod
def from_filter(cls, **kwargs):
"""
Get data for the class by filtering by allowed arguments.
The allowed filters are cls.ALLOWED_QRY_KWARGS.
"""
cls.check_for_single_dataset(**kwargs)
with db_session_with_credentials() as (_engine, session):
try:
# Rebuild the query and form the raster
base_query = cls.MODEL.raster
qry = session.query(
func.ST_AsTiff(func.ST_Union(base_query, type_=Raster))
)
qry = cls.extend_qry(qry, **kwargs)
rasters = qry.all()
# Get the rasterio object of the raster
datasets = raster_to_rasterio(rasters)
except Exception as e:
LOG.error("Failed query for Raster Data")
raise e
return datasets
@classmethod
def from_area(cls, shp=None, pt=None, buffer=None, crs=26912, **kwargs):
if shp is None and pt is None:
raise ValueError("We need a shape description or a point and buffer")
if (pt is not None and buffer is None) or (buffer is not None and pt is None):
raise ValueError("pt and buffer must be given together")
# Check if geometric operations are available
if not HAS_CORE_GIS:
raise ImportError(
"geoalchemy2 and shapely are required for "
"geometric filtering. Install with: "
"pip install geoalchemy2 shapely"
)
with db_session_with_credentials() as (_engine, session):
try:
# Get shape ready for cropping with rasters
if shp:
db_shp = from_shape(shp, srid=crs)
else:
qry_pt = from_shape(pt)
qry = session.query(
gfunc.ST_SetSRID(func.ST_Buffer(qry_pt, buffer), crs)
)
db_shp = qry.all()[0][0]
# Grab the rasters, union and clip them
base_query = func.ST_AsTiff(
func.ST_Clip(
func.ST_Union(ImageData.raster, type_=Raster), db_shp, True
)
)
q = session.query(base_query)
# Find all the tiles that
q = q.filter(gfunc.ST_Intersects(ImageData.raster, db_shp))
limit = kwargs.get("limit")
if limit:
kwargs.pop("limit")
q = cls.extend_qry(q, check_size=False, **kwargs)
rasters = q.all()
# Get the rasterio object of the raster
datasets = raster_to_rasterio(rasters)
if len(datasets) > 0:
dataset = datasets[0]
else:
dataset = datasets
return dataset
except Exception as e:
LOG.error("Failed query for Raster Data")
raise e
