#!/usr/bin/env python
# coding: utf-8

# # GeoPandas: Pandas + geometry data type + custom geo goodness
# [Emilio Mayorga](https://github.com/emiliom/), 2017-2-6

# ## 1. Background

# [GeoPandas](http://geopandas.org) adds a spatial geometry data type to `Pandas` and enables spatial operations on these types, using [shapely](http://toblerity.org/shapely/). GeoPandas leverages Pandas together with several core open source geospatial packages and practices to provide a uniquely simple and convenient framework for handling geospatial feature data, operating on both geometries and attributes jointly, and as with Pandas, largely eliminating the need to iterate over features (rows). Also as with Pandas, it adds a very convenient and fine-tuned plotting method, and read/write methods that handle multiple file and "serialization" formats.

# **_NOTES:_**
# - Like `shapely`, these spatial data types are limited to discrete entities/features and do not address continuously varying rasters or fields.
# - While GeoPandas spatial objects can be assigned a Coordinate Reference System (`CRS`), operations can not be performed across CRS's. Plus, geodetic ("unprojected", lat-lon) CRS are not handled in a special way; the area of a geodetic polygon will be in degrees.

# GeoPandas is still young, but it builds on mature and stable and widely used packages (Pandas, shapely, etc). Expect kinks and continued growth!
# 
# **When should you use GeoPandas?**
# - For exploratory data analysis, including in Jupyter notebooks.
# - For highly compact and readable code. Which in turn improves reproducibility.
# - If you're comfortable with Pandas, R dataframes, or tabular/relational approaches.
# 
# **When it may not be the best tool?**
# - If you need high performance (though I'm not completely sure -- it uses a nice `rtree` index).
# - For polished map creation and multi-layer, interactive visualization; if you're comfortable with GIS, use a desktop GIS like QGIS! You can generate intermediate GIS files and plots with GeoPandas, then shift over to QGIS. Or refine the plots in Python with matplotlib or additional packages.

# ## 2. Set up packages and data file path
# We'll use these throughout the rest of the tutorial.

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')

import os

import matplotlib.pyplot as plt
# The two statemens below are used mainly to set up a plotting
# default style that's better than the default from matplotlib
import seaborn as sns
plt.style.use('bmh')

from shapely.geometry import Point
import pandas as pd
import geopandas as gpd
from geopandas import GeoSeries, GeoDataFrame

data_pth = "../data"


# ## 3. GeoSeries: The geometry building block

# Like a Pandas `Series`, a `GeoSeries` is the building block for the more broadly useful and powerful `GeoDataFrame` that we'll focus on in this tutorial. Here we'll take a bit of time to examine a `GeoSeries`.

# A `GeoSeries` is made up of an index and a GeoPandas `geometry` data type. This data type is a [shapely.geometry object](http://toblerity.org/shapely/manual.html#geometric-objects), and therefore inherits their attributes and methods such as `area`, `bounds`, `distance`, etc.
# 
# GeoPandas has six classes of **geometric objects**, corresponding to the three basic single-entity geometric types and their associated homogeneous collections of multiples entities:
# - **Single entity (core, basic types):**
#   - Point
#   - Line (*formally known as a LineString*)
#   - Polygon
# - **Homogeneous entity collections:**
#   - Multi-Point
#   - Multi-Line (*MultiLineString*)
#   - Multi-Polygon
# 
# A GeoSeries is then a list of geometry objects and their associated index values.

# **_NOTE/WATCH:_**   
# Entries (rows) in a GeoSeries can store different geometry types; GeoPandas does not constrain the geometry column to be of the same geometry type. This can lead to unexpected problems if you're not careful! Specially if you're used to thinking of a GIS file format like shape files, which store a single geometry type. Also beware that certain export operations (say, to shape files ...) will fail if the list of geometry objects is heterogeneous.

# But enough theory! Let's get our hands dirty (so to speak) with code. We'll start by illustrating how GeoSeries are constructured.

# ### Create a `GeoSeries` from a list of `shapely Point` objects constructed directly from `WKT` text (though you will rarely need this raw approach)

# In[2]:


from shapely.wkt import loads

GeoSeries([loads('POINT(1 2)'), loads('POINT(1.5 2.5)'), loads('POINT(2 3)')])


# ### Create a `GeoSeries` from a list of `shapely Point` objects
# Then enhance it with a crs and plot it.

# In[3]:


gs = GeoSeries([Point(-120, 45), Point(-121.2, 46), Point(-122.9, 47.5)])
gs


# In[4]:


type(gs), len(gs)


# A GeoSeries (and a GeoDataframe) can store a CRS implicitly associated with the geometry column. This is useful as essential spatial metadata and for transformation (reprojection) to another CRS.

# In[5]:


gs.crs = {'init': 'epsg:4326'}


# The `plot` method accepts standard `matplotlib.pyplot` style options, and can be tweaked like any other `matplotlib` figure.

# In[6]:


