These code snippets illustrate how to convert coordinate data from plain text files into polygonal boundaries for administrative districts, then calculate the proportion of agricultural land cover within each polygon for two years. By leveraging Shapely and GeoPandas, it demonstrates straightforward ways to generate, analyze, and report vector and raster data metrics (e.g., number of agricultural pixels) in a cohesive workflow. 

import pandas as pd

import geopandas as gpd

from shapely.geometry import Polygon

import rasterio.features

import numpy as np

# Part I: Read District Files and Create Polygons

    # Read coordinates from text files and create polygons

    # Store district number, polygon geometry, and vertex count

# Read the district files, skipping the first row

df1 = pd.read_csv("district01.txt", delimiter="\t", names=["X", "Y"], skiprows=1)

df2 = pd.read_csv("district05.txt", delimiter="\t", names=["X", "Y"], skiprows=1)

df3 = pd.read_csv("district06.txt", delimiter="\t", names=["X", "Y"], skiprows=1)

# Create polygons from the extracted coordinates

polygon1 = Polygon(zip(df1["X"], df1["Y"]))

polygon2 = Polygon(zip(df2["X"], df2["Y"]))

polygon3 = Polygon(zip(df3["X"], df3["Y"]))

# Create a GeoDataFrame with 'num_coords' and 'district'

gdf = gpd.GeoDataFrame(

    {

        "district": ["01", "05", "06"],  # District identifiers

        "num_coords": [len(df1), len(df2), len(df3)],  # Number of coordinates per polygon

        "geometry": [polygon1, polygon2, polygon3]  # Polygon geometries

    },

    crs="EPSG:4326"  # Assign WGS 84 coordinate reference system

)

# Part II: Load Raster Data and Compute Agricultural Land Percentage

    # Load GLOBCOVER rasters for 2004 and 2009

    # Calculate agricultural land coverage within each district

# Define raster file paths

raster_2004 = "GLOBCOVER_2004_lab2.tif"

raster_2009 = "GLOBCOVER_2009_lab2.tif"

# Function to count total and agricultural pixels in each district

def count_pixels(raster_path, gdf):

    with rasterio.open(raster_path) as src:  # Open raster file

        affine = src.transform

        array = src.read(1)  # Read first band (land cover classification)

        total_pixel_counts = []

        agri_pixel_counts = []

        for _, row in gdf.iterrows():

            # Create a mask for the polygon

            mask = rasterio.features.geometry_mask(

                [row["geometry"]], transform=affine, out_shape=src.shape, invert=True

            )

            # Count total pixels inside the polygon

            total_pixels = np.count_nonzero(mask)

            # Count agricultural pixels (assuming agriculture is coded as "1")

            agri_pixels = np.count_nonzero((array == 1) & mask)

            total_pixel_counts.append(total_pixels)

            agri_pixel_counts.append(agri_pixels)

    return total_pixel_counts, agri_pixel_counts

# Compute total and agricultural pixels for 2004

gdf["total_pixels_2004"], gdf["agri_pixels_2004"] = count_pixels(raster_2004, gdf)

# Compute total and agricultural pixels for 2009

gdf["total_pixels_2009"], gdf["agri_pixels_2009"] = count_pixels(raster_2009, gdf)

# Calculate agricultural land percentage for each district

gdf["agri_percent_2004"] = (gdf["agri_pixels_2004"] / gdf["total_pixels_2004"]) * 100

gdf["agri_percent_2009"] = (gdf["agri_pixels_2009"] / gdf["total_pixels_2009"]) * 100

# Print pixel counts and agricultural percentages

# print("\nPixel Totals and Agricultural Pixel Totals by District\n")

# for _, row in gdf.iterrows():

#     print(f"District {row['district']}:")

#     print(f"  Total Pixels: \n    {row['total_pixels_2004']:.0f} in 2004 \n    {row['total_pixels_2009']:.0f} in 2009")

#     print(f"  Agricultural Pixels: \n    {row['agri_pixels_2004']:.0f} in 2004 \n    {row['agri_pixels_2009']:.0f} in 2009\n")

print("\nFinal Agricultural Land Percentages by District\n")

for _, row in gdf.iterrows():

    print(f"District {row['district']}:")

    print(f"  {row['agri_percent_2004']:.2f}% in 2004")

    print(f"  {row['agri_percent_2009']:.2f}% in 2009\n")