Geospatial Python

# Numeric
integer = 42
float_num = 3.14
cloud_cover = 15

# Strings
text = "Hello EO"
multiline = """Multiple
lines"""

# Lists (mutable)
coords = [14.5995, 120.9842]  # Manila lat/lon
sensors = ["Sentinel-1", "Sentinel-2", "Landsat-8"]

# Tuples (immutable)
bbox = (120.0, 14.0, 121.0, 15.0)

# Dictionaries
metadata = {
    "satellite": "Sentinel-2",
    "date": "2025-01-15",
    "cloud_cover": 15
}

# If statements
if cloud_cover < 20:
    print("Good quality image")
elif cloud_cover < 50:
    print("Moderate quality")
else:
    print("Too cloudy")

# For loops
for sensor in sensors:
    print(f"Processing {sensor} data")

# List comprehension
times_two = [x * 2 for x in range(11)]
print(times_two)

# While loop
count = 1
while count <= 10:
    print(count)
    count += 1

class Dog:
    def __init__(self, name, age):
        self.name = name
        self.age = age

    def bark(self):
        print(f"{self.name} says Woof!")

    def wag_tail(self):
        print(f"{self.name} is wagging their tail!")

my_dog = Dog("Buddy", 3)
print(my_dog.name)  # prints "Buddy"
my_dog.bark()  # prints "Buddy says Woof!"
my_dog.wag_tail()  # prints "Buddy is wagging their tail!"

class Image:

    def __init__(self, metadata):
        self.metadata = metadata

img1 = Image(metadata)
img2 = Image({
        "satellite": "Sentinel-2",
        "date": "2025-01-16",
        "cloud_cover": 32
        })

images = [img1, img2]

for img in images:
    print(img.metadata)

import numpy as np

# From lists
arr = np.array([1, 2, 3, 4, 5])
matrix = np.array([[1, 2], [3, 4]])

# Special arrays
zeros = np.zeros((3, 3))           # 3x3 array of zeros
ones = np.ones((2, 4))             # 2x4 array of ones
identity = np.eye(4)               # 4x4 identity matrix
random = np.random.rand(3, 3)      # 3x3 random [0, 1)
range_arr = np.arange(0, 10, 2)    # [0, 2, 4, 6, 8]
linspace = np.linspace(0, 1, 5)    # 5 values from 0 to 1

import numpy as np

# Define the arrays
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print("Original Arrays:")
print("a:", a)
print("b:", b)

# Arithmetic Operations
print("\nArithmetic Operations:")

# Sum
sum = a + b
print("Sum:", sum)

# Difference
diff = a - b
print("Difference:", diff)

# Product
prod = a * b
print("Product:", prod)

# Quotient
quotient = a / b
print("Quotient:", quotient)

# Power
power = a ** 2
print("Power:", power)

arr = np.array([1, 2, 3, 4, 5])
print("\nStatistics Operations:")

# Mean
mean = arr.mean()
print("Mean:", mean)

# Standard Deviation
std_dev = arr.std()
print("Standard Deviation:", std_dev)

# Minimum Value
min_val = arr.min()
print("Minimum Value:", min_val)

# Maximum Value
max_val = arr.max()
print("Maximum Value:", max_val)

# Sum of Elements
sum_val = arr.sum()
print("Sum:", sum_val)

# Indexing Operations
matrix = np.array([[1, 2], [3, 4]])
print("\nIndexing Operations:")

# First Row
first_row = matrix[0, :]
print("First Row:", first_row)

# Second Column
second_col = matrix[:, 1]
print("Second Column:", second_col)

# Boolean indexing
print("\nBoolean Indexing Operation:")
arr = np.array([1, 2, 3, 4, 5])
greater_than_5 = arr[arr > 5]
print("Elements Greater than 5:", greater_than_5)

# Shape Operations
arr = np.array([[1, 2, 3], [4, 5, 6]])

print("Original Array:\n", arr)

# Reshape to 3x2
reshaped_arr = arr.reshape(3, 2)
print("Reshaped Array (3x2):\n", reshaped_arr)

# Flatten to 1D
flattened_arr = arr.flatten()
print("Flattened Array:\n", flattened_arr)

# Transpose
transposed_arr = arr.transpose()
print("Transposed Array:\n", transposed_arr)

# Stacking Operations
# Vertical Stack
vertical_stack = np.vstack([a, b])
print("Vertical Stack:\n", vertical_stack)

# Horizontal Stack
horizontal_stack = np.hstack([a, b])
print("Horizontal Stack:\n", horizontal_stack)

# Stack along axis
stacked_arr = np.stack([a, b], axis=0)
print("Stacked Array:\n", stacked_arr)

# Concatenation Operations
# Concatenate
concatenated_arr = np.concatenate([a, b])
print("Concatenated Array:\n", concatenated_arr)

import geopandas as gpd

gdf_municities_gpkg = gpd.read_file("data/phl.gpkg", layer='ncr_municities_pop') # Vector data inside of a geopackage
gdf_municities_gpkg

df_mcdo_csv = gpd.read_file("data/NCR_McDonalds.csv") # CSV w/ a geometry field but with fields for latitude and longitude

df_mcdo_csv

# Let's add a geometry field using the known latitudes and longitudes
gdf_mcdo_csv = gpd.GeoDataFrame(df_mcdo_csv, geometry=gpd.points_from_xy(df_mcdo_csv.longitude, df_mcdo_csv.latitude))

gdf_mcdo_csv

import geopandas as gpd

gdf_municities = gpd.read_file("data/phl.gpkg", layer='ncr_municities_pop')

gdf_municities.crs # Print current CRS

gdf_municities_3123 = gdf_municities.to_crs(3123) # Create new GDF that's a copy of gdf_municities but w/ CRS as EPSG:3123

gdf_municities_3123 # Notice the coordinates used in the geometry field

# Filtering by attributes
import geopandas as gpd

gdf_municities = gpd.read_file("data/phl.gpkg", layer='ncr_municities_pop')

gdf_municities[gdf_municities.pop2020 > 1000000] # Get municities with pop2020 > 1000000

dist2 = gdf_municities[gdf_municities["province"] == "Second District, NCR"] # Get only the municities in the Second District of NCR
dist2

increase = gdf_municities[(gdf_municities.pop2020 - gdf_municities.pop2015) > 50000] # Get only the municities whose populatin increased by at least 50000 from 2015 to 2020
increase

# Add the population from 1990 to 2010 to the data

import pandas as pd

ncr_pop1990_2010 = pd.read_csv("data/ncr_pop1990_2010.csv") # dataframe of 1990-2010 population
ncr_pop1990_2010 

gdf_municities_pop1990_2020 = gdf_municities_3123.merge(ncr_pop1990_2010, on="psgc_municity") #join the data using the psgc_municity field 

gdf_municities_pop1990_2020

gdf_municities_pop1990_2020_fix1 = gdf_municities_pop1990_2020.copy() # IMPORTANT: It's always a good idea to make a copy of your dataframes
gdf_municities_pop1990_2020_fix1 = gdf_municities_pop1990_2020_fix1.drop("municity_y", axis=1) # drop the municity_y
gdf_municities_pop1990_2020_fix1.rename(columns={'municity_x':'municity'}, inplace=True)
gdf_municities_pop1990_2020_fix1 

gdf_municities_pop1990_2020_fix2 = gdf_municities_pop1990_2020.copy() # IMPORTANT: It's always a good idea to make a copy of your dataframes
del gdf_municities_pop1990_2020_fix2["municity_y"] # delete the municity_y
gdf_municities_pop1990_2020_fix2.rename(columns={'municity_x':'municity'}, inplace=True)
gdf_municities_pop1990_2020_fix2 

gdf_municities_3123.area # Area in sqm

# Try it here
# gdf_municities.area # Area in sqm

gdf_municities_3123.area/10000 # Area in hectares 

gdf_ncr_floods = gpd.read_file("data/ncr.gpkg", layer="manila_flood_hazard_lipad")
gdf_ncr_hospitals = gpd.read_file("data/ncr.gpkg", layer="manila_hospitals_osm")

gdf_ncr_hospitals

hospitals_join = gdf_ncr_hospitals.sjoin(gdf_ncr_floods, how="left") # Tag/assign the Var and hazard fields for each hospital
hospitals_join

gdf_municities.plot("pop2020", legend="True") # Simple map using GeoPandas

