.ipynb .pdf

Contents

1. Extracting vegetation area

Written by Men Vuthy, 2021

---

Import modules

import os
import numpy as np
import pandas as pd

import rasterio
from rasterio.plot import show

import matplotlib.pyplot as plt
%matplotlib inline

# import data
raster = rasterio.open('image_indo.tiff')

# visualize red band
show(raster, 3, cmap = 'Reds')

<AxesSubplot:>

# Read the grid values into numpy arrays
nir = raster.read(4)
red = raster.read(3)
green = raster.read(2)
blue = raster.read(1)

Define normalize function to reduce scale to 0.0 - 1.0

# Function to normalize the grid values
def normalize(array):
    """Normalizes numpy arrays into scale 0.0 - 1.0"""
    array_min, array_max = np.nanmin(array), np.nanmax(array)
    return ((array - array_min)/(array_max - array_min))

Normalize all the bands

# Normalize the bands
nirn = normalize(nir)
redn = normalize(red)
greenn = normalize(green)
bluen = normalize(blue)

a. Calculate NDVI

Create a function to calculate NDVI

# Function to calculate ndvi
def get_ndvi(nir, red):
    # By default numpy will complain about dividing with zero values.
    # We need to change that behaviour because we have a lot of 0 values in our data.
    np.seterr(divide='ignore', invalid='ignore')
    # NDVI formula
    ndvi = (nir - red) / (nir + red)
    return ndvi

# Calculate ndvi
NDVI = get_ndvi(nirn, redn)

Visualize ndvi in a plot NDVI

# Plot the NDVI
plt.figure(figsize=(7,7))
plt.imshow(NDVI, cmap='RdYlGn')

# Add colorbar to show the index
plt.colorbar(fraction=0.03, pad=0.04)
plt.title('NDVI')

# plt.savefig('NDVI.png', dpi = 150)

plt.show();

Export ndvi to GeoTiff

# Data dir
data_dir = "output"

# Output raster
out_tif = os.path.join(data_dir, "NDVI.tif")

# Copy the metadata
out_meta = raster.meta.copy()
out_meta

# Update the metadata
out_meta.update({'driver': 'GTiff',
                 'dtype': 'float32',
                 'nodata': None,
                 'width': NDVI.shape[1],
                 'height': NDVI.shape[0],
                 'crs': raster.crs,
                 'count':1,
                 'transform': raster.transform
                })

with rasterio.open(out_tif, "w", **out_meta) as dest:
    dest.write(NDVI.astype(np.float32), indexes = 1)

b. Extract vegetation area

# Make a copy of above NDVI
import copy
vegt_area = copy.copy(NDVI)

# Set Threshold value (it's dependent on area. The value must be manually checked against NRG image in QGIS.)
vegt_area[NDVI<0.25] = np.nan

# Plot the NDVI
plt.figure(figsize=(7,7))
plt.imshow(NDVI, cmap='RdYlGn')

# Add colorbar to show the index
plt.colorbar(fraction=0.03, pad=0.04)
plt.title('NDVI')

# plt.savefig('NDVI.png', dpi = 150)

plt.show();

# Data dir
data_dir = "output"

# Output raster
out_tif = os.path.join(data_dir, "veg_area.tif")

# Copy the metadata
out_meta = raster.meta.copy()
out_meta

# Update the metadata
out_meta.update({'driver': 'GTiff',
                 'dtype': 'float32',
                 'nodata': None,
                 'width': vegt_area.shape[1],
                 'height': vegt_area.shape[0],
                 'crs': raster.crs,
                 'count':1,
                 'transform': raster.transform
                })

with rasterio.open(out_tif, "w", **out_meta) as dest:
    dest.write(vegt_area.astype(np.float32), indexes = 1)

c. Visualize result

import hvplot
import hvplot.xarray
import geoviews as gv
import xarray as xr

# Load data
image = xr.open_rasterio('image_indo.tiff')
ndvi = xr.open_rasterio('output/NDVI.tif')
veg = xr.open_rasterio('output/veg_area.tif')

# Convert datatype from Uint16 to Uint8
dnan = image.where(image>0)
norm = lambda arr : (arr-np.nanmin(arr)) * 255/(np.nanmax(arr)-np.nanmin(arr))
val = np.array([norm(dnan.values[0]), norm(dnan.values[1]), norm(dnan.values[2]), norm(dnan.values[3])])
val = np.where(np.isnan(val), 0, val)
val = val.astype(np.uint8)
dnan.values = val

# Add basemap
geomap = gv.tile_sources.EsriImagery

# Create composite images for natural and false color
natural = geomap *dnan.isel(band=[2,3,1]).hvplot.rgb(x='x',y='y',z='value', bands='band',geo=True, frame_height=400).options(title='natural color')
false = geomap *dnan.isel(band=[3,2,1]).hvplot.rgb(x='x',y='y',z='value', bands='band',geo=True, frame_height=400).options(title='false color')

natural + false

C:\Users\a9418\AppData\Roaming\Python\Python37\site-packages\pyproj\crs\crs.py:292: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  projstring = _prepare_from_string(projparams)
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
C:\Users\a9418\AppData\Roaming\Python\Python37\site-packages\pyproj\crs\crs.py:292: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  projstring = _prepare_from_string(projparams)
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace

Data type cannot be displayed:

# visualize ndvi
ndvi_plot = geomap * ndvi[0].hvplot.image(x='x',y='y',cmap='PRGn', geo=True, frame_height=400,clim=(-1,1)).options(title='NDVI', clipping_colors={'NaN':'transparent'})
ndvi_plot

C:\Users\a9418\AppData\Roaming\Python\Python37\site-packages\pyproj\crs\crs.py:292: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  projstring = _prepare_from_string(projparams)
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace

Data type cannot be displayed:

# visualize vegetation area
veg_plot = false * veg[0].hvplot.image(x='x',y='y',cmap='PRGn', geo=True, frame_height=400,clim=(-1,1)).options(title='Vegetation area', clipping_colors={'NaN':'transparent', '-3':'transparent'})
veg_plot

C:\Users\a9418\AppData\Roaming\Python\Python37\site-packages\pyproj\crs\crs.py:292: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
  projstring = _prepare_from_string(projparams)
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace
WARNING:param.ParameterizedMetaclass: Use method 'params' via param namespace

Data type cannot be displayed:

previous

Indonesia project: high-resolution satellite image analysis

next

2. Extracting water area