Post #10: Advanced Geospatial Analysis with Remote Sensing Data in Python

 

Title: Unlocking Insights with Remote Sensing Data for Geospatial Analysis in Python

---

📍Introduction

Remote sensing involves capturing data from a distance, often using satellite imagery or aerial sensors. This data is incredibly valuable for analyzing land cover, vegetation health, urban development, and many other aspects of the Earth’s surface.

In this post, we’ll:

• Learn how to process and analyze satellite imagery using Python.

• Apply advanced geospatial techniques to analyze land cover, vegetation, and urban expansion using Remote Sensing Data.

• Explore tools like Sentinel-2 imagery with Rasterio and Geospatial Libraries.

---

🧰 Step 1: Install Required Libraries

To get started with remote sensing data analysis, you’ll need the following libraries:

bash
pip install rasterio geopandas matplotlib numpy scikit-image

• Rasterio: For reading and writing geospatial raster data (e.g., satellite imagery).

• GeoPandas: For handling vector data.

• Matplotlib: For plotting the data.

• Scikit-image: For image processing.

---

🗺️ Step 2: Load and Inspect Satellite Imagery

We’ll start by loading satellite imagery data. For example, we might use Sentinel-2 imagery, which provides multi-spectral data, including visible, infrared, and other bands.

python
import rasterio
import matplotlib.pyplot as plt

# Open the Sentinel-2 image (adjust path to your file)
img_path = "data/sentinel_image.tif"
src = rasterio.open(img_path)

# Inspect the image's metadata (e.g., number of bands, CRS)
print(src.meta)

# Read the bands (Sentinel-2 usually has 13 bands, but we'll use just a few for visualization)
red_band = src.read(4)  # Red band
green_band = src.read(3)  # Green band
blue_band = src.read(2)  # Blue band

# Visualize the RGB image
rgb = np.dstack((red_band, green_band, blue_band))

plt.figure(figsize=(10, 10))
plt.imshow(rgb)
plt.title("RGB Composition of Sentinel-2 Image")
plt.axis("off")
plt.show()

---

🧠 What Just Happened?

• We opened a GeoTIFF file (a popular format for satellite imagery).

• We accessed the red, green, and blue bands and combined them to create an RGB image, which is a common way of visualizing satellite imagery.

• Metadata inspection allows you to understand the structure of the image, including the number of bands, spatial resolution, and coordinate reference system (CRS).

---

📍 Step 3: Preprocess Satellite Data (Normalization and Clipping)

Satellite images often require preprocessing, like normalization or clipping, to enhance analysis.

Normalization: Adjust the pixel values (typically between 0 and 1).

python
import numpy as np

# Normalize the bands (min-max normalization)
red_band_norm = (red_band - red_band.min()) / (red_band.max() - red_band.min())
green_band_norm = (green_band - green_band.min()) / (green_band.max() - green_band.min())
blue_band_norm = (blue_band - blue_band.min()) / (blue_band.max() - blue_band.min())

# Create a normalized RGB image
rgb_norm = np.dstack((red_band_norm, green_band_norm, blue_band_norm))

# Plot the normalized image
plt.figure(figsize=(10, 10))
plt.imshow(rgb_norm)
plt.title("Normalized RGB Composition of Sentinel-2 Image")
plt.axis("off")
plt.show()

Clipping: Focus on a region of interest (ROI) by clipping the image.

python
from rasterio.mask import mask

# Define a bounding box or polygon for your region of interest (ROI)
# (this can come from your GIS dataset or a manually defined shape)
polygon = [
    {
        "type": "Polygon",
        "coordinates": [[
            [-5.0, 42.0],
            [-4.5, 42.0],
            [-4.5, 42.5],
            [-5.0, 42.5],
            [-5.0, 42.0]
        ]]
    }
]

# Clip the image with the ROI
out_image, out_transform = mask(src, polygon, crop=True)

# Display the clipped image (just the red band for simplicity)
plt.imshow(out_image[0], cmap='Reds')
plt.title("Clipped Red Band")
plt.axis("off")
plt.show()

---

🧠 What Just Happened?

