pandergwu / Masters-Thesis-Code

Code to visualize dNBR data for Montecito landslides
0 stars 0 forks source link

Issue in "KF Factors_LCF.py" - same issue as in the Montecito case, where the data is being mirrored and rotated. #2

Open pandergwu opened 2 weeks ago

pandergwu commented 2 weeks ago

For reference, this is what I'm getting (rotated and mirrored to make it easier to match the real data): Figure 2024-10-10 124339

When this is what I am expecting to see (I traced the general outline to help the viewer identify the burned areas I'm looking at) : image

The code is included in the .py file mentioned in the title, specifically Step 6, which I'm also including below:

%% Step 6

Visualize the dNBR

import pandas as pd import geopandas as gpd import matplotlib.pyplot as plt from shapely.geometry import Polygon from pyproj import CRS

Load the Excel file with dNBR data

file_path = "D:\Masters\Thesis Data\LCF Data\3_LCF_dNBR.xlsx" df = pd.read_excel(file_path)

Print the min and max UTM X and Y coordinates

print("UTM_X range:", df['UTM_X'].min(), "-", df['UTM_X'].max()) print("UTM_Y range:", df['UTM_Y'].min(), "-", df['UTM_Y'].max())

Ensure that the DataFrame has the correct CRS (EPSG:32611 for UTM Zone 11N)

crs = CRS.from_epsg(32611)

Function to create a polygon for each grid cell based on the SW corner

def create_polygon(utm_x, utm_y, grid_size=30): return Polygon([ (utm_x, utm_y), (utm_x, utm_y + grid_size), (utm_x + grid_size, utm_y + grid_size), (utm_x + grid_size, utm_y), (utm_x, utm_y) ])

Convert the DataFrame to a GeoDataFrame with polygon geometries

df['geometry'] = df.apply(lambda row: create_polygon(row['UTM_X'], row['UTM_Y']), axis=1) gdf = gpd.GeoDataFrame(df, geometry='geometry', crs=crs)

Print a few sample geometries to verify correctness

print(gdf.geometry.head())

Plot the data using GeoPandas without Cartopy

fig, ax = plt.subplots(figsize=(12, 12))

Plot the geometry and dNBR values

gdf.plot(column='dNBR', cmap='RdYlBu', ax=ax, edgecolor='none')

Set the title for the plot

ax.set_title('dNBR Values without Cartopy')

Ensure equal aspect ratio

ax.set_aspect('equal')

Show the plot

plt.show()