Introductie tot geospatiale data in Python

In deze tutorial maak je kennis met twee populaire pakketten om met geospatiale data te werken: geopandas en Shapely. Daarna pas je deze pakketten toe om de geospatiale data met Python in te lezen en het spoor van orkaan Florence te plotten van 30 augustus tot en met 18 september.

Wat is geospatiale data?

Ruimtelijke data, geospatiale data, GIS-data of geodata zijn benamingen voor numerieke data die de geografische locatie van een fysiek object identificeert, zoals een gebouw, een straat, een dorp, een stad, een land, enzovoort, volgens een geografisch coördinatenstelsel. Uit ruimtelijke data kun je niet alleen de locatie afleiden, maar ook de lengte, omvang, oppervlakte of vorm van een object. Een voorbeeld van ruimtelijke data die je kunt verkrijgen zijn coördinaten van een object, zoals breedtegraad, lengtegraad en hoogte. Geografische informatiesystemen (GIS) of andere gespecialiseerde softwareapplicaties kunnen worden gebruikt om geospatiale data te benaderen, te visualiseren, te bewerken en te analyseren.

Waarom geospatiale data?

Geospatiale data heeft talloze toepassingen in het dagelijks leven. Een voorbeeld is de Kaarten-app die je gebruikt om van de ene naar de andere locatie te navigeren. Een ander voorbeeld dat je misschien elke dag ziet, is op het weerkanaal.

Ja, geospatiale data wordt gebruikt om posities weer te geven van iets ten opzichte van de omgeving: jouw huis op de stadskaart, de orkaan op de wereldkaart, enzovoort. In deze tutorial kijk je naar de verwoestende orkaan Florence en volg je zijn locatie.

Pakketvereisten

Niet overslaan. Dit is een belangrijk onderdeel. Als je niet zeker weet of je aan de vereisten voldoet, check het dan nog eens. Allereerst moet je alle onderstaande pakketten geïnstalleerd hebben.

• Pandas: biedt datastructuren en tools voor data-analyse
• Numpy: een fundamenteel pakket voor wetenschappelijk rekenen met Python
• SciPy:(uitspraak “Saj Pai”) is een Python-gebaseerd ecosysteem van open-sourcesoftware voor wiskunde, wetenschap en engineering
• RTree: een ctypes Python-wrapper van libspatialindex die geavanceerde ruimtelijke indexeermogelijkheden biedt
• GDAL: translatiebibliotheek voor raster- en vectorformaten voor geospatiale data
• Fiona: Fiona leest en schrijft ruimtelijke databestanden
• Shapely: Geometrische objecten, predicaten en bewerkingen
• GeoPandas: breidt de datatypes van pandas uit om ruimtelijke operaties op geometrische types mogelijk te maken.
• PySAL: een bibliotheek met functies voor ruimtelijke analyse in Python ter ondersteuning van de ontwikkeling van high-level applicaties.
• Matplotlib: Python 2D-plotbibliotheek
• Missingno: Visualisatiemodule voor missende data in Python

Installeer alle pakketten in bovenstaande volgorde om te zorgen dat alles werkt. Sommige pakketten zijn vereisten voor andere; bijvoorbeeld: voor GeoPandas is Shapely vereist, en voor Shapely moeten RTree, GDAL en Fiona geïnstalleerd zijn. De eenvoudigste manier om een pakket te installeren is: pip install PACKAGE_NAME

Als je Windows gebruikt en het pakket zo niet kunt installeren, bekijk dan deze website om de betreffende pakketten te downloaden en voer pip install PATH_TO_PACKAGE uit.

Over de data

Zoals je wellicht weet, is de angstaanjagende orkaan Florence net over een deel van de oostkust van de Verenigde Staten getrokken, met een geschatte schade van 17 miljard USD. Deze tutorial helpt je uit te zoeken waar hij vandaan kwam, wanneer en waar hij krachtiger werd, en om deze natuurramp beter te begrijpen en te analyseren in Python. Er zijn veel websites met informatie over deze orkaan. Zo biedt deze website gegevens over verschillende stormen en orkanen in de VS van 1902 tot 2018, dus er is een berg data om later mee te werken. Voor deze tutorial gebruik je alleen de data van orkaan Florence

Je gebruikt ook de geospatiale data van de kaart van de VS van internet. De blogpost van Eric Celeste bevat allerlei grensbestanden voor Amerikaanse counties en staten. De data die in deze tutorial wordt gebruikt is het GeoJSON-bestand US States, 5m.

 # Load all importance packages
import geopandas
import numpy as np
import pandas as pd
from shapely.geometry import Point

import missingno as msn

import seaborn as sns
import matplotlib.pyplot as plt

% matplotlib inline

Laten we eerst kijken naar de eerste geospatiale dataframe: US States Geodata

# Getting to know GEOJSON file:
country = geopandas.read_file("data/gz_2010_us_040_00_5m.json")
country.head()

Door het type van de dataframe te controleren die je zojuist hebt ingeladen, zie je dat het een GeoDataFrame is, met alle reguliere eigenschappen van een Pandas DataFrame.

