Rasterio - inspection, plotting raster data¶

The open geospatial Python ecosystem has powerful libraries for processing vector data. We have seen the likes of geopandas, pyshp, shapely, folium, fiona, pysal etc. How about for raster data processing? This notebook explores a few of them, starting with rasterio. My understanding so far is, in general digital image and signal processing via Python is a well established field. Thus if a library allows you to read geospatial image formats into numpy arrays and write the results back, you can easily accomplish your processing using standard Pythonic data analysis process of reading into numpy arrays, doing filtering, enhancements, analysis and presisting into the native file format of choice.

Let us see how far this theory holds up.

In [1]:

45+677

45+677

Out[1]:

722

Accessing raster files¶

In [1]:

import rasterio
import matplotlib.pyplot as plt
%matplotlib inline

import rasterio import matplotlib.pyplot as plt %matplotlib inline

In [2]:

raster = rasterio.open(path='/Users/atma6951/Documents/GIS_data/Imagery/open-geo/Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif')
type(raster)

raster = rasterio.open(path='/Users/atma6951/Documents/GIS_data/Imagery/open-geo/Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif') type(raster)

Out[2]:

rasterio._io.RasterReader

Let us query a few important properties to understand the image better

In [3]:

raster.crs

raster.crs

Out[3]:

CRS({'proj': 'tmerc', 'lat_0': 0, 'lon_0': -183, 'k': 0.9996, 'x_0': 500000, 'y_0': 0, 'datum': 'WGS84', 'units': 'm', 'no_defs': True})

In [4]:

raster.shape

raster.shape

Out[4]:

(1439, 1288)

In [5]:

raster.bounds

raster.bounds

Out[5]:

BoundingBox(left=698592.0, bottom=6656859.0, right=735300.0, top=6697870.5)

In [6]:

raster.count

raster.count

Out[6]:

7

In [8]:

raster.profile

raster.profile

Out[8]:

{'driver': 'GTiff',
 'dtype': 'uint8',
 'nodata': None,
 'width': 1288,
 'height': 1439,
 'count': 7,
 'crs': CRS({'proj': 'tmerc', 'lat_0': 0, 'lon_0': -183, 'k': 0.9996, 'x_0': 500000, 'y_0': 0, 'datum': 'WGS84', 'units': 'm', 'no_defs': True}),
 'transform': (698592.0, 28.5, 0.0, 6697870.5, 0.0, -28.5),
 'affine': Affine(28.5, 0.0, 698592.0,
        0.0, -28.5, 6697870.5),
 'tiled': False,
 'interleave': 'pixel'}

Read individual bands¶

In [10]:

band1 = raster.read([1])
type(band1)

band1 = raster.read([1]) type(band1)

Out[10]:

numpy.ndarray

In [11]:

band1.shape, band1.ndim

band1.shape, band1.ndim

Out[11]:

((1, 1439, 1288), 3)

In [19]:

band1.min(), band1.max(), band1.mean()

band1.min(), band1.max(), band1.mean()

Out[19]:

(0, 255, 59.63132232528628)

Read into a specified shape¶

I am unable to get this into a 2D array no matter what. Could be a bug

In [14]:

band1a = raster.read([1], out_shape=(raster.profile['width'], raster.profile['height']))
band1a.reshape(1288, 1439)
band1a.shape

band1a = raster.read([1], out_shape=(raster.profile['width'], raster.profile['height'])) band1a.reshape(1288, 1439) band1a.shape

Out[14]:

(1, 1288, 1439)

Preview the band¶

In [15]:

band1a

band1a

Out[15]:

array([[[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ...,
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]]], dtype=uint8)

In [16]:

from rasterio.plot import show
show(band1a)

from rasterio.plot import show show(band1a)

Out[16]:

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

Making a FCC or TCC

In [18]:

fig, ax = plt.subplots(figsize=(8,8))
show(band1a, cmap='Blues', ax=ax, alpha=0.33)
show(raster.read(3), cmap='Greens', ax=ax, alpha=0.33)
show(raster.read(4), cmap="Reds", ax=ax, alpha=0.33)

fig, ax = plt.subplots(figsize=(8,8)) show(band1a, cmap='Blues', ax=ax, alpha=0.33) show(raster.read(3), cmap='Greens', ax=ax, alpha=0.33) show(raster.read(4), cmap="Reds", ax=ax, alpha=0.33)

Out[18]:

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