gs.plot(marker='*', color='red', markersize=12, figsize=(4, 4))
plt.xlim([-123, -119.8])
plt.ylim([44.8, 47.7]);


# **Let's get a bit fancier, as a stepping stone to GeoDataFrames.** First, we'll define a simple dictionary of lists, that we'll use again later.

# In[7]:


data = {'name': ['a', 'b', 'c'],
        'lat': [45, 46, 47.5],
        'lon': [-120, -121.2, -122.9]}


# Note this convenient, compact approach to create a list of `Point` shapely objects out of X & Y coordinate lists:

# In[8]:


geometry = [Point(xy) for xy in zip(data['lon'], data['lat'])]
geometry


# We'll wrap up by creating a GeoSeries where we explicitly define the index values.

# In[9]:


gs = GeoSeries(geometry, index=data['name'])
gs


# ## 4. GeoDataFrames: The real power tool.

# **_NOTE/HIGHLIGHT:_**   
# - It's worth noting that a GeoDataFrame can be described as a *Feature Collection*, where each row is a *Feature*, a *geometry* column is defined (thought the name of the column doesn't have to be "geometry"), and the attribute *Properties* are simply the other columns (the Pandas DataFrame part, if you will).
# - More than one column can store geometry objects! We won't explore this capability in this tutorial.

# ### Start with a simple, manually constructed illustration

# We'll build on the GeoSeries examples. Let's reuse the `data` dictionary we defined earlier, this time to create a DataFrame.

# In[10]:


df = pd.DataFrame(data)
df


# Now we use the DataFrame and the "list-of-shapely-Point-objects" approach to create a GeoDataFrame. Note the use of two GeoDataFrame attribute columns, which are just two simple Pandas Series.

# In[11]:


geometry = [Point(xy) for xy in zip(df['lon'], df['lat'])]
gdf = GeoDataFrame(df, geometry=geometry)


# There's nothing new to visualize, but this time we're using the `plot` method from a GeoDataFrame, *not* from a GeoSeries. They're not exactly the same thing under the hood.

# In[12]:


gdf.plot(marker='*', color='green', markersize=6, figsize=(2, 2));


# ### FINALLY, we get to work with real data! Load and examine the simple "oceans" shape file

# `gpd.read_file` is the workhorse for reading GIS files. It leverages the [fiona](http://toblerity.org/fiona/README.html) package.

# In[13]:


oceans = gpd.read_file(os.path.join(data_pth, "oceans.shp"))


# In[14]:


oceans.head()


# The `crs` was read from the shape file's `prj` file:

# In[15]:


oceans.crs


# Now we finally plot a real map (or blobs, depending on your aesthetics), from a dataset that's global and stored in "geographic" (latitude & longitude) coordinates. It'snot *quite* the actual ocean shapes defined by coastal boundaries, but bear with me.

# In[16]:


oceans.plot();


# `oceans.shp` stores both `Polygon` and `Multi-Polygon` geometry types (but a `Polygon` may be viewed as a `Multi-Polygon` with 1 member). We can get at the geometry types and other geometry properties easily.

# In[17]:


oceans.geom_type


# In[18]:


# Beware that these area calculations are in degrees, which is fairly useless
oceans.geometry.area


# In[19]:


oceans.geometry.bounds


# The `envelope` method returns the bounding box for each polygon. This could be used to create a new spatial column or GeoSeries; directly for plotting; etc.

# In[20]:


oceans.envelope.plot();


# Does it seem weird that some envelope bounding boxes, such as the North Pacific Ocean, span all longitudes? That's because they're Multi-Polygons with edges at the ends of the -180 and +180 degree coordinate range.

# In[21]:


oceans[oceans['Oceans'] == 'North Pacific Ocean'].plot();


# ### Load "Natural Earth" countries dataset, bundled with GeoPandas
# "[Natural Earth](http://www.naturalearthdata.com) is a public domain map dataset available at 1:10m, 1:50m, and 1:110 million scales. Featuring tightly integrated vector and raster data, with Natural Earth you can make a variety of visually pleasing, well-crafted maps with cartography or GIS software." It (a subset?) comes bundled with GeoPandas and is accessible from the `gpd.datasets` module. We'll use it as a helpful global base layer map.

# In[22]:


world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.head(2)


# Its CRS is also EPSG:4326:

# In[23]:


world.crs


# In[24]:


world.plot();


# ### Map plot overlays: Plotting multiple spatial layers

# Here's a compact, quick way of using GeoDataFrame plot method to overlay two GeoDataFrame, while style customizing the styles for each layer.

# In[25]:


world.plot(ax=oceans.plot(cmap='Set2', alpha=1), alpha=1);


# **_NOTE/WATCH:_**   
# The oceans polygon boundaries are coming through the world dataset! I've looked into this and can't find a solution. I think it's a bug in the GeoPandas underlying plotting library or plotting approach.

# We can also compose the plot using conventional `matplotlib` steps and options that give us more control.

# In[26]:


