RASTERIO¶

By: Danraj Lamichhane¶

1.Introduction to Rasterio¶

Before Rasterio there was one Python option for accessing the many different kind of raster data files used in the GIS field: the Python bindings distributed with the Geospatial Data Abstraction Library, GDAL. These bindings extend Python, but provide little abstraction for GDAL’s C API. This means that Python programs using them tend to read and run like C programs.

For more info: https://rasterio.readthedocs.io/en/latest/intro.html

2. Installation and Setup¶

Rasterio works with Python 3.8+, Numpy 1.18+, and GDAL 3.1+.

• pip install rasterio

• conda install -c conda-forge rasterio

3. Using Rasterio¶

In [1]:

# Import Rasterio
import rasterio
from rasterio.plot import show

# Import Rasterio import rasterio from rasterio.plot import show

In [2]:

# Read a DEM file
src=rasterio.open('DEM.tif')

# Read a DEM file src=rasterio.open('DEM.tif')

In [3]:

show(src)

show(src)

Out[3]:

<Axes: >

In [4]:

src.meta #Reads the metadata

src.meta #Reads the metadata

Out[4]:

{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -9999.0,
 'width': 10814,
 'height': 21613,
 'count': 1,
 'crs': CRS.from_epsg(4326),
 'transform': Affine(0.0002777777780359812, 0.0, -76.00194669160615,
        0.0, -0.00027777777803598107, 43.00195317827816)}

In [5]:

src.bounds # Boounding Box

src.bounds # Boounding Box

Out[5]:

BoundingBox(left=-76.00194669160615, bottom=36.9983420615865, right=-72.99805779992505, top=43.00195317827816)

3.1 Cropping¶

In [6]:

from rasterio.windows import Window
#Cropping the DEM
window = Window(500,500 ,3000, 3000) # parameters are (col_off, row_off, width, height)
crp = src.read(window=window)
show(crp)

from rasterio.windows import Window #Cropping the DEM window = Window(500,500 ,3000, 3000) # parameters are (col_off, row_off, width, height) crp = src.read(window=window) show(crp)

Out[6]:

<Axes: >

You can crop a raster dataset using geographic coordinates in Rasterio, but it involves a few more steps compared to cropping with pixel coordinates. Since raster datasets are inherently grid-based and use row and column indices, you'll first need to convert the geographic coordinates to pixel coordinates.

3.2 Resampling¶

Resampling involves changing the spatial resolution of the raster data, either increasing (upsampling) or decreasing (downsampling) the number of pixels.

In [47]:

from rasterio.enums import Resampling
data = src.read(
     out_shape=(
        src.count,
        int(src.height * 0.01),  # decreasing the resolution
        int(src.width *0.01)
    ),
        resampling=Resampling.bilinear
)

    # Update the metadata to reflect the new shape
transform = src.transform * src.transform.scale(
    (src.width / data.shape[-2]),
    (src.height / data.shape[-1])
)

from rasterio.enums import Resampling data = src.read( out_shape=( src.count, int(src.height * 0.01), # decreasing the resolution int(src.width *0.01) ), resampling=Resampling.bilinear ) # Update the metadata to reflect the new shape transform = src.transform * src.transform.scale( (src.width / data.shape[-2]), (src.height / data.shape[-1]) )

In [48]:

show(data)

show(data)

Out[48]:

<Axes: >

In [50]:

resampled_height = data.shape[-1]  # Height of the resampled image
resampled_width = data.shape[-2]

resampled_height = data.shape[-1] # Height of the resampled image resampled_width = data.shape[-2]

In [51]:

resampled_height

resampled_height

Out[51]:

108

In [52]:

resampled_width

resampled_width

Out[52]:

216

3.3 Raster Calculation and Aanlysis¶

3.3.1 Slope Calculation:¶

Slope represents the rate of change in elevation at each pixel. It is a crucial factor in many geospatial analyses.

In [9]:

import numpy as np

import numpy as np

In [34]:

def calculate_slope(dem, pixel_size):
    """Calculate slope from a DEM."""
    dx, dy = np.gradient(dem, pixel_size[0], pixel_size[1])
    slope = np.arctan(np.sqrt(dx**2 + dy**2)) * (180 / np.pi)  # Convert to degrees
    return slope
affine=src.transform
dem=src.read(1)

def calculate_slope(dem, pixel_size): """Calculate slope from a DEM.""" dx, dy = np.gradient(dem, pixel_size[0], pixel_size[1]) slope = np.arctan(np.sqrt(dx**2 + dy**2)) * (180 / np.pi) # Convert to degrees return slope affine=src.transform dem=src.read(1)

In [35]:

slope = calculate_slope(dem, affine)

slope = calculate_slope(dem, affine)

/home/dl1197/miniconda3/envs/rcaes_env/lib/python3.9/site-packages/numpy/lib/function_base.py:1238: RuntimeWarning: divide by zero encountered in divide
  out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
/home/dl1197/miniconda3/envs/rcaes_env/lib/python3.9/site-packages/numpy/lib/function_base.py:1238: RuntimeWarning: invalid value encountered in divide
  out[tuple(slice1)] = (f[tuple(slice4)] - f[tuple(slice2)]) / (2. * ax_dx)
/home/dl1197/miniconda3/envs/rcaes_env/lib/python3.9/site-packages/numpy/lib/function_base.py:1259: RuntimeWarning: divide by zero encountered in divide
  out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
/home/dl1197/miniconda3/envs/rcaes_env/lib/python3.9/site-packages/numpy/lib/function_base.py:1259: RuntimeWarning: invalid value encountered in divide
  out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_0
/home/dl1197/miniconda3/envs/rcaes_env/lib/python3.9/site-packages/numpy/lib/function_base.py:1266: RuntimeWarning: divide by zero encountered in divide
  out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n
/home/dl1197/miniconda3/envs/rcaes_env/lib/python3.9/site-packages/numpy/lib/function_base.py:1266: RuntimeWarning: invalid value encountered in divide
  out[tuple(slice1)] = (f[tuple(slice2)] - f[tuple(slice3)]) / dx_n

In [36]:

print(slope) #Values are in degrees

print(slope) #Values are in degrees

[[nan nan nan ... nan nan nan]
 [90. 90. 90. ... nan nan nan]
 [90. 90. 90. ... nan nan nan]
 ...
 [90. 90. nan ... nan nan nan]
 [90. 90. nan ... nan nan nan]
 [90. 90. nan ... nan nan nan]]

3.3.2 Histogram¶

In [13]:

from rasterio.plot import show_hist

from rasterio.plot import show_hist

In [14]:

x=show_hist(crp, bins=10,stacked='False',edgecolor='black', alpha=0.7)

x=show_hist(crp, bins=10,stacked='False',edgecolor='black', alpha=0.7)

Other Functions¶

Cloud-Optimized GeoTIFF Support Read from Remote Sources: Directly read from and write to cloud-optimized GeoTIFFs, a standard format for cloud-based geospatial data storage.

If you are using gdal you can use most of the functions of gdal using rasterio.

Creating Rasters from arryas.

In [ ]: