In [1]:

```
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from IPython.display import Image
import pandas as pd
import cartopy.crs as ccrs
from cartopy.feature import NaturalEarthFeature, LAND, COASTLINE, OCEAN, LAKES, BORDERS
from cartopy import feature as cfeature
import warnings
# this will allow us to use the code in peace :)
warnings.filterwarnings("ignore")
```

- Learn about geoplot and geopandas
- Learn a bit about coordinate systems (UTM versus WGS84, as examples)
- Learn something about Hawaiian volcanism

By now we have learned a great deal about plotting and map making. In this lecture we will pull together a number of concepts and add some new tools (geoplot and geopandas). **Geopandas** is a software package built on top of **Pandas** that is made for dealing with spatial data. **Geoplot** is built on top of **matplotlib** and has a number of mapping tools that we can use. **Geoplot** is to **cartopy** as **seaborn** is to **matplotlib**.

If you haven't installed these, you can use **conda** on your command line or within this notebook (on macs).

In [2]:

```
#!conda install geoplot -c conda-forge
# and
#!conda install geopandas
```

In [3]:

```
import geopandas as gpd
import geoplot as gplt
```

By now you should be familiar with the rudiments of **cartopy**. We will now learn how to make more complicated plots using the volcanoes active over the last 6 million years as an example. We downloaded the data from this website: https://www.navdat.org/NavdatSearch/Search.cfm by setting the Age filter to 6 (Ma) and exporting the data as an excel spreadsheet which in turn got saved as a tab delimited text file in *Datafiles/navdat.txt*)

The first new thing we need to grapple with is the **geopandas** overlay on **Pandas**. This is a very powerful, but also complicated software package. The first thing we will use it for is wrangling our volcano dataset into a form that we can plot with **geoplot**.

A **Geopandas DataFrame** differs from a **Pandas DataFrame** in that it has a **geometry** field which can be points, polygons or other spatial data types.

We'll approach this as follows:

- we need to first read in the data set with
**Pandas** - make a map as we did in Lecture 18
- turn the
**Pandas DataFrame**into a**geopandas**one by putting in the*geometry*information required by**geoplot**. - Use
**geoplot**to do some new stuff.

In [4]:

```
#Read in the data file
df=pd.read_csv('Datasets/navdat.txt',sep='\t')
df.head()
```

Out[4]:

Make the map (cribbing from Lecture 18).

In [5]:

```
proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33)
ax = plt.axes(projection=proj)
ax.set_extent([-130, -70, 15, 52], crs=ccrs.PlateCarree())
gl=ax.gridlines(ylocs=np.arange(0,90,15.),xlocs=np.arange(-180.,180,15.),\
linewidth=2, linestyle="dotted")
ax.coastlines();
ax.add_feature(BORDERS,linestyle='-',linewidth=2)
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
edgecolor='purple',
facecolor='none',
linestyle='dotted')
ax.add_feature(states_provinces);
```

Now for the new stuff. We nedd to use the **.GeoDataFrame** method of **geopandas** to put in the geometry field required by **geoplots**.

In [6]:

```
gdf=gpd.GeoDataFrame(df,geometry=gpd.points_from_xy(df.LONGITUDE,df.LATITUDE))
gdf.head()
```

Out[6]:

You can see that **Geopandas** put in a new column, called *geometry* with the lat/lon as a tuple with the geometry type of *POINT*.

Now we can just plot the points with **pointplot** from **geoplot** to make the heat map! As with **seaborn** we can use one of the columns to determine color (**hue**), so in this example, we color the points by age:

In [7]:

```
# same as before:
proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33)
ax = plt.axes(projection=proj)
ax.set_extent([-130, -70, 15, 52], crs=ccrs.PlateCarree())
ax.coastlines();
ax.add_feature(BORDERS,linestyle='-',linewidth=2)
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
edgecolor='purple',
facecolor='none',
linestyle='dotted')
ax.add_feature(states_provinces);# The NEW stuff!
gplt.pointplot(gdf,ax=ax,hue='CALCULATED AGE',legend=True);
```

There are many more plots that we can make with **geoplots**. See this website for all the options: https://residentmario.github.io/geoplot/plot_references/plot_reference.html

Let's go back to the data set you look at in Lectures 15 and 18 and make a better heatmap that with **plt.hist2d()**.

In [8]:

