🔄 Batch Append Feature Classes from Multiple GDBs into a Final Data Model (ArcPy)

 

🔄 Batch Append Feature Classes from Multiple GDBs into a Final Data Model (ArcPy)

In GIS data integration workflows, especially when working with multiple sources contributing to the same data model, it's common to append data from various geodatabases into a final unified schema.

This script automates that process using ArcPy, respecting both:

• Dataset-based feature classes, and

• Root-level feature classes.

---

📁 Folder Setup

• final_gdb: The destination GDB that contains your field-validated feature classes (the data model).

• folder_path: A folder with multiple source .gdb files. These should match the structure of the final GDB (names must match for auto-append to work).

---

🐍 ArcPy Script

python
import arcpy
import os

# Path to the geodatabase containing the final feature class models (with correct fields)
final_gdb = r"C:\Work\Projects\Zatca\UG_Data\5Th\Merge\Batha_Total.gdb"

# Folder containing the GDB files to append from
folder_path = r"C:\Work\Projects\Zatca\UG_Data\5Th\DataBaseBathaNew"

# Get a list of all datasets and feature classes in the final GDB (these are the models)
arcpy.env.workspace = final_gdb
final_datasets = arcpy.ListDatasets() or []
final_feature_classes = arcpy.ListFeatureClasses()

# Iterate through each GDB in the folder containing the source feature classes
for gdb_name in os.listdir(folder_path):
    gdb_path = os.path.join(folder_path, gdb_name)

    # Check if the item is a valid File Geodatabase
    if os.path.isdir(gdb_path) and gdb_name.endswith('.gdb'):
        arcpy.env.workspace = gdb_path

        # List all datasets and feature classes in the current source GDB
        source_datasets = arcpy.ListDatasets() or []
        source_feature_classes = arcpy.ListFeatureClasses()

        # Process datasets
        for final_dataset in final_datasets:
            if final_dataset in source_datasets:
                final_dataset_path = os.path.join(final_gdb, final_dataset)
                source_dataset_path = os.path.join(gdb_path, final_dataset)

                # Set workspace to dataset for feature class listing
                arcpy.env.workspace = source_dataset_path
                source_fc_list = arcpy.ListFeatureClasses()

                arcpy.env.workspace = final_dataset_path
                final_fc_list = arcpy.ListFeatureClasses()

                for final_fc in final_fc_list:
                    if final_fc in source_fc_list:
                        source_fc_path = os.path.join(source_dataset_path, final_fc)
                        final_fc_path = os.path.join(final_dataset_path, final_fc)

                        arcpy.Append_management([source_fc_path], final_fc_path, "NO_TEST")
                        print(f"✅ Appended data from '{source_fc_path}' to '{final_fc_path}'.")

        # Process individual feature classes not in datasets
        arcpy.env.workspace = gdb_path
        for final_fc in final_feature_classes:
            if final_fc in source_feature_classes:
                source_fc_path = os.path.join(gdb_path, final_fc)
                final_fc_path = os.path.join(final_gdb, final_fc)

                arcpy.Append_management([source_fc_path], final_fc_path, "NO_TEST")
                print(f"✅ Appended data from '{source_fc_path}' to '{final_fc_path}'.")

---

✅ Good To Know

• This script assumes schema consistency between source and target FCs.

