pandergwu / Masters-Thesis-Code

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

Issue in "KF Factor_Montecito.py" where data is transformed (mirrored and rotated 90 degrees CW) during visualization #1

Open pandergwu opened 1 month ago

pandergwu commented 1 month ago

Whenever I try to visualize the data, it appears to look correct, but as if it was mirrored and then rotated CW 90 degrees. I am very new to coding and have no idea what is causing this transformation of the data. For reference, this is what my code is producing:

image

Note the shoreline going North-South on the right side of the figure.

This is the shape of the data I am expecting to see (red square on the left):

image

Note that the coastline is East-West and not North-South.

Here is my code:

%% 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\Montecito Data\4_Mntc_KF&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()

Any help is appreciated!

alex-wzm commented 1 month ago

The response I got from copilot:

Thank you for providing the code. I've reviewed it, and I think I understand the issue your friend is experiencing. The problem is likely in the way the grid is being created and plotted. Let's focus on the visualization part in Step 6, which is where the issue seems to be occurring.

Here are a few observations and suggestions:

  1. Grid Creation: In Step 2, the grid is created correctly with UTM coordinates. The southwest corner of each cell is used as the reference point.

  2. Polygon Creation: In Step 6, the create_polygon function is defining the polygons correctly, starting from the southwest corner and moving counterclockwise.

  3. Plotting: The issue is likely in how matplotlib is interpreting and plotting the data. By default, matplotlib uses a coordinate system where the origin (0,0) is at the bottom-left corner, and the y-axis increases upwards. However, in image processing and some GIS applications, the origin is often at the top-left corner, with the y-axis increasing downwards.

To fix this issue, we can try the following:

  1. Flip the y-axis in the plot.
  2. Adjust the extent of the plot to match the data's bounds.

Here's a modified version of the plotting code in Step 6:

# ... (previous code remains the same)

# Now let's 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')

# Ensure equal aspect ratio
ax.set_aspect('equal')

# Flip the y-axis
ax.invert_yaxis()

# Set the extent of the plot to match the data's bounds
bounds = gdf.total_bounds
ax.set_xlim(bounds[0], bounds[2])
ax.set_ylim(bounds[1], bounds[3])

# Show the plot
plt.show()

These changes should correct the orientation of the plot. The ax.invert_yaxis() call flips the y-axis, and setting the extent ensures that the plot covers the correct area.

If this doesn't fully resolve the issue, we might need to investigate further:

  1. Check if the UTM coordinates are correct and consistent throughout the process.
  2. Verify that the dNBR calculations are correct and associated with the right grid cells.
  3. Consider using a different plotting library or approach, such as rasterio for raster-based visualization.

Let me know if this helps or if you need further assistance!

ai

pandergwu commented 1 month ago

Didn't work :(

image

Instead, it just added some weird white space to the right. I've been using CoPilot and ChatGPT to help me with this code, so in the past they have suggested other codes that have also created the additional space on the right. ChatGPT and CoPilot have also said that, perhaps, the issue is that the raster data (satellite data where I got the measurements to calculate the dNBR, or difference normalized burn ratio) is smaller than the study area, and so some spaces are getting cut off when I clip the data down to my study area. However, the raster data encompasses so much more of California and the study area sits squarely inside of it, so I'm not sure that is the problem.

As for your/CoPilot's response: I am using UTM as my coordinate system, which means that, like a regular Cartesian system graph, the values begin in the bottom left (i.e., [0,0]) and would increase as you move to the East (right) and North (up) along the graph. So I'm okay with transposing/transforming the data along a horizontal axis in a way where it would be mirrored, but in order for the figure to make sense, the Y-axis would still need to increase as you move north. Is there a way to artificially flip and rotate the data but keep the axes as they are? I'm like 90% sure that the data is being shown correctly and that the issue is simply with how I am trying to visualize the data with Cartopy (or whatever environment I'm using, idk).