Chapter 2 — Spatial Raster Data in Python cont.

Objectives are:

• Create new raster objects.
• Assign the correct projection or CRS.

A raster data model is an array of cells, or pixels, to represent real-world objects.

Raster datasets are commonly used for representing and managing imagery, surface temperatures, digital elevation models, and numerous other entities.

Here we create an example of two ndarray objects one X spans [-90°,90°] longitude, and Y covers [-90°,90°] latitude.

import numpy as np
x = np.linspace(-90, 90, 6)
y = np.linspace(90, -90, 6)
X, Y = np.meshgrid(x, y)
X

array([[-90., -54., -18.,  18.,  54.,  90.],
       [-90., -54., -18.,  18.,  54.,  90.],
       [-90., -54., -18.,  18.,  54.,  90.],
       [-90., -54., -18.,  18.,  54.,  90.],
       [-90., -54., -18.,  18.,  54.,  90.],
       [-90., -54., -18.,  18.,  54.,  90.]])

Let’s generate data representing temperature (Z )

import matplotlib.pyplot as plt

Z1 =  np.abs(((X - 10) ** 2 + (Y - 10) ** 2) / 1 ** 2)
Z2 =  np.abs(((X + 10) ** 2 + (Y + 10) ** 2) / 2.5 ** 2)
Z =  (Z1 - Z2)

plt.imshow(Z)
plt.title("Temperature")
plt.show()

Assigning spatial data to an array in python

How do we create a spatial dataset from this data?

We need three components:

1. An array of data and the xy coordinates.
2. A coordinate reference system.
3. A transform defining the coordinate of the upper left hand corner and the cell resolution.

Now, we can write out a working spatial raster dataset in python in a few lines of code. We just need to provide the information listed above in a format that rasterio understands.

from rasterio.transform import Affine
import rasterio as rio

res = (x[-1] - x[0]) / 240.0
transform = Affine.translation(x[0] - res / 2, y[0] - res / 2) * Affine.scale(res, res)

# open in 'write' mode, unpack profile info to dst
with rio.open(
   "../temp/new_raster.tif",
   "w",
   driver="GTiff",         # output file type
   height=Z.shape[0],      # shape of array
   width=Z.shape[1],
   count=1,                # number of bands
   dtype=Z.dtype,          # output datatype
   crs="+proj=latlong",    # CRS
   transform=transform,    # location and resolution of upper left cell
) as dst:
   # check for number of bands
   if dst.count == 1:
       # write single band
       dst.write(Z, 1)
   else:
       # write each band individually
       for band in range(len(Z)):
           # write data, band # (starting from 1)
           dst.write(Z[band], band + 1)

Note:

Raster data is often ‘multiband’ (e.g. red, green, blue), so I provided code that works for both multiband and single-band data (dst.count).

If you are storing multiband data the dimensions should be stored as (band, y, x ).

Open tif file and plot it:

from rasterio.plot import show
%matplotlib inline

# Open the file:
raster = rio.open("Data/new_raster.tif")

# Plot band 1
show((raster, 1))

Link GitHub for this chapter here.

References: pyGis, Geopandas, rasterio, and Wiki.