• "NO_TEST" option tells ArcPy not to check schema compatibility (use only if you've pre-validated).

• You can switch to "TEST" if you want to enforce field-by-field schema checks.

📝 Bulk Update Field Aliases in All Feature Classes (ArcPy + CSV)

 

📝 Bulk Update Field Aliases in All Feature Classes (ArcPy + CSV)

Keeping your GIS field names meaningful is important—but sometimes they can be cryptic due to system constraints (e.g., LU_CODE, ADM_NAME). That’s where aliases come in: they provide readable labels without changing the field's actual name.

Here’s a simple ArcPy script that updates field aliases across all feature classes—including those inside datasets—using a CSV file.

---

🧾 Input: CSV Format

Make sure your CSV looks like this:

csv
FieldName,AliasName
LU_CODE,Land Use Code
ADM_NAME,Administrative Name
CITY_ID,City Identifier

This approach is great for global alias updates when you know the field names but don’t want to repeat entries for every feature class.

---

💻 Python Script (ArcPy)

python
import arcpy
import csv
import os

# === INPUT PARAMETERS ===
gdb_path = r"C:\Users\Testuser\OneDrive - GPC Global Information Solutions LLC\Projects\MPDA\Data\MPDA_MigratedData_25_03\MPDA_MigratedData_1.gdb"
csv_path = r"C:\Users\Testuser\Downloads\Updated_Alias_Names.csv"

# === LOAD CSV TO DICTIONARY ===
field_alias_dict = {}
with open(csv_path, newline='', encoding='utf-8') as csvfile:
    reader = csv.DictReader(csvfile)
    for row in reader:
        field_name = row['FieldName'].strip()
        alias_name = row['AliasName'].strip()
        field_alias_dict[field_name] = alias_name

# === PROCESS FEATURE CLASSES IN DATASETS ===
arcpy.env.workspace = gdb_path
datasets = arcpy.ListDatasets(feature_type='feature') or []

# Include root-level feature classes
datasets.append("")

for ds in datasets:
    feature_classes = arcpy.ListFeatureClasses(feature_dataset=ds)
    for fc in feature_classes:
        fc_path = os.path.join(gdb_path, ds, fc) if ds else os.path.join(gdb_path, fc)
        print(f"🔄 Processing: {fc_path}")

        fields = arcpy.ListFields(fc_path)
        for field in fields:
            if field.name in field_alias_dict:
                new_alias = field_alias_dict[field.name]
                try:
                    print(f"  ➤ Updating alias for '{field.name}' to '{new_alias}'")
                    arcpy.AlterField_management(
                        in_table=fc_path,
                        field=field.name,
                        new_field_alias=new_alias
                    )
                except Exception as e:
                    print(f"  ⚠️ Failed to update field '{field.name}': {e}")

print("\n✅ Alias name update completed.")

---

🔧 Use Cases

• Renaming confusing or auto-generated field names in bulk

• Preparing datasets for external sharing (with readable labels)

• Standardizing field aliases across hundreds of layers

✏️ Bulk Update Field Aliases in a Geodatabase Using ArcPy and CSV

 

✏️ Bulk Update Field Aliases in a Geodatabase Using ArcPy and CSV

Maintaining clear and readable aliases for fields in your GIS layers is vital for both data clarity and user experience. Whether you’re preparing datasets for a client, public release, or internal documentation, consistent aliases make a big difference.

Here’s a handy Python script using ArcPy that lets you batch-update field aliases using a CSV file—no more clicking through each field in ArcGIS Pro!

---

🧾 What the Script Does

• Reads a CSV file containing:

  • Dataset name (optional)

  • Feature class name

  • Field name

  • New alias

• Locates the feature class in the specified .gdb

• Applies the new alias using AlterField_management

---

📦 Folder Setup

Your project folder might look like this:

objectivec
MPDA_MigratedData_25_03/
│
├── MPDA_MigratedData.gdb/
│     └── (your feature classes & datasets)
│
└── BackUP/
      └── AliasNamesToUpdate.csv

---

📄 Sample CSV Structure

csv
DatasetName,FeatureClassName,FieldName,AliasName
,CityBoundaries,City_ID,City Identifier
LandUse,Parcels,LU_CODE,Land Use Code
Admin,Regions,ADM_NAME,Administrative Name

If the feature class is at the root level of the GDB (not inside a dataset), leave DatasetName blank.

---

💻 Python Script (ArcPy)

python
import arcpy
import csv
import os
import codecs

# === INPUTS ===
gdb_path = r"C:\Path\To\Your.gdb"
csv_path = r"C:\Path\To\AliasNamesToUpdate.csv"

# Read CSV with UTF-8-SIG (handles BOM from Excel)
with codecs.open(csv_path, 'r', encoding='utf-8-sig') as csvfile:
    reader = csv.DictReader(csvfile)
    reader.fieldnames = [col.strip() for col in reader.fieldnames]
    print("🟢 Headers:", reader.fieldnames)

    for row in reader:
        dataset_name = row['DatasetName'].strip()
        feature_class = row['FeatureClassName'].strip()
        field_name = row['FieldName'].strip()
        alias_name = row['AliasName'].strip()

        # Get full path to feature class
        full_fc_path = os.path.join(gdb_path, dataset_name, feature_class) if dataset_name else os.path.join(gdb_path, feature_class)

        if arcpy.Exists(full_fc_path):
            print(f"🔄 Processing: {full_fc_path} ➤ Field: {field_name}")
            try:
                arcpy.AlterField_management(
                    in_table=full_fc_path,
                    field=field_name,
                    new_field_alias=alias_name
                )
                print(f"✅ Updated alias to: {alias_name}")
            except Exception as e:
                print(f"⚠️ Failed to update alias for {field_name}: {e}")
        else:
            print(f"❌ Feature class not found: {full_fc_path}")

---

🔄 Why Use This?

• Speeds up metadata cleanup

• Enables non-GIS staff to manage aliases via Excel

• Reduces manual errors

• Keeps alias naming consistent across large datasets

---

🧠 Pro Tips

• Test on a copy of your geodatabase first.

• To also rename fields, use new_field_name= in the same function.

• Works in ArcGIS Pro (Python 3) or ArcMap (Python 2)—but always use the version that matches your .gdb.

🧩 Add Standard Metadata Fields to All Feature Classes in Multiple Geodatabases Using ArcPy

 

🧩 Add Standard Metadata Fields to All Feature Classes in Multiple Geodatabases Using ArcPy

When working on large GIS projects—especially in urban planning, infrastructure, or government contexts—it's essential to maintain consistent metadata across all your datasets.

This script automates the process of adding predefined fields to every feature class in multiple file geodatabases (.gdb) inside a folder. It's perfect for enforcing a standard schema and ensuring that downstream users and systems have the metadata they need.

---

🎯 What the Script Does

• Scans a folder for all .gdb files (file geodatabases).

• For each .gdb:

  • Adds a set of user-defined fields to all standalone feature classes.

  • Adds the same fields to all feature classes inside datasets.

• Skips adding a field if it already exists.

---

📜 The Python Script (ArcPy)

python
import arcpy
import os

# 📁 Folder containing your .gdb files
gdb_folder = r'C:\Path\To\Your\GeodatabaseFolder'

# 🧾 Fields to be added to each feature class
fields_to_add = [
    ("OriginalName", "TEXT", 255),
    ("OriginalPath", "TEXT", 1000),
    ("DIA_ID", "TEXT", 255),
    ("DIA_Date", "Date", 255)
]

# 🧰 Function to add fields to a single feature class
def add_fields_to_fc(feature_class):
    for field_name, field_type, field_length in fields_to_add:
        fields = arcpy.ListFields(feature_class, field_name)
        if not fields:
            if field_length and field_type.upper() == "TEXT":
                arcpy.AddField_management(feature_class, field_name, field_type, field_length=field_length)
            else:
                arcpy.AddField_management(feature_class, field_name, field_type)
            print(f"✅ Added field '{field_name}' of type {field_type} to {feature_class}")
        else:
            print(f"ℹ️ Field '{field_name}' already exists in {feature_class}")

# 🔁 Loop through all geodatabases in the folder
for folder in os.listdir(gdb_folder):
    if folder.endswith(".gdb"):
        gdb_path = os.path.join(gdb_folder, folder)
        arcpy.env.workspace = gdb_path
        print(f"\n📦 Processing Geodatabase: {gdb_path}")

        # ➕ Standalone Feature Classes
        feature_classes = arcpy.ListFeatureClasses()
        for fc in feature_classes:
            print(f"🔹 Standalone FC: {fc}")
            add_fields_to_fc(fc)

        # 📁 Feature Datasets
        datasets = arcpy.ListDatasets(feature_type='Feature')
        for dataset in datasets:
            dataset_fcs = arcpy.ListFeatureClasses(feature_dataset=dataset)
            for fc in dataset_fcs:
                print(f"🔸 FC in Dataset '{dataset}': {fc}")
                add_fields_to_fc(fc)

print("\n🎉 Metadata field addition complete.")

---

🧠 Use Cases

• Standardizing data before submission to a national GIS database.

• Adding internal tracking fields to all project layers.

• Preparing spatial data for external collaboration or audit.

---

🔧 Customizing the Script

You can change or expand the fields_to_add list as needed:

python
fields_to_add = [
    ("NewField", "TEXT", 50),
    ("ImportDate", "DATE", None),
    ("Status", "TEXT", 25)
]

---

💡 Tips

• Want to calculate values too? Use arcpy.CalculateField_management() after adding the field.

• Running in ArcGIS Pro? Make sure your Python environment is set to arcgispro-py3.

• For batch runs, wrap this into a .bat file and run on schedule using Task Scheduler.

📏 Auto-Add Area and Length Fields in a Geodatabase Using ArcPy

 

📏 Auto-Add Area and Length Fields in a Geodatabase Using ArcPy

When managing spatial data in ArcGIS, calculating geometric properties like area or length is a common requirement—especially for analytics, reporting, or quality control.

In this post, you'll learn how to write a Python script using ArcPy that automatically adds:

• a CALC_AREA field to all Polygon feature classes

• a CALC_LENGTH field to all Polyline feature classes

All within a given File Geodatabase (.gdb)—whether the feature classes are inside datasets or at the root level.

---

🧾 The Script

python
import arcpy
import os

# 📍 Set your input File Geodatabase path
input_gdb = r"C:\Path\To\Your\DataModel.gdb"

# Set the workspace
arcpy.env.workspace = input_gdb

# 🔍 List all feature datasets (or include root-level FCS by using [''] if none)
datasets = arcpy.ListDatasets(feature_type='feature') or ['']

# Loop through all datasets and standalone feature classes
for dataset in datasets:
    dataset_path = os.path.join(input_gdb, dataset) if dataset else input_gdb
    feature_classes = arcpy.ListFeatureClasses(feature_dataset=dataset)

    for fc in feature_classes:
        fc_path = os.path.join(dataset_path, fc) if dataset else os.path.join(input_gdb, fc)

        # 🧠 Get geometry type (Polygon or Polyline)
        desc = arcpy.Describe(fc_path)
        geom_type = desc.shapeType.upper()

        # 🏗️ Define field to add based on geometry type
        if geom_type == "POLYGON":
            field_name, field_type = "CALC_AREA", "DOUBLE"
        elif geom_type == "POLYLINE":
            field_name, field_type = "CALC_LENGTH", "DOUBLE"
        else:
            continue  # Skip Points or unsupported types

        # ✅ Check if field exists before adding
        existing_fields = [f.name for f in arcpy.ListFields(fc_path)]
        if field_name not in existing_fields:
            arcpy.AddField_management(fc_path, field_name, field_type)
            print(f"✅ Added field '{field_name}' to '{fc_path}'")
        else:
            print(f"ℹ️ Field '{field_name}' already exists in '{fc_path}', skipping.")

print("\n🎉 Geometry field addition complete.")

---

🛠️ How It Works

• ListDatasets() gets all feature datasets in your geodatabase.

• ListFeatureClasses() gets each feature class inside them (or in the root).

• Describe() helps us detect if a feature class is a Polygon or Polyline.

• If the appropriate field (CALC_AREA or CALC_LENGTH) is missing, it gets added as a DOUBLE.

---

📦 Why Use This?

• You're preparing data for area/length calculations.

• You want to standardize schema across multiple datasets.

• You're doing bulk data migration or preparing for a data audit.

---

✍️ Bonus Tip: Calculate Values

If you'd like to calculate values right after adding the field, insert this below the AddField_management line:

python
if field_name == "CALC_AREA":
    arcpy.CalculateField_management(fc_path, field_name, "!shape.area!", "PYTHON3")
elif field_name == "CALC_LENGTH":
    arcpy.CalculateField_management(fc_path, field_name, "!shape.length!", "PYTHON3")

🔍 Automating Field Addition in File Geodatabases using ArcPy

 

🔍 Automating Field Addition in File Geodatabases using ArcPy

Managing large GIS data repositories often requires performing the same operation across multiple datasets and feature classes. Whether you're maintaining metadata, tracking source information, or preparing datasets for analysis, automating repetitive tasks can save a lot of time.

In this post, I’ll show you a Python script using ArcPy that automatically adds a predefined set of fields to all feature classes in all File Geodatabases (.gdb) within a given folder. This is especially useful when you're dealing with standardized data models or working on large-scale projects with multiple data sources.

---

✅ What This Script Does

• Searches for all .gdb (File Geodatabases) in a specified folder.

• For each geodatabase:

  • Adds fields to standalone feature classes.

  • Adds fields to feature classes inside feature datasets.

• Checks if each field already exists before trying to add it.

• Handles both TEXT and DATE field types.

---

🧠 How to Use It

You’ll need:

• ArcGIS Pro or ArcMap with ArcPy available.

• A folder path where your .gdb files are stored.

Update the gdb_folder and fields_to_add as needed.

---

🧾 The Script

python
import arcpy
import os

# 📁 Set the directory where your .gdb files are stored
gdb_folder = r'C:\Path\To\Your\Geodatabases'

# 📝 Define the fields to be added: (field name, type, optional length)
fields_to_add = [
    ("OriginalName", "TEXT", 255),
    ("OriginalPath", "TEXT", 1000),
    ("DIA_ID", "TEXT", 255),
    ("DIA_Date", "DATE", None)
]

# 🚀 Function to add fields to a given feature class
def add_fields_to_feature_class(feature_class):
    for field_name, field_type, field_length in fields_to_add:
        try:
            existing_fields = arcpy.ListFields(feature_class, field_name)
            if not existing_fields:
                if field_length and field_type.upper() == "TEXT":
                    arcpy.AddField_management(feature_class, field_name, field_type, field_length=field_length)
                else:
                    arcpy.AddField_management(feature_class, field_name, field_type)
                print(f"✅ Added field '{field_name}' to {feature_class}")
            else:
                print(f"ℹ️ Field '{field_name}' already exists in {feature_class}")
        except Exception as e:
            print(f"❌ Error adding field '{field_name}' to {feature_class}: {e}")

# 🔄 Loop through each .gdb in the folder
for item in os.listdir(gdb_folder):
    if item.endswith(".gdb"):
        gdb_path = os.path.join(gdb_folder, item)
        arcpy.env.workspace = gdb_path
        print(f"\n📂 Processing Geodatabase: {gdb_path}")

        # Process standalone feature classes
        fcs = arcpy.ListFeatureClasses()
        if fcs:
            for fc in fcs:
                print(f"🔹 Feature Class: {fc}")
                add_fields_to_feature_class(fc)
        else:
            print("⚠️ No standalone feature classes found.")

        # Process feature datasets
        datasets = arcpy.ListDatasets(feature_type='Feature')
        if datasets:
            for ds in datasets:
                print(f"\n📁 Dataset: {ds}")
                ds_fcs = arcpy.ListFeatureClasses(feature_dataset=ds)
                if ds_fcs:
                    for fc in ds_fcs:
                        print(f"🔸 Feature Class in Dataset: {fc}")
                        add_fields_to_feature_class(fc)
                else:
                    print(f"⚠️ No feature classes in dataset '{ds}'")
        else:
            print("⚠️ No feature datasets found.")

print("\n🎉 Field addition completed for all geodatabases.")

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