type(country)

geopandas.geodataframe.GeoDataFrame

Het datatype van de kolom met coördinaten controleren: dat is een GeoSeries.

type(country.geometry)

geopandas.geoseries.GeoSeries

Elke waarde in de GeoSeries is een Shapely-object. Dat kan zijn:

• Point
• Line
• Polygon
• MultiPolygon

Elk object kan gebruikt worden voor een ander soort fysiek object, zoals: Point voor een gebouw, Line voor een straat, Polygon voor een stad en MultiPolygon voor een land met meerdere steden erin. Voor meer informatie over elk geometrisch object, lees dit artikel: https://shapely.readthedocs.io/en/stable/manual.html#geometric-objects

type(country.geometry[0])

shapely.geometry.multipolygon.MultiPolygon

Net als een Pandas DataFrame heeft een GeoDataFrame ook het attribuut plot, dat de geometrie in de dataframe gebruikt om een kaart te plotten:

country.plot()

<matplotlib.axes._subplots.AxesSubplot at 0x1cfe68c1358>

Zoals je ziet is de kaart van de VS relatief klein ten opzichte van het frame. Dat komt doordat de informatie Alaska, Hawaï en Puerto Rico bevat, die verspreid liggen. Voor deze tutorial kun je Alaska en Hawaï uitsluiten, aangezien de orkaan niet in de buurt van deze staten kwam. Je kunt ook de figuurgrootte en kleur toevoegen om je eigen plot aan te passen:

# Exclude Alaska and Hawaii for now
country[country['NAME'].isin(['Alaska','Hawaii']) == False].plot(figsize=(30,20), color='#3B3C6E');

Niet zo moeilijk, toch! Nu je de kaart van de VS hebt, laden we de orkaandata in:

florence = pd.read_csv('data/florence.csv')
florence.head()

Exploratory Data Analysis

Dit is altijd het eerste wat je doet wanneer je een dataset laadt:

• Informatie en datatypes controleren
• Eventuele missende waarden
• Statistische gegevens

florence.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 105 entries, 0 to 104
Data columns (total 11 columns):
AdvisoryNumber    105 non-null object
Date              105 non-null object
Lat               105 non-null float64
Long              105 non-null float64
Wind              105 non-null int64
Pres              105 non-null int64
Movement          105 non-null object
Type              105 non-null object
Name              105 non-null object
Received          105 non-null object
Forecaster        104 non-null object
dtypes: float64(2), int64(2), object(7)
memory usage: 9.1+ KB

Missende waarden controleren met het pakket missingno. Dit is een handig pakket dat visualisatie gebruikt om missende data te tonen. Zoals je hieronder ziet, is er slechts één missende waarde in de kolom "Forecaster" die je voor deze tutorial niet nodig hebt. Die kun je dus negeren.

# Notice you can always adjust the color of the visualization
msn.bar(florence, color='darkolivegreen');

Kijk naar wat statistische informatie; sommige dingen kunnen erg nuttig zijn, zoals de gemiddelde windsnelheid, maximale en minimale windsnelheid van deze orkaan, enz.

# Statistical information
florence.describe()

Voor de meeste data moet je opschonen en alleen meenemen wat je nodig hebt. Hier heb je alleen de tijd, de coördinaten: breedte- en lengtegraad, windsnelheid, luchtdruk en naam nodig. Movement en Type zijn optioneel, de rest kan worden gedropt.

# dropping all unused features:
florence = florence.drop(['AdvisoryNumber', 'Forecaster', 'Received'], axis=1)
florence.head()

Normaal gesproken hoef je bij het plotten van de data op zichzelf geen extra aandacht te besteden aan de coördinaten. Maar als je wilt dat het eruitziet zoals op de kaart, is het belangrijk om de lengte- en breedtegraad te controleren. Hier is de lengtegraad westelijk; je moet een "-" voor het getal zetten om de data correct te plotten:

# Add "-" in front of the number to correctly plot the data:
florence['Long'] = 0 - florence['Long']
florence.head()

Vervolgens kun je breedte- en lengtegraad combineren om orkaancoördinaten te maken, die je daarna omzet in GeoPoint voor visualisatie.

# Combining Lattitude and Longitude to create hurricane coordinates:
florence['coordinates'] = florence[['Long', 'Lat']].values.tolist()
florence.head()

# Change the coordinates to a geoPoint
florence['coordinates'] = florence['coordinates'].apply(Point)
florence.head()

Het type van de florence-dataframe en de kolom coordinates van de florence-data controleren. Dat zijn respectievelijk een pandas DataFrame en een pandas Series.

type(florence)

pandas.core.frame.DataFrame

type(florence['coordinates'])

pandas.core.series.Series

Nadat we de data naar geospatiale data hebben omgezet, controleren we het type van de florence-dataframe en de kolom coordinates van de Florence-data opnieuw. Nu zijn het een GeoDataFrame en GeoSeries.

# Convert the count df to geodf
florence = geopandas.GeoDataFrame(florence, geometry='coordinates')
florence.head()

type(florence)

geopandas.geodataframe.GeoDataFrame

type(florence['coordinates'])

geopandas.geoseries.GeoSeries

Let op: ook al is het nu een GeoDataFrame en GeoSeries, het gedraagt zich nog steeds als een normale DataFrame en Series. Dat betekent dat je nog steeds kunt filteren, groupby kunt gebruiken voor de dataframe, of de min-, max- of gemiddelde waarden van de kolom kunt bepalen.

# Filtering from before the hurricane was named.
florence[florence['Name']=='Six']

# Groupping by name to see how many names it has in the data set:
florence.groupby('Name').Type.count()

Name
FLORENCE     6
Florence    85
SIX          4
Six         10
Name: Type, dtype: int64

De gemiddelde windsnelheid van orkaan Florence bepalen:

print("Mean wind speed of Hurricane Florence is {} mph and it can go up to {} mph maximum".format(round(florence.Wind.mean(),4),
                                                                                         florence.Wind.max()))

Mean wind speed of Hurricane Florence is 74.4286 mph and it can go up to 140 mph maximum

De gemiddelde windsnelheid van orkaan Florence is dus 74,43 mijl per uur (119,78 km per uur) en de maximumsnelheid is 140 mijl per uur (225,308 km per uur). Om je een idee te geven hoe beangstigend deze windsnelheid is: de website Beaufort Wind Scale, ontwikkeld door de Britse Royal Navy, laat zien hoe windeffecten eruitzien op water en land. Bij snelheden van 48 tot 55 mijl per uur kunnen bomen al afbreken of ontwortelen en ontstaat er "aanzienlijke structurele schade".

Daar wil je op dat moment niet zijn.

Visualisatie

Net als bij een pandas DataFrame heeft een GeoDataFrame ook het attribuut .plot. Dit attribuut gebruikt echter de coördinaten binnen het GeoDataFrame om te tekenen. Kijk maar:

florence.plot(figsize=(20,10));

Wat is er gebeurd? Je ziet alleen een stel punten zonder referentie. Is er iets mis?
Nee, alles is in orde. Omdat deze dataframe alleen coördinateninformatie (locatie) bevat van orkaan Florence op elk tijdstip, kunnen we alleen de positie op een blanco kaart plotten.

De volgende stap is dus de orkaanpositie op de kaart van de VS plotten om te zien waar hij kwam en hoe sterk hij op dat moment was. Hiervoor gebruik je de kaartcoördinaten van de VS (de data die we in het begin hebben geladen) als basis en plot je de positie van orkaan Florence erbovenop.

 # Plotting to see the hurricane overlay the US map:
fig, ax = plt.subplots(1, figsize=(30,20))
base = country[country['NAME'].isin(['Alaska','Hawaii']) == False].plot(ax=ax, color='#3B3C6E')

# plotting the hurricane position on top with red color to stand out:
florence.plot(ax=base, color='darkred', marker="*", markersize=10);

Ziet er goed uit! Nu maken we het af met meer details zoals:

• Een titel toevoegen
• De orkaanposities kleuren op basis van de windsnelheid om te zien hoe sterk de orkaan was toen hij elke stad bereikte.
• Assen verwijderen
• Legenda toevoegen
• Het resultaat opslaan als afbeeldingsbestand voor later gebruik

 fig, ax = plt.subplots(1, figsize=(20,20))
base = country[country['NAME'].isin(['Alaska','Hawaii']) == False].plot(ax=ax, color='#3B3C6E')
florence.plot(ax=base, column='Wind', marker="<", markersize=10, cmap='cool', label="Wind speed(mph)")
_ = ax.axis('off')
plt.legend()
ax.set_title("Hurricane Florence in US Map", fontsize=25)
plt.savefig('Hurricane_footage.png',bbox_inches='tight');

De orkaan was dus het sterkst toen hij voor de kust bij de oostkust lag. Toen hij het land naderde, begon de orkaan aan kracht in te boeten, maar met windsnelheden tussen 60 en 77 mijl per uur kan er nog steeds enorme schade ontstaan.

Conclusie

Goed gedaan! Je hebt een aantal essentiële stappen geleerd om met geospatiale data in Python te werken. Je hebt ook geleerd hoe je geospatiale data plot en de vorm, kleur en overlay van plots aanpast om een verhaal te vertellen. Wil je oefenen, dan is er online een berg geospatiale data beschikbaar om mee aan de slag te gaan. Een voorbeeld is de weerwebsite die ik hierboven heb genoemd.

Voor bronnen, kijk in dit boek: “Python Geospatial Development Essentials” van Karim Bahgat (2015) voor diepgaandere begeleiding.

Als je contact met me wilt opnemen, kun je me een e-mail sturen op dqvu.ubc@gmail.com of verbinden via LinkedIn. Veel leerplezier en blijf veilig!

Wil je meer leren over geospatiale data in Python, volg dan DataCamp's cursus Visualizing Geospatial Data in Python.

Hieronder staat de data die in deze tutorial is gebruikt:

• export.csv
• florence.csv
• gz_2010_us_040_00_5m
• hurricane_data