intro-geopands.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to geopandas for analyzing spatial data\n",
    "\n",
    "* NOTE: this is a slight 2024 update to material originally from https://github.com/uwescience/dssg-pandas-sql-tutorial\n",
    "\n",
    "[*Geopandas*](https://github.com/geopandas/geopandas) is an excellent extension to the [*Pandas*](https://github.com/pandas-dev/pandas) package for working with tabular data that is georeferenced - for example, Points, Lines, and Polygons that have associated attributes. Geopandas builds upon other great libraries including [*Shapely*](https://github.com/Toblerity/Shapely) and [*Fiona*](https://github.com/Toblerity/Fiona).\n",
    "\n",
    "There are a lot of tutorials out there: Check out this lesson from a [eScience DSSG tutorial](https://uwescience.github.io/SQL-geospatial-tutorial/) or the [eScience Geohackweek tutorial](https://geohackweek.github.io/vector/).\n",
    "\n",
    "## This notebook covers the following topics:\n",
    "\n",
    "* dataset exploration\n",
    "* basic statistical analysis\n",
    "* simple geometric operations (convex_hull)\n",
    "* exporting data for GIS software\n",
    "* converting coordinates and calculating distance and area\n",
    "* attribute and spatial joins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import all the libraries we are going to use\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "\n",
    "# plotting configuration\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.rcParams['font.size'] = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import warnings\n",
    "#warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First thing to note is that many geospatial python packages are under active development\n",
    "# and updated regularly. So the version you're using today will likely be updated in the near \n",
    "# future!\n",
    "print('Pandas version: ', pd.__version__) #0.24.2\n",
    "print('Geopandas version: ', gpd.__version__) #0.5.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dataset\n",
    "\n",
    "We'll work with 2 CSV files from here\n",
    "https://github.com/uwescience/dssg-pandas-sql-tutorial/tree/master/data\n",
    "\n",
    "\n",
    "**Note that the 'places' information has 9 columns with labels in the first row**\n",
    "!head ./data/Places_Full.csv\n",
    "\n",
    "name,address,city,lat,lng,place_id,rating,class,type\n",
    "Trader Joe's,1700 East Madison Street,Seattle,47.6158665,-122.3099133,ChIJx0M1ztNqkFQRtgspEllQxk8,4.5,supermarket,\n",
    "\n",
    "**distance data has information about transit between two places**\n",
    "!head ./data/Dist_Out.csv\n",
    "class,departure_time,distance.value,duration.value,duration_in_traffic.value,end_lat,end_lon,fare.value,mode,origin,pair,place_id,rank,start_lat,start_lon,status\n",
    "hospital,2018-06-06 12:41:26,1.47545059,32.18333333,,47.7160344,-122.3122296,,walking,530330001003,530330001003-ChIJEerXpl0RkFQRDOjfwBQDzlA,ChIJEerXpl0RkFQRDOjfwBQDzlA,1,47.72281,-122.28945,OK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pandas review"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's work with the smaller 'Places_Full.csv' first\n",
    "# All Pandas methods are accessed via the 'pd' attribute.\n",
    "# Since the file is well-formatted, it is easily read into memory:\n",
    "# We can read small datasets directly from a URL\n",
    "filePath = 'https://github.com/uwescience/dssg-pandas-sql-tutorial/raw/master/data/Places_Full.csv'\n",
    "df = pd.read_csv(filePath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 'df' stands for 'Data Frame'. It is essentially a spreadsheet:\n",
    "# HTML formatted output makes for easier reading\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The dataframe has a lot of convenient methods for fast data exploration\n",
    "# Start with info to confirm that things were read in correctly\n",
    "df.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simple statistics are obtained for numerical columns\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Columns can be accessed as dictionary items:\n",
    "print(df['city'].unique())\n",
    "# Or accessed as attributes for faster typing:\n",
    "print(df.city.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gotcha: \n",
    "# 'class' attribute is reserved for internal use, but it is also a column heading\n",
    "#print(df.class.unique()) # This causes an error\n",
    "print(df['class'].unique()) # This is ok"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's change the name of class to avoid that error\n",
    "df.rename(columns={'class':'place_class'}, inplace=True)\n",
    "df.place_class.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Another common issue with tabular data - \n",
    "# certain measurements don't always fit into the defined columns \n",
    "# or are missing data, and therefore filled with 'not-a-number (nan)'\n",
    "# For example, some entries don't have a listed city:\n",
    "print(df.city.unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract data entries without a city\n",
    "dfNan = df[df.city.isna()]\n",
    "dfNan.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise: print top 10 places, sorted by rating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pandas has some convenient built-in plotting methods using matplotlib\n",
    "# NOTE: semi-colon prevents bins and data objects from being returned\n",
    "df.rating.hist(bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pandas has some convenient built-in plotting methods using matplotlib\n",
    "# NOTE: semi-colon prevents bins and data objects from being returned\n",
    "df.boxplot('rating');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise: plot box-plots showing 'rating' for each 'place_class'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can work with latitude and longitude infrmation in pandas, but it doesn't integrate well\n",
    "# with other geospatial tools for analysis and visualization, for example:\n",
    "df.plot.scatter(x='lng', y='lat', c='rating', cmap='viridis');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From Pandas to Geopandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Commonly, there at latitude, longitude columns for places, but we must tell geopandas \n",
    "# explicity what the geometries (Points) and coordinate reference is (EPSG:4326)\n",
    "# https://geopandas.readthedocs.io/en/latest/projections.html#coordinate-reference-systems\n",
    "\n",
    "crs = 'EPSG:4326'\n",
    "gf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lng, df.lat), crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: we've added the special 'geometry' columns\n",
    "gf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization\n",
    "\n",
    "* GeoJSON files are rendered on maps on GitHub.org:\n",
    "https://github.com/uwescience/dssg2018-geopandasSQL-tutorial/blob/master/data/places.geojson\n",
    "\n",
    "* There are also some Python plotting libraries that work well with Geopandas including [*cartopy*](https://github.com/SciTools/cartopy), and [*ipyleaflet*](https://github.com/jupyter-widgets/ipyleaflet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Geopandas also has some convenient built-in plotting methods\n",
    "# Note that we no longer have to specify lat and lon\n",
    "gf.plot(c=gf.rating);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Geopands 'explore' is awesome\n",
    "gf.explore(column='rating')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometric operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We also now have access to lots of quantifiable spatial information. For example, how large\n",
    "# is the area containing all these points?\n",
    "#http://geopandas.org/geometric_manipulations.html\n",
    "point_collection = shapely.geometry.MultiPoint(gf.geometry.tolist())\n",
    "polygon = point_collection.convex_hull\n",
    "polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perimeter = polygon.boundary\n",
    "perimeter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can save these geometries in a separte geodataframe as you go:\n",
    "gfShape = gpd.GeoDataFrame(geometry=[polygon], crs = 'EPSG:4326')\n",
    "gfShape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subsetting is very easy in Geopandas. Often we only want points in a certain bounding box\n",
    "ymin, ymax, xmin, xmax = [47.55, 47.65, -122.4, -122.3]\n",
    "subset = gf.cx[xmin:xmax, ymin:ymax]\n",
    "print(f'original entries = {len(gf)}, subsetted entries = {len(subset)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subset.explore()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other free geospatial tools\n",
    "\n",
    "Geographic Information Systems (GIS) are designed for working with geospatial data. If geopandas is lacking, consider using [QGIS](https://qgis.org/en/site/). Also, [Google Earth Pro](https://www.google.com/earth/desktop/) is now free and is a great visualization tool. Finally, [GDAL/OGR](http://www.gdal.org) is a powerful library with command line tools that many other software packages are based upon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you want to share the geometry / export it to a GIS program for further analysis:\n",
    "\n",
    "# There are much more powerful ways to explore geospatial data on interactive maps:\n",
    "# Uploading this file to GitHub will automatically render on a map\n",
    "gf.to_file('./data/places.geojson', driver='GeoJSON')\n",
    "\n",
    "# Default is ESRI shapefile, but Geojson or Geopackage, etc should work to\n",
    "#http://geopandas.org/io.html\n",
    "#http://www.geopackage.org\n",
    "\n",
    "#gfShape.to_file('./data/myshape.shp') #ESRI shapefile is default\n",
    "gfShape.to_file('./data/myshape.gpkg', driver='GPKG')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distance, area, etc.\n",
    " \n",
    "We commonly have points in latitude and longitude, but want to know distances on the ground in meters. These conversions are non-trivial and require some knowledge of different coordinate systems and map projections. Be carefule when making these calculations!\n",
    "https://support.esri.com/en/technical-article/000011356"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Excercise: Determine the area of gfShape, what are the unites?\n",
    "gfShape.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1) convert to a local coordinate system in metric units of distance!\n",
    "\n",
    "# Google Mercator (EPSG:3857 - https://epsg.io/3857) is another popular \n",
    "# projection for web maps, but it is \"direction preserving\", not \"area preserving!\"\n",
    "\n",
    "# A good choice for local distances is the Universal Transverse Mercator\n",
    "# projection. For Seattle it's EPSG:32610\n",
    "# http://spatialreference.org/ref/epsg/wgs-84-utm-zone-10n/\n",
    "\n",
    "gfShape.to_crs('epsg:32610', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# will be in units of km^2\n",
    "area = gfShape.area.values[0] * 1e-6\n",
    "print(f'Points are within an area of {area:.1f} km^2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exercise:\n",
    "# Calculate the length of the perimeter\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load a second database of trip distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfD = gpd.pd.read_csv('https://github.com/uwescience/dssg-pandas-sql-tutorial/raw/master/data/Dist_Out.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfD.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfD.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at attributes (columns) of first entry (row)\n",
    "dfD.iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 'distance.value in the table is route distance, what about straight-line distance?\n",
    "# we already have destination (place_id) lon,lat points in the earlier geodataframe 'gf'\n",
    "row = 0\n",
    "dest = dfD.place_id.iloc[row]\n",
    "\n",
    "gf.query('place_id == @dest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# So let's use the start_lat, start_lon for the points in this second geodataframe\n",
    "crs = 'EPSG:4326'\n",
    "gfD = gpd.GeoDataFrame(dfD, geometry=gpd.points_from_xy(dfD.start_lon, dfD.start_lat), crs=crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To get distance let's convert once again to UTM\n",
    "# Or UTM EPSG:32610\n",
    "gf_merc = gf.to_crs('epsg:32610')\n",
    "gfD_merc = gfD.to_crs('epsg:32610')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's do a sanity check for distance for this first pair:\n",
    "pointDest = gf_merc.query('place_id == @dest').geometry.values[0]\n",
    "pointDest.coords.xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pointOrig = gfD_merc.geometry.iloc[0]\n",
    "pointOrig.coords.xy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_m = pointDest.distance(pointOrig)\n",
    "print(f'Distance in meters: {dist_m:.1f}')\n",
    "dist_mi = dist_m/1609.34\n",
    "print(f'Distance in miles: {dist_mi:.1f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# measuring the straight line distance in google earth i get (1.87 km, 1.16 miles)\n",
    "# Does this make sense? if this is a walking trip that the reported distances would be a bit larger!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Excercise: save the Start Point, End Point, and Straight Line between\n",
    "# them to view in Google Earth or QGIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# That was kind of complicated for one point, but it's now easy to get all the \n",
    "# straight-line distances between origins and a particular destination:\n",
    "pId = 'ChIJEerXpl0RkFQRDOjfwBQDzlA'\n",
    "pointDest = gf_merc.query('place_id == @pId').geometry.values[0]\n",
    "gfDest = gfD_merc.query('place_id == @pId')\n",
    "gfDest\n",
    "# NOTE: lots of the same origin-id. Are these duplicates? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfDest.loc[:,'straightDist'] = gfDest.distance(pointDest).values/1609.34\n",
    "gfDest.loc[:, ['distance.value','straightDist']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Joining databases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Given these two databases share a common column (place_id), we could do an\n",
    "# \"attribute\" join\n",
    "# NOTE that the spatial information doesn't really matter in this case!\n",
    "# http://geopandas.org/mergingdata.html\n",
    "\n",
    "# Maybe we're only interested in 'departure time, distance and duration'\n",
    "# for a given destination, but want to keep the location and other attributes:\n",
    "\n",
    "columns = ('place_id','departure_time','distance.value','duration.value')\n",
    "gfSubset = gfD.loc[:,columns]\n",
    "gfSubset\n",
    "gfPlace = gfSubset.merge(gf, on='place_id')\n",
    "gfPlace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spatial merging uses spatial relationships (“intersects”, “within” or “contains”)\n",
    "# For example, let's extract points from the database for a particular neighborhood\n",
    "# See also:\n",
    "#https://github.com/seattleio\n",
    "#https://github.com/geopandas/geopandas/blob/master/examples/spatial_joins.ipynb\n",
    "\n",
    "URL = 'https://raw.githubusercontent.com/seattleio/seattle-boundaries-data/master/data/neighborhoods.geojson'\n",
    "gfN = gpd.read_file(URL)\n",
    "gfN.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gfN.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neighborhood = 'University District'\n",
    "gfNeighborhood = gfN.query('nhood == @neighborhood')\n",
    "gfNeighborhood.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nhood_places = gpd.sjoin(gf, gfNeighborhood, how=\"inner\", predicate='within')\n",
    "nhood_places"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE how columns with the same row were automatically dealt with by \n",
    "# appending a _left or _right suffix!\n",
    "print('{} Places in {}\\n'.format(len(nhood_places), neighborhood))\n",
    "print(nhood_places.name_left)\n",
    "\n",
    "# Excercise:\n",
    "# Sort by top 50 ratings in university district (or another district) and save results to both csv and json files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lots of directions to take this in!\n",
    "- Which neighborhood has the longest travel times?\n",
    "- What about ratio of straight-line distance to traveled distance?\n",
    "- Export data for a particular neighborhood\n",
    "- Incorporate another GIS dataset of your choosing!\n",
    "\n",
    "# Good luck from here!"
   ]
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}