```
#Read in the data file
df=pd.read_csv('Datasets/WUS_navdat.txt',sep='\t')
gdf=gpd.GeoDataFrame(df,geometry=gpd.points_from_xy(df.LONGITUDE,df.LATITUDE))
```

In [9]:

```
df.head()
```

Out[9]:

In [10]:

```
# same as before:
proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33)
ax = plt.axes(projection=proj)
ax.add_feature(BORDERS,linestyle='-',linewidth=2)
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
edgecolor='purple',
facecolor='none',
linestyle='dotted')
ax.add_feature(states_provinces);# The NEW stuff!
ax.coastlines()
# but with the KDE plot.
gplt.kdeplot(gdf,ax=ax);
```

That made a contour plot for the locations of volcanoes in the Western US. What we were after was a heatmap in which color indicates density of volcanoes. We can do that with the **shade** argument and we can also specify a colormap with the familiar **cmap** argument. Maybe **cmap=inferno** is appropriate?

In [11]:

```
# same as before:
proj = ccrs.LambertConformal(central_longitude=260, central_latitude=33)
ax = plt.axes(projection=proj)
ax.add_feature(BORDERS,linestyle='-',linewidth=2)
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
edgecolor='purple',
facecolor='none',
linestyle='dotted')
ax.add_feature(states_provinces);# The NEW stuff!
ax.coastlines()
# but with the KDE plot.
gplt.kdeplot(gdf,ax=ax,shade=True,cmap='inferno');
```

**geoplot** also allows us to access open source datasets like OpenStreetMap using **gplt.webmap()**.

In [12]:

```
help(gplt.webmap)
```

In [13]:

```
gplt.webmap(gdf,figsize=(10,10),provider='OSM_A');
```

The big advantage of **geopandas** is the ability to plot polygons. Polygons outline geological features such as formations, age units, etc. Coastlines are really polygons too. In the rest of this lecture, we will explore making maps with some geologically interesting polygons - historical lava flows in Hawaii.

We will proceed as follows:

- get the coastline of the Big Islanad of Hawaii into a
**geopandas.DataFrame**.- make a quick plot of Hawaii as a lightgreen polygon in a sea of lightblue.

- find the elevation data for the Big Island of Hawaii.
- make a heatmap of the elevation data
- make a contour plot of the elevation data (a topographic map).

- find the lava flow outlines from https://pubs.usgs.gov/of/2007/1089/.
- these data are in Datafiles/BigIslandData.shp and are in Universal Transverse Mercator and we really want them in lat/lon (in fact, on the WGS84 ellipsoid). Converting between different coordinate systems just means we to know the transforms. Here is a big list of them: http://spatialreference.org.
- plot the converted polygons on the topographic map.

