CODEWITHSUNDEEP
pythonpandasmatplotlibplotgeopandas
BenCaldwell
BenCaldwell

Reputation: 443

How to overlay a quiver plot on a geodataframe plot in python

I would like to overlay a quiver plot of wind direction on a basemap in my jupyter notebook. I have a pandas dataframe that includes columns: | Latitude | Longitude | True Wind Inferred |

I have already used geopandas to create a geodataframe and plot gps track data on an osm basemap using contextily (code below). I have also been able to bin the latitude and longitude to get average True Wind Inferred (wind direction) for a "box" on the map. However, I haven't found any examples on how to plot a quiver plot of the binned True Wind Inferred in the boxes. I have only plotted as a scatter plot so far but the colour map does not visualize directional data well.

Imports:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# plot inline graphics
%pylab inline
import os.path
from shapely.geometry import Point
import geopandas as gpd
import contextily as ctx

Sample dataframe:

df[['Latitude', 'Longitude', 'True Wind Inferred', 'coords']].head()

    Latitude    Longitude   True Wind Inferred  coords
0   -31.991899  115.848825  173.835559  POINT (115.848825 -31.991899)
1   -31.992036  115.848873  182.620880  POINT (115.848873 -31.992036)
2   -31.992181  115.848895  192.140276  POINT (115.848895 -31.992181)
3   -31.992308  115.848832  206.655730  POINT (115.848832 -31.992308)
4   -31.992430  115.848784  218.656646  POINT (115.848784 -31.99243)

Binning the dataframe:

step = 0.005
to_bin = lambda x: np.floor(x / step) * step
dfLocBin['latbin'] = df['Latitude'].map(to_bin)
dfLocBin['lonbin'] = df['Latitude'].map(to_bin)
dfLocBin = df.groupby(['latbin', 'lonbin'])[['True Wind Inferred']].mean()
dfLocBin.reset_index(inplace=True)
dfLocBin['coords'] = list(zip(dfLocBin['lonbin'], dfLocBin['latbin']))
dfLocBin['coords'] = dfLocBin['coords'].apply(Point)
dfLocBin.head()

    latbin  lonbin  True Wind Inferred  coords
0   -32.015 115.790 223.149075  POINT (115.79 -32.015)
1   -32.015 115.795 222.242870  POINT (115.795 -32.015)
2   -32.015 115.800 223.710092  POINT (115.8 -32.015)
3   -32.015 115.805 225.887096  POINT (115.805 -32.015)
4   -32.015 115.810 225.298059  POINT (115.81 -32.015)

And plotting:

def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

geo_df = gpd.GeoDataFrame(
    dfLocBin, crs  ={'init': 'epsg:4326'},
    geometry = dfLocBin['coords']
).to_crs(epsg=3857)

ax = geo_df.plot(
    figsize= (20, 20),
    alpha  = 1,
    c=dfLocBin['True Wind Inferred']
)

add_basemap(ax, zoom=15, url=ctx.tile_providers.ST_TONER)
ax.set_axis_off()
plt.title('Binned True Wind Direction')
plt.show()

scatter plot

I would like to change the type of plot from a scatter with colours to a quiver plot with arrows representing the compass direction of the wind.

Upvotes: 1

Views: 2088

Answers (1)

BenCaldwell
BenCaldwell

Reputation: 443

I worked it out. The X,Y for the quiver need to come from the geometry of the geodataframe to plot properly on the same axis. The geodataframe columns look like:

geo_df.head()

latbin  lonbin  True Wind Inferred  coords  geometry
0   -32.014 115.798 220.492453  POINT (115.798 -32.014) POINT (12890574.39487949 -3765148.48502445)
1   -32.014 115.800 225.718756  POINT (115.8 -32.014)   POINT (12890797.03386108 -3765148.48502445)

Working Code:

# bin the coordinates and plot a vector field
step = 0.002
to_bin = lambda x: np.floor(x / step) * step
df['latbin'] = df['Latitude'].map(to_bin)
df['lonbin'] = df['Longitude'].map(to_bin)
dfLocBin = df.groupby(['latbin', 'lonbin'])[['True Wind Inferred']].mean()
dfLocBin.reset_index(inplace=True)
dfLocBin['coords'] = list(zip(dfLocBin['lonbin'], dfLocBin['latbin']))
dfLocBin['coords'] = dfLocBin['coords'].apply(Point)

# ... turn them into geodataframe, and convert our
# epsg into 3857, since web map tiles are typically
# provided as such.
geo_df = gpd.GeoDataFrame(
    dfLocBin, crs  ={'init': 'epsg:4326'},
    geometry = dfLocBin['coords']
).to_crs(epsg=3857)

# ... and make the plot
ax = geo_df.plot(
    figsize= (20, 20),
    alpha  = 1
)

geo_df['X'] = geo_df['geometry'].x
geo_df['Y'] = geo_df['geometry'].y

geo_df['U'] = np.cos(np.radians(geo_df['True Wind Inferred']))
geo_df['V'] = np.sin(np.radians(geo_df['True Wind Inferred']))

ax.quiver(geo_df['X'], 
          geo_df['Y'], 
          geo_df['U'], 
          geo_df['V'],
         color='deepskyblue')


add_basemap(ax, zoom=15, url=ctx.tile_providers.ST_TONER)

ax.set_axis_off()
plt.title('Binned True Wind Direction')
plt.show()

quiver on map

Upvotes: 3

Related Questions