Brayden Noh

Here is a GIF of the change in Zonal District Map in Pasadena.

var region = ee.Geometry.Polygon([
  [
    [-118.22, 34.08],
    [-118.04, 34.08],
    [-118.04, 34.28],
    [-118.22, 34.28],
    [-118.22, 34.08]
  ]
]);

var maskClouds = function(image) {
  var qa = image.select('pixel_qa');
  var cloud = qa.bitwiseAnd(1 << 5).or(qa.bitwiseAnd(1 << 3));
  var mask = cloud.not();
  return image.updateMask(mask);
};

var calculateNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  return image.addBands(ndvi);
};

var years = [2000,2010,2020];

years.forEach(function(year) {
  var startDate = year + '-08-01'; // Start date in August
  var endDate = year + '-08-20';   // End date in August

  var l7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
            .filterBounds(region)
            .filterDate(startDate, endDate)
            .map(maskClouds)
            .map(calculateNDVI);

  var meanNDVI = l7.select('NDVI').mean();

  Export.image.toDrive({
    image: meanNDVI,
    description: 'NDVI_L7_August_' + year,
    region: region,
    scale: 30,
    maxPixels: 1e10
  });
});

import matplotlib.pyplot as plt
import rasterio
from matplotlib.colors import Normalize, LinearSegmentedColormap
import geopandas as gpd
from rasterstats import zonal_stats
import numpy as np

colors = ["#640000", "#ff0000", "#ffff00", "#00c800", "#006400"]
cmap = LinearSegmentedColormap.from_list("ndvi", colors, N=256)

tiff_path = '/Users/braydennoh/Downloads/drive-download-20240514T154154Z-001/NDVI_L7_August_2000.tif'
satellite_photo_path = '/Users/braydennoh/Downloads/drive-download-20240514T154154Z-001/satellite_photo.tif'
shapefile_path = '/Users/braydennoh/Downloads/pasadenandvi/Zoning_-4568161463289048110/Zoning.shp'

with rasterio.open(tiff_path) as src:
    ndvi = src.read(1, masked=True)
    affine = src.transform
    bounds = src.bounds
    extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]
    width = src.width
    height = src.height

with rasterio.open(satellite_photo_path) as src_sat:
    satellite_photo = src_sat.read([1, 2, 3])
    satellite_affine = src_sat.transform

zones = gpd.read_file(shapefile_path)
if zones.crs != src.crs:
    zones = zones.to_crs(src.crs)

zonal_stats_results = zonal_stats(zones, ndvi, affine=affine, stats="mean")
zones['mean_ndvi'] = [stat['mean'] for stat in zonal_stats_results]

fig_width = 10
fig_height = fig_width * (height / width)
fig, ax = plt.subplots(figsize=(fig_width, fig_height))

ax.imshow(satellite_photo.transpose(1, 2, 0), extent=extent)
cax = ax.imshow(ndvi, cmap=cmap, norm=Normalize(vmin=0, vmax=0.5), extent=extent, alpha=0.6)
zones.boundary.plot(ax=ax, edgecolor='white', linewidth=1)
zones.plot(ax=ax, column='mean_ndvi', cmap='RdBu', alpha=0.6)

cbar = fig.colorbar(cax, ax=ax, orientation='vertical')
cbar.set_label('NDVI')
ax.set_title('Average NDVI in Pasadena Zonal Districts - 2000')

x_ticks = np.arange(-118.22, -118.04, 0.04)
y_ticks = np.arange(34.08, 34.28, 0.02)
ax.set_xticks(x_ticks)
ax.set_yticks(y_ticks)
ax.set_xlim([-118.21, -118.06])
ax.set_ylim([34.11, 34.26])
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
plt.axis("off")
plt.show()