{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Forming Queries: PostGIS Functions\n", "\n", "PostGIS offer a host of functions that we can access through python using special functions to utilize them \n", "\n", "Don't forget your [cheat sheets](https://snowexsql.readthedocs.io/en/latest/cheat_sheet.html)! \n", "\n", "\n", "In general they follow the convention\n", "``` sql\n", "ST_\n", "```\n", "\n", "They also tend to fall into (generally) 2 categories, points and rasters. \n", "\n", "\n", "## Process \n", "### Get Connected" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import the function to get connect to the db\n", "from snowexsql.db import get_db\n", "\n", "# This is what you will use for all of hackweek to access the db\n", "db_name = 'snow:hackweek@db.snowexdata.org/snowex'\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Let's get a single raster tile\n", "\n", "Checkout the documentation for [`ST_AsTiff`](https://postgis.net/docs/RT_ST_AsTIFF.html)\n", "\n", "Raster data in the database is stored in Well Known binary format so to make it useful to us we convert to geotiff format. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from snowexsql.data import ImageData\n", "\n", "# Grab a session\n", "engine, session = get_db(db_name)\n", "\n", "# What will this return?\n", "result = session.query(ImageData.raster).limit(1).all()\n", "\n", "print(type(result[0][0]))\n", "\n", "session.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# import this to use define sql functions (e.g. postgis!)\n", "from sqlalchemy.sql import func \n", "\n", "# Import this to convert to a rasterio object for easy plotting\n", "from snowexsql.conversions import raster_to_rasterio \n", "\n", "# Import a convenient function to plot with \n", "from rasterio.plot import show" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Grab a session\n", "engine, session = get_db(db_name)\n", "\n", "# Remember in the query parentheses is what we get back, in this case were asking for the raster data as a geotiff\n", "result = session.query(func.ST_AsTiff(ImageData.raster)).filter(ImageData.type == 'depth').limit(1).all()\n", "\n", "# Now make it more available as a python object \n", "datasets = raster_to_rasterio(session, result)\n", "\n", "# Plot the georeferenced image \n", "show(datasets[0], vmax=1.2, vmin=0, cmap='winter')\n", "\n", "# Close the dataset\n", "datasets[0].close()\n", "session.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lets use a few more. \n", "\n", "Lets try to get a raster tile on a pit!\n", "\n", "Checkout the documentation for \n", "\n", "* [`ST_Union`](https://postgis.net/docs/RT_ST_Union.html)\n", "* [`ST_Intersects`](https://postgis.net/docs/RT_ST_Intersects.html)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# import our pits metadata table class\n", "from snowexsql.data import SiteData\n", "from geoalchemy2.types import Raster\n", "import geoalchemy2.functions as gfunc\n", "\n", "# Grab a session\n", "engine, session = get_db(db_name)\n", "\n", "# session.rollback()\n", "\n", "# 1. Lets choose a site we want to grab a raster tile\n", "site_id = '5S31'\n", "\n", "# 2. Get the location of the pit, POSTGIS functions like to work in the text format of things so convert the point geom to text which is also in binary in the db \n", "point = session.query(SiteData.geom).filter(SiteData.site_id == site_id).distinct().all()[0][0]\n", "\n", "# 3. Merge all the tiles together, note gfunc vs func. This is because ST_Union exists in two places in postgis for geom and rasters!\n", "base = gfunc.ST_Union(ImageData.raster, _type=Raster)\n", "\n", "# 4. Get the merged result as a geotiff! \n", "base = func.ST_AsTiff(base)\n", "\n", "# 5. Filter by uavsar interferogram data\n", "qry = session.query(base).filter(ImageData.type == 'insar interferogram real')\n", "\n", "# 6. Filter by a polarization in the description \n", "qry = qry.filter(ImageData.description.contains('Polarization = HH'))\n", "\n", "# 7. Isolate tiles touching the pit location\n", "qry = qry.filter(func.ST_Intersects(ImageData.raster, point))\n", "\n", "print(qry.count())\n", "\n", "# 8. Execute, convert and plot! \n", "result = qry.all()\n", "datasets = raster_to_rasterio(session, result)\n", "show(datasets[0], vmin=-0.02, vmax=0.02, cmap='Purples')\n", "\n", "session.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Discussion**\n", "\n", "* What is a fundamental difference do you see in using `ST_Union` vs `ST_Intersects`\n", "* Did you notice `ST_Union` used `gfunc.` instead of `func.` ? How many `ST_Union`'s exist? \n", "\n", "\n", "* [`RT_ST_Union`](https://postgis.net/docs/RT_ST_Union.html)\n", "* [`ST_Union`](https://postgis.net/docs/ST_Union.html)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Lets work with some points and Postgis\n", "These functions are critical to rasters use with the database. But there are plenty of very useful functions for non-raster data too! A common use is to grab points in a certain geometry of a locations like a pit.\n", "\n", "Lets pick a pit and grab data with a certain radius of that pit using postgis functions.\n", "\n", "Checkout the documentation on:\n", "\n", "[`ST_Buffer`](https://postgis.net/docs/ST_Buffer.html)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from snowexsql.data import PointData\n", "from snowexsql.conversions import query_to_geopandas\n", "import matplotlib.pyplot as plt \n", "\n", "\n", "# Pick a pit ID\n", "site_id = '1N3'\n", "\n", "# Pick a distance around the pit to collect data in meters\n", "buffer_dist = 50\n", "\n", "# Grab a session\n", "engine, session = get_db(db_name)\n", "\n", "# Grab our pit location by provided site id from the site details table\n", "qry = session.query(SiteData.geom).filter(SiteData.site_id == site_id)\n", "\n", "# convert qry to df for easy plotting \n", "site_df = query_to_geopandas(qry, engine)\n", "\n", "# Also execute it for the normal db usage\n", "site_geom = qry.all()[0][0]\n", "\n", "# Create a polygon buffered by our distance centered on the pit\n", "qry = session.query(func.ST_Buffer(site_geom, buffer_dist))\n", "\n", "# Execute for other querying\n", "buffered_pit = qry.all()[0][0]\n", "\n", "# Filter by the dataset type depth\n", "qry = session.query(PointData).filter(PointData.type == 'depth').filter(PointData.instrument.in_(['magnaprobe','mesa']))\n", "\n", "# Grab all the point data in the buffer\n", "qry = qry.filter(func.ST_Within(PointData.geom, buffered_pit))\n", "df = query_to_geopandas(qry, engine)\n", "\n", "session.close()\n", "\n", "# plot it with style!\n", "fig, ax = plt.subplots(figsize=(8,8))\n", "ax = df.plot(ax=ax, column='value', cmap='cool')\n", "site_df.plot(ax=ax, marker='^', markersize=100, color='green')\n", "ax.legend([\"Depth\", \"Pit\"])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Recap\n", "\n", "Postgis functions are awesome but can be finicky. So go slow with them.\n", "\n", "**You should know**\n", "* Where to find PostGIS functions \n", "* When to use geoalchemy2 over sqlachemy functions call \n", "* How to chain together commands " ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.4" } }, "nbformat": 4, "nbformat_minor": 4 }