Adding basics

Author: sghelichkhani

basics_01.py

@@ -0,0 +1,190 @@
+# Tutorial: Working with Shapefiles and NetCDF in Python
+# ========
+# This tutorial introduces key geospatial analysis workflows using Python:
+# - Reading shapefiles and converting to `shapely` geometries
+# - Testing if coordinates fall inside polygons
+# - Visualising spatial data and geographic boundaries
+# - Reading and analysing gridded NetCDF rainfall datasets
+# - Creating spatial masks and integrating values over regions
+#
+# We’ll focus on the Murray-Darling Basin (MDB) as a case study.
+
+# + [markdown]
+# ### Import Required Libraries
+#
+# We begin by importing essential Python libraries:
+#
+# - `numpy` for numerical operations
+# - `netCDF4` for handling NetCDF files
+# - `matplotlib.pyplot` for plotting
+# - `pyshp` (via `shapefile.Reader`) for shapefiles
+# - `shapely.geometry` to convert and operate on vector shapes
+# + tags=["empty-cell"]
+# Import numpy, netCDF4, matplotlib.pyplot, shapefile Reader, and shapely.geometry Point and shape
+# -
+# + tags=["solution"]
+import numpy as np
+import netCDF4 as nc
+import matplotlib.pyplot as plt
+from shapefile import Reader
+from shapely.geometry import Point, shape
+# -
+# + [markdown]
+# ###  Load and Convert Shapefiles to Geometry
+# For this tutorial, we will use the MDB boundaries shapefile. Download the shapefiles from <a href="https://data.gadopt.org/water-course/MDB_boundaries.zip" download>here</a> and unzip it in the current directory.  Use `Reader()` to load two shapefiles for the MDB: one for the **north** and one for the **south**.
+# Convert each to a `shapely` geometry object so we can do spatial queries.
+#
+# This enables operations like containment tests (`Point.within(polygon)`, i.e., if a point is inside a polygon).
+#
+# + tags=["empty-cell"]
+# Load north and south MDB shapefiles and convert each to a shapely shape
+# -
+# + tags=["solution"]
+# Load the north and south MDB shapefiles and convert each to a shapely shape
+NMDB = Reader("MDB_boundaries/MDB_north_boundary.shp")
+SMDB = Reader("MDB_boundaries/MDB_south_boundary.shp")
+# Convert each to a shapely shape
+NMDB_shape = shape(NMDB.shape())
+SMDB_shape = shape(SMDB.shape())
+# -
+# + [markdown]
+# ### Check if a Point Lies Within the Basin
+#
+# Check if the point (151°E, 29°S) lies in either basin.
+# `Point().within(polygon)` returns `True` if the point lies *inside* the shape.
+# -
+# + tags=["empty-cell"]
+# Use Point.within() to test if (151, -29) is in the north or south MDB
+# -
+# + tags=["solution"]
+print(Point([151, -29]).within(NMDB_shape))
+print(Point([151, -29]).within(SMDB_shape))
+# -
+# + [markdown]
+# ### Plot Basin Boundaries and Test Coordinate
+# Extract the coordinate arrays from the shapefiles and plot them.
+# Overlay the test point to confirm visually.
+
+# + tags=["empty-cell"]
+# Extract shapefile coordinates, plot north and south boundaries, and add test point
+# -
+# + tags=["solution"]
+# Extract the coordinate arrays from the shapefiles and plot them.
+NMDB_coords = np.array(NMDB.shape().points)
+SMDB_coords = np.array(SMDB.shape().points)
+# Plot the north and south boundaries
+plt.figure(figsize=(5, 5))
+plt.plot(SMDB_coords[:, 0], SMDB_coords[:, 1])
+plt.plot(NMDB_coords[:, 0], NMDB_coords[:, 1])
+plt.scatter(151, -29, c='red')
+plt.title("North and South MDB Boundaries with Test Point")
+plt.show()
+# -
+# + [markdown]
+# ### Load and Explore Rainfall Data
+# Download the `rain_day_2025.nc`, which is the amount of rain on 27th of July 2025 as downloaded from BoM, from <a href="https://data.gadopt.org/water-course/rain_day_2025.nc" download>here</a>. This file is in NetCDF format, which is a common format for gridded data. If you want to learn more about NetCDF, you can read the <a href="https://pro.arcgis.com/en/pro-app/latest/help/data/multidimensional/what-is-netcdf-data.htm" target="_blank">NetCDF Quick Start Guide</a>.
+# Open the `rain_day_2025.nc` NetCDF file and extract:
+# - Rainfall values
+# - Latitude/longitude grids
+# - Time axis (convert to years)
+#
+# NetCDF is ideal for gridded data like daily rainfall.
+# -
+# + tags=["empty-cell"]
+# Load NetCDF file and extract rain_day, lats, lons, and converted time
+# -
+# + tags=["solution"]
+# Load the NetCDF file and extract the rainfall data, latitude, longitude, and time
+data = nc.Dataset("rain_day_2025.nc")
+# Extract the rainfall data, latitude, longitude, and time
+rain_day = data["rain_day"][:]
+# Extract the latitude and longitude data
+lats = data["latitude"][:]
+lons = data["longitude"][:]
+time = data["time"][:] / 365.25 + 1900  # convert from days since 1900
+# -
+# + [markdown]
+# ### Create Spatial Masks for Each Basin
+# We want to isolate grid cells that fall within the MDB.
+# Loop through the lat/lon grid and use `Point.within()` to assign `True` to cells inside the north or south basin.
+#
+# + tags=["empty-cell"]
+# Build boolean masks for NMDB and SMDB based on which grid points fall inside each
+# -
+# + tags=["solution"]
+# Build boolean masks for NMDB and SMDB based on which grid points fall inside each
+NMDB_mask = np.zeros((len(lats), len(lons)))
+SMDB_mask = np.zeros((len(lats), len(lons)))
+# Loop through the lat/lon grid and use `Point.within()` to assign `True` to cells inside the north or south basin.
+for ilat in range(len(lats)):
+    for ilon in range(len(lons)):
+        pt = Point([lons[ilon], lats[ilat]])
+        if pt.within(NMDB_shape):
+            NMDB_mask[ilat, ilon] = True
+        elif pt.within(SMDB_shape):
+            SMDB_mask[ilat, ilon] = True
+# -
+# + [markdown]
+# ### Plot the Mask to Verify Coverage
+# Combine the north and south masks, and plot to confirm the shape of the coverage.
+#
+# + tags=["empty-cell"]
+# Combine north/south masks and plot the result
+# -
+# + tags=["solution"]
+# Combine the north and south masks, and plot to confirm the shape of the coverage.
+MDB_mask = NMDB_mask + SMDB_mask
+plt.figure(figsize=(5, 5))
+plt.imshow(MDB_mask)
+plt.title("Spatial Mask for North + South MDB")
+plt.show()
+# -
+# + [markdown]
+# ### Visualise Rainfall for a Single Day
+# Mask the rainfall data at one time step (e.g. day 0) using the MDB mask.
+# This reveals only the rainfall within the basin.
+# -
+# + tags=["empty-cell"]
+# Apply spatial mask to rain_day[0] and plot the result
+# -
+# + tags=["solution"]
+plt.figure(figsize=(5, 5))
+plt.imshow(rain_day[0] * MDB_mask)
+plt.title("Masked Rainfall at Time = 0")
+plt.show()
+# -
+# + [markdown]
+# ### Integrate Rainfall Over Time
+# Calculate the spatially averaged rainfall over the MDB at each time step.
+# This gives a time series of total rainfall.
+# + tags=["empty-cell"]
+# Loop over time to compute daily average rainfall over MDB
+# -
+# + tags=["solution"]
+MDB_total_rain = np.zeros(len(time))
+for i in range(len(time)):
+    MDB_total_rain[i] = np.sum(MDB_mask * rain_day[i]) / np.sum(MDB_mask)
+# -
+# + [markdown]
+# ### Plot Time Series of Integrated Rainfall
+# Finally, plot the time series of total MDB rainfall to observe temporal patterns.
+# -
+# + tags=["empty-cell"]
+# Plot rainfall time series vs time
+# -
+# + tags=["solution"]
+plt.figure(figsize=(6, 2))
+# Plot the time series of total MDB rainfall
+plt.plot(time, MDB_total_rain, linewidth=1)
+# Set the x and y labels
+plt.xlabel("Year")
+plt.ylabel("Rainfall (mm)")
+plt.title("Total Rainfall in MDB over Time")
+plt.show()
+# -
+# + [markdown]
+# ### ✅ Extensions
+# To go further, try:
+# - Integrating rainfall separately for NMDB and SMDB
+# - Comparing trends or seasonal cycles between them
+# - Exporting the rainfall time series to CSV for reporting