• Normalization: We normalized the satellite bands to scale the pixel values between 0 and 1, which is useful for many analyses and visualizations.

• Clipping: We focused on a specific region of interest by applying a polygonal mask to the satellite image, limiting the analysis to just the relevant area.

---

📍 Step 4: Perform Land Cover Classification with NDVI

A common analysis in remote sensing is land cover classification. For example, we can use the Normalized Difference Vegetation Index (NDVI) to classify areas based on vegetation.

The NDVI is calculated using the red and near-infrared (NIR) bands, and its formula is:

$$\text{NDVI} = \frac{\text{NIR} - \text{Red}}{\text{NIR} + \text{Red}}$$

For this example, we’ll use the NIR band (Band 8) from Sentinel-2.

python
# Read the NIR band (usually Band 8 for Sentinel-2)
nir_band = src.read(8)

# Calculate NDVI
ndvi = (nir_band - red_band) / (nir_band + red_band)

# Plot NDVI
plt.figure(figsize=(10, 10))
plt.imshow(ndvi, cmap='RdYlGn')
plt.colorbar(label="NDVI")
plt.title("NDVI - Vegetation Index")
plt.axis("off")
plt.show()

---

🧠 What Just Happened?

• We calculated the NDVI to classify vegetation, which is a widely used index in remote sensing.

• Areas with high NDVI values represent vegetation (green), while low values represent non-vegetated areas.

---

📍 Step 5: Change Detection with Remote Sensing

Another common technique in remote sensing is change detection, where we compare two images taken at different times to identify changes in land cover, vegetation, or urbanization.

To perform change detection, we can calculate the difference in NDVI between two images taken at different times.

python
# Load the second image (e.g., a more recent Sentinel-2 image)
img_path_2 = "data/sentinel_image_2.tif"
src_2 = rasterio.open(img_path_2)

# Read the red and NIR bands from the second image
red_band_2 = src_2.read(4)
nir_band_2 = src_2.read(8)

# Calculate NDVI for the second image
ndvi_2 = (nir_band_2 - red_band_2) / (nir_band_2 + red_band_2)

# Calculate the change in NDVI between the two images
ndvi_change = ndvi_2 - ndvi

# Plot the NDVI change
plt.figure(figsize=(10, 10))
plt.imshow(ndvi_change, cmap='coolwarm')
plt.colorbar(label="NDVI Change")
plt.title("NDVI Change Detection")
plt.axis("off")
plt.show()

---

🧠 What Just Happened?

• Change Detection: We compared two images from different times and calculated the change in NDVI values to detect vegetation loss, urban expansion, or other changes.

---

📍 Step 6: Advanced Analysis – Object-Based Image Analysis (OBIA)

For even more advanced analysis, we can use Object-Based Image Analysis (OBIA), which segments the image into objects (e.g., land parcels, water bodies) and then classifies these objects based on various characteristics like shape, texture, and spectral properties.

python
from skimage.measure import label, regionprops

# Convert the NDVI values to binary for a basic example (vegetated or not)
binary_ndvi = ndvi > 0.3  # Vegetated areas

# Label the connected components (objects)
labeled = label(binary_ndvi)

# Get properties of each object (region)
regions = regionprops(labeled)

# Visualize the labeled regions
plt.imshow(labeled, cmap='tab20')
plt.title("Object-Based Segmentation (Labeled Regions)")
plt.axis("off")
plt.show()

---

🧠 What Just Happened?

• Object-Based Segmentation: We segmented the image into distinct objects (e.g., groups of pixels that form land features) and classified them based on NDVI.

• This can be used for more detailed analysis, like identifying urban areas, water bodies, or agricultural land.

---

🎯 Conclusion

Remote sensing opens up a wide range of possibilities for geospatial analysis, allowing you to:

• Analyze land cover and vegetation health using indices like NDVI.

• Detect spatial changes over time through change detection.

• Perform advanced image analysis using OBIA to identify distinct land features.

These techniques are invaluable for applications such as urban planning, environmental monitoring, and agriculture.

---

📌 Next Up:

➡️ Post 11: Working with LIDAR Data for 3D Geospatial Analysis