49 lines
1.8 KiB
Python
49 lines
1.8 KiB
Python
import geopandas as gpd
|
|
import rasterio
|
|
from rasterio.mask import mask
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from matplotlib.ticker import ScalarFormatter
|
|
|
|
# Read the original shapefile
|
|
gdf = gpd.read_file('continent_overlay/Continent_Outlines.shp')
|
|
|
|
# Path to the raster file
|
|
raster_path = 'gtopo_world_dem_6m.tif'
|
|
|
|
# Open the raster file
|
|
with rasterio.open(raster_path) as src:
|
|
plt.figure(figsize=(14, 8)) # Create a larger figure for better visualization
|
|
|
|
# Loop through each unique continent
|
|
for continent in gdf['CONTINENT'].unique():
|
|
# Filter GeoDataFrame for the current continent
|
|
gdf_continent = gdf[gdf['CONTINENT'] == continent]
|
|
|
|
# Clip the raster with the geometry of the current continent
|
|
out_image, out_transform = mask(src, gdf_continent.geometry, crop=True)
|
|
|
|
# Flatten the array and remove masked values
|
|
data = out_image[0].flatten() # Assuming a single-band raster
|
|
data = data[(data != src.nodata) & (data >= 0) & (data <= 7000)] # Remove nodata values and filter the range
|
|
|
|
# Calculate histogram data (bin counts)
|
|
counts, bin_edges = np.histogram(data, bins=50, range=(0, 7000))
|
|
|
|
# Plot the line for this continent
|
|
plt.plot(bin_edges[:-1], counts, label=continent, linewidth=1.5)
|
|
|
|
# Add x-axis and y-axis limits, legend, title, and labels
|
|
plt.xlim(0, 7000)
|
|
plt.ylim(0, 40000)
|
|
plt.gca().yaxis.set_major_formatter(ScalarFormatter(useOffset=False)) # Ensure y-axis shows whole numbers
|
|
plt.ticklabel_format(style='plain', axis='y') # Disable scientific notation
|
|
plt.legend(loc='upper right')
|
|
plt.xlabel('Elevation (m)')
|
|
plt.ylabel('Count')
|
|
|
|
# Save the plot to a file
|
|
plt.savefig('plot1.png')
|
|
|
|
print('Line plot saved as "plot1.png"')
|