f, ax = plt.subplots(1, figsize=(10, 5))
# Other nice categorical color maps (cmap) include 'Set2' and 'Set3'
oceans.plot(cmap='Paired', alpha=1, linewidth=0.2, ax=ax)
world.plot(alpha=1, ax=ax)
ax.set_ylim([-100, 100])
ax.set_title('Countries and Ocean Basins')
plt.axis('equal');


# ## 5. Extras: Reading from other data source types; fancier plotting
# - Read another shape file and explore filtering and plotting. In the [original UW GeoHack Week tutorial](https://geohackweek.github.io/vector/04-geopandas-intro/), this exercise involved reading geospatial data from a remote PostGIS source.
# - Read from a remote OGC WFS service.

# ### Read from Shapefile with GeoDataFrame

# In[27]:


seas = GeoDataFrame.from_file(os.path.join(data_pth, "World_Seas.shp"))


# Let's take a look at the GeoDataFrame.

# In[28]:


seas.head()


# ### More advanced plotting and data filtering

# Color the layer based on one column that aggregates individual polygons; using a categorical map, as before, but explicitly selecting the column (`column='oceans'`) and categorical mapping (`categorical=True`); dispplaying an auto-generated legend; while displaying all polygon boundaries. "oceans" (ocean basins, actually) contain one or more 'seas'.

# In[29]:


seas.plot(column='oceans', categorical=True, legend=True, figsize=(14,6));


# **_NOTE/COOL:_**  
# See http://darribas.org/gds15/content/labs/lab_04.html for great examples of lots of other cool GeoPandas map plotting tips.

# Combine what we've learned. A map overlay, using `world` as a background layer, and filtering `seas` based on an attribute value (from `oceans` column) and an auto-derived GeoPandas geometry attribute (`area`). **`world` is in gray scale, while the filtered `seas` is in color.**

# In[30]:


seas_na_arealt1000 = seas[(seas['oceans'] == 'North Atlantic Ocean') 
                          & (seas.geometry.area < 1000)]


# In[31]:


seas_na_arealt1000.plot(ax=world.plot(alpha=0.1), cmap='Paired')

# Use the bounds geometry attribute to set a nice
# geographical extent for the plot, based on the filtered GDF
bounds = seas_na_arealt1000.geometry.bounds

plt.xlim([bounds.minx.min()-5, bounds.maxx.max()+5])
plt.ylim([bounds.miny.min()-5, bounds.maxy.max()+5]);


# ### Save the filtered seas GeoDataFrame to a shape file
# The `to_file` method uses the [fiona](http://toblerity.org/fiona/README.html) package to write to a GIS file. The default `driver` for output file format is 'ESRI Shapefile', but many others are available because `fiona` leverages [GDAL/OGR](http://www.gdal.org).

# In[32]:


seas_na_arealt1000.to_file(os.path.join(data_pth, "seas_na_arealt1000.shp"))


# ### Read from OGC WFS GeoJSON response into a GeoDataFrame
# Use an [Open Geospatial Consortium](http://www.opengeospatial.org) (OGC) [Web Feature Service](https://en.wikipedia.org/wiki/Web_Feature_Service) (WFS) request to obtain geospatial data from a remote source. OGC WFS is an open geospatial standard.

# We won't go into all details here about what's going on. Suffice it to say that we issue an OGC WFS request for all features from the layer named "oa:goainv" found in a [GeoServer](http://geoserver.org) instance from [NANOOS](http://nanoos.org), requesting the response in `GeoJSON` format. Then we "load" it into a `geojson` feature object (basically a dictionary) using the `geojson` package.
# 
# The "oa:goainv" layer is a global dataset of monitoring sites and cruises where data relevant to ocean acidification is collected. It's a work in progress from the [Global Ocean Acidification Observation Network (GOA-ON)](http://www.goa-on.org); for additional information see the [GOA-ON Data Portal](http://portal.goa-on.org).

# In[33]:


import requests
import geojson

wfs_url = "http://data.nanoos.org/geoserver/ows"
params = dict(service='WFS', version='1.0.0', request='GetFeature',
              typeName='oa:goaoninv', outputFormat='json')

r = requests.get(wfs_url, params=params)
wfs_geo = geojson.loads(r.content)


# Let's examine the general characteristics of this GeoJSON object. We'll take advantage of the `__geo_interface__` interface we discussed earlier.

# In[34]:


print(type(wfs_geo))
print(wfs_geo.keys())
print(len(wfs_geo.__geo_interface__['features']))


# Now we use the `from_features` constructor method to create a GeoDataFrame, passing to it the `features` from the `__geo_interface__` method.

# In[35]:


wfs_gdf = GeoDataFrame.from_features(wfs_geo.__geo_interface__['features'])


# Display the values for the last feature, as an example.

# In[36]:


wfs_gdf.iloc[-1]


# Finally, a simple map overlay plot.

# In[37]:


wfs_gdf.plot(ax=world.plot(alpha=1), figsize=(10, 6),
             marker='o', color='red', markersize=4);


# In[ ]: