 # Spatial Data With Python - Part 1

Last updated: February 24th, 2020
In [ ]:
import rasterio
import matplotlib.pyplot as plt
from rasterio.plot import show


## Raster Files¶

We'll start exploring Raster files. The one that we're using for our example comes from the NASA's Visible Earth project (download link). It's already downloaded, so there's nothing you need to worry about 😁.

The rasterio library is used to read and process Raster files, and it should be already installed.

In [ ]:
dataset = rasterio.open("land_shallow_topo_2048.tif")


We can check the most important properties:

In [ ]:
print("count: {}".format(dataset.count)) # Raster bands
print("height: {}".format(dataset.height)) # Rows
print("width: {}".format(dataset.width)) # Columns


And let's explore one band as an example:

In [ ]:
band_1 = dataset.read(1) # Band 1

print(type(band_1))
print(band_1.shape)
print(band_1)


We can show the entire Raster file with the help of matplotlib:

In [ ]:
fig, ax = plt.subplots(figsize=(15, 10));
show(dataset, ax=ax)


## Vectors¶

Vector files are already downloaded, and we'll use GeoPandas to import and process them:

In [ ]:
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

In [ ]:
# Read the shapefiles into GeoDataFrames
world_gdf = gpd.read_file("TM_WORLD_BORDERS-0.3/TM_WORLD_BORDERS-0.3.shp")
volcanoes_gdf = gpd.read_file("harvard-glb-volc-shapefile/GLB_VOLC.shp")

In [ ]:
print("world_gdf shape {}".format(world_gdf.shape))
print("volcanoes_gdf shape {}".format(volcanoes_gdf.shape))


This is what each one of the files contain.

World data:

In [ ]:
world_gdf.head()

In [ ]:
world_gdf.shape


We see that the world map contains a variety of polygons, multipolygons and others. When plotted, we can see the whole world represented in this shape file:

In [ ]:
fig, ax = plt.subplots(figsize=(15, 10));
world_gdf.plot(ax=ax)


Volcanoes:

In [ ]:
volcanoes_gdf.head()

In [ ]:
volcanoes_gdf.shape


Volcanoes will be just a list of points

In [ ]:
fig, ax = plt.subplots(figsize=(15, 10));
volcanoes_gdf.plot(color="red", ax=ax)


### Putting everything together¶

We can combine our world data with the information of individual volcanoes to see everything together:

In [ ]:
fig, ax = plt.subplots(figsize=(15, 10));
world_gdf.plot(color="white", edgecolor="black", ax=ax)
volcanoes_gdf.plot(color="red", ax=ax)


Impressive, right! 🎉