gdf_municities.explore("pop2020", legend="True") # Interactive version

hospitals_join.explore("hazard", legend="True")

import rasterio

ncr_dem = rasterio.open("data/ncr_srtm_dem.tif")

ncr_dem.meta

ncr_dem.shape # (rows, cols)

ncr_dem.count # number of bands

ncr_dem.dtypes # data type of pixel value per band

ncr_dem.bounds # bounding box in data CRS

%matplotlib inline
from rasterio import plot
plot.show(ncr_dem)

%matplotlib inline
plot.show_hist(ncr_dem)

from rasterio.plot import show

with rasterio.open("data/pan.tiff") as src:
    show(src, cmap="gray")

with rasterio.open("data/multi.tiff") as src:
    show((src, [3, 2, 1]))  # True color (R, G, B)

up_multiband = rasterio.open("data/multi.tiff")
up_pan = rasterio.open("data/pan.tiff")

from matplotlib import pyplot as plt

fig, (axr, axg, axb) = plt.subplots(1,3, figsize=(21,7))
plot.show((up_multiband, 1), ax=axr, cmap='Reds', title='R')
plot.show((up_multiband, 2), ax=axg, cmap='Greens', title='G')
plot.show((up_multiband, 3), ax=axb, cmap='Blues', title='B')

import numpy as np

rasterpath = 'data/sample_landsat.tif'

satellite_image = rasterio.open(rasterpath)

# Read Band 3 (Red) and Band 4 (NIR)
red = satellite_image.read(3).astype('float32')
nir = satellite_image.read(4).astype('float32')

# Create an empty array of zeros that is the exact same size as our image
ndvi = np.zeros_like(red)

# Create a True/False map of all the pixels where (nir + red) is NOT zero
valid_pixels = (nir + red) != 0

# Only calculate NDVI on those valid pixels 
# The background pixels will just stay as the 0s in step 2.
ndvi[valid_pixels] = (nir[valid_pixels] - red[valid_pixels]) / (nir[valid_pixels] + red[valid_pixels])

print(f"Max NDVI (Healthiest vegetation): {np.max(ndvi):.2f}")
print(f"Min NDVI (Water or Built-up): {np.min(ndvi):.2f}")

import matplotlib.pyplot as plt

# Initialize the plot
fig, ax = plt.subplots(figsize=(10, 8))

# Plot the NDVI array
# cmap='RdYlGn' stands for Red-Yellow-Green.
# vmin and vmax lock our color scale strictly between -1 and 1
im = ax.imshow(ndvi, cmap='RdYlGn', vmin=-1, vmax=1)

# Add a colorbar to the side so we know what the colors mean
plt.colorbar(im, ax=ax, label="NDVI Value")

plt.title("Landsat NDVI (Crop Health Index)", fontsize=16)
ax.set_axis_off() # Hides the X/Y pixel numbers on the borders
plt.show()

import geopandas as gpd
from shapely.geometry import box
import rasterio.mask
import matplotlib.pyplot as plt

# Get the boundaries of original satellite image
bounds = satellite_image.bounds
minx, miny, maxx, maxy = bounds.left, bounds.bottom, bounds.right, bounds.top

# Create a smaller bounding box in the center of the image
center_x = (minx + maxx) / 2
center_y = (miny + maxy) / 2
aoi_geom = box(center_x - 5000, center_y - 5000, center_x + 5000, center_y + 5000)

# Convert Shapely shape into a GeoDataFrame (so it has a CRS)
aoi_gdf = gpd.GeoDataFrame({'geometry': [aoi_geom]}, crs=satellite_image.crs)

# Perform the Clip
clipped_image_array, clipped_transform = rasterio.mask.mask(
    satellite_image, 
    aoi_gdf.geometry.values, 
    crop=True
)

print(f"Original shape: {satellite_image.shape}")
print(f"Clipped shape: {clipped_image_array.shape[1:]}")

# Visualize the clipped result
fig, ax = plt.subplots(figsize=(6, 6))
ax.imshow(clipped_image_array[2], cmap='gray') # Band 3 (Red) is at index 2
plt.title("Clipped Satellite Image (Red Band)")
ax.set_axis_off()
plt.show()