Read in the coastlines for the State of Hawaii, filter for the Island Hawaii (it's the big one with the active lava flows) and plot the outline with the **plot** method on **geopandas** **DataFrames**.

In [14]:

```
# read in the data as a geopandas dataframe
hawaii=gpd.read_file('Datasets/Coastline.shp')
# take a quick look at the format:
hawaii.head()
```

Out[14]:

There are data from all the Hawaiian islands, so we should filter for the Big Island, which is actually named "Hawaii" too.

In [15]:

```
#Select only the data from island called "Hawaii".
bigIsland=hawaii[(hawaii.isle=='Hawaii')]
fig, ax = plt.subplots(figsize=(8,6)) # make fig and ax objects
#Set the ocean to blue.
ax.set_facecolor('lightblue')
# plot the polygons from the geometry column
bigIsland.plot(ax=ax,edgecolor='black',linewidth=1,facecolor='lightgreen');
```

Ok, that was pretty easy. Now let's put on the topography. First, let's make a heatmap of the elevation data. In Lecture 19 we learned how to make **meshgrids** to use **plt.contour( )** and **plt.contourf( )**. There is a cool way to make the same plot using **tricontour( )** and **tricontourf** from the **matplotlib.tri** module.

To make life easier, we can use the **matplotlib.tri.tricontourf()** function. This takes unstructured x,y,z data and constructs all the meshes for us.

But we have to import the function first, as this is a new one. It has a number of functions which you can explore (try **help(tri)** after importing the module.

In [16]:

```
import matplotlib.tri as tri
```

In [17]:

```
fig, ax = plt.subplots(figsize=(8,6))
# read in the elevation data as a Pandas DatFrame.
bigIslandElev=pd.read_csv('Datasets/BigIslandElev.csv')
# plot longitude as x, latitude as y and elevation as z:
x=bigIslandElev.lon.values
y=bigIslandElev.lat.values
z=bigIslandElev.elevation.values
ax.tricontourf(x,y,z);
```

Well that was relatively painless. Now let's plot the same data, but as a contour plot.

In [18]:

```
fig, ax = plt.subplots(figsize=(8,6))
#Setting the ocean to blue.
ax.set_facecolor('lightblue')
bigIsland.plot(ax=ax,edgecolor='black',linewidth=1,facecolor='lightgreen')
ax.tricontour(x,y,z);
```

and a few whistles and bells:

In [19]:

```
fig, ax = plt.subplots(figsize=(8,6))
#Setting the ocean to blue.
ax.set_facecolor('lightblue')
bigIsland.plot(ax=ax,edgecolor='black',linewidth=1,facecolor='lightgreen')
ax.tricontour(x,y,z,colors='k',linewidths=1,levels=np.arange(-5000,5000,1000));
```

Now we can read in the flow polygons and convert them from UTM to WGS84.

But first what are 'UTM' and 'WGS84'?

There are many coordinate systems, but the two used with GPS are geodetic latitude/longitude/elevation and UTM (Universal Transverse Mercator).

Lat/Lon/Elevation: The prime meridian and the equator are the reference planes used to define latitude and longitude. The geodetic latitude of a point is the angle from the equatorial plane to the vertical direction of a line normal to the reference ellipsoid. The geodetic longitude of a point is the angle between a reference plane and a plane passing through the point, both planes being perpendicular to the equatorial plane. The geodetic elevation is at a point is the distance from the reference ellipsoid to the point in a direction normal to the ellipsoid. We will be using the reference ellipsoid known as WGS84 (although there are many so be careful!)

The Universal Transverse Mercator (UTM) coordinates define map locations (2D) within zones. Zone numbers designate 6 degree longitudinal strips. The letters go from A (south pole) to Z (north pole). So Each zone has a central meridian (CM). Northings/eastings are measured from the equator and the CM. San Diego is in zone S11, for example.

It turns out that with **Geopandas**, conversion is pretty straightforward, once you know the details of what coordinate system you are in. For this exercise, the data are from Hawaii and are in zone 4 based on the NAD83 ellipsoid in units of meters. We want them in longitude/latitude on ellipsoid WGS84.

In [20]:

```
Image('Figures/coordinateSystems.png')
```

Out[20]:

Lets set our conversion strings, read in the data and first take a look at the geometry field in UTM.

In [21]:

```
#Set Conversion Strings
UTM4='+proj=utm +zone=4 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'
WGS84='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
#Get Big Island Flows and convert to WGS84
bigIslandFlows=gpd.read_file('Datasets/BigIslandData.shp',crs=UTM4)
# in UTM4
bigIslandFlows.head()
```

Out[21]:

Looking at the polygon for the first flow (ID=5522), we can see how the POLYGON column has points like (862851.7650805543 2185118.765394682,...). These are the Easting and Northing of a point on the flow with ID 5522. How do I know which is easting and which is northing??? (HINT: think about the definitions).

In [22]:

```
bigIslandFlows.loc[bigIslandFlows.ID==5522]['geometry']
```

Out[22]:

We can convert between coordinate systems when we read the file in like this:

In [23]:

```
bigIslandFlows=gpd.read_file('Datasets/BigIslandData.shp',crs=UTM4).to_crs(WGS84)
bigIslandFlows.head()
```

Out[23]:

And now notice how the polygon is in longitude/latitude:

In [24]:

```
bigIslandFlows.loc[bigIslandFlows.ID==5522]['geometry']
```

Out[24]:

And now for the plot!

In [26]:

```
# from before
fig, ax = plt.subplots(figsize=(8,6))
#Setting the ocean to blue.
ax.set_facecolor('lightblue')
bigIsland.plot(ax=ax,edgecolor='black',linewidth=1,facecolor='lightgreen')
ax.tricontour(x,y,z,colors='k',linewidths=1,levels=np.arange(-5000,5000,500))
# now the flow polygons:
bigIslandFlows.plot(ax=ax,column='Ages',cmap='magma',legend=True,legend_kwds={'label':'Date of Eruption (A.D.)'})
plt.xlim(-156.2,-154.8)
plt.savefig('hawaii.jpg');
```