In [ ]:
from google.colab import files
uploaded = files.upload() # This will prompt you to upload files
Saving CWA-EQ-Catalog-1990.xml to CWA-EQ-Catalog-1990 (1).xml Saving CWA-EQ-Catalog-1991.xml to CWA-EQ-Catalog-1991.xml Saving CWA-EQ-Catalog-1992.xml to CWA-EQ-Catalog-1992.xml Saving CWA-EQ-Catalog-1993.xml to CWA-EQ-Catalog-1993.xml Saving CWA-EQ-Catalog-1994.xml to CWA-EQ-Catalog-1994.xml Saving CWA-EQ-Catalog-1995.xml to CWA-EQ-Catalog-1995.xml Saving CWA-EQ-Catalog-1996.xml to CWA-EQ-Catalog-1996.xml Saving CWA-EQ-Catalog-1997.xml to CWA-EQ-Catalog-1997.xml Saving CWA-EQ-Catalog-1998.xml to CWA-EQ-Catalog-1998.xml Saving CWA-EQ-Catalog-1999.xml to CWA-EQ-Catalog-1999.xml Saving CWA-EQ-Catalog-2000.xml to CWA-EQ-Catalog-2000.xml Saving CWA-EQ-Catalog-2001.xml to CWA-EQ-Catalog-2001.xml Saving CWA-EQ-Catalog-2002.xml to CWA-EQ-Catalog-2002.xml Saving CWA-EQ-Catalog-2003.xml to CWA-EQ-Catalog-2003.xml Saving CWA-EQ-Catalog-2004.xml to CWA-EQ-Catalog-2004.xml Saving CWA-EQ-Catalog-2005.xml to CWA-EQ-Catalog-2005.xml Saving CWA-EQ-Catalog-2006.xml to CWA-EQ-Catalog-2006.xml Saving CWA-EQ-Catalog-2007.xml to CWA-EQ-Catalog-2007.xml Saving CWA-EQ-Catalog-2008.xml to CWA-EQ-Catalog-2008.xml Saving CWA-EQ-Catalog-2009.xml to CWA-EQ-Catalog-2009.xml Saving CWA-EQ-Catalog-2010.xml to CWA-EQ-Catalog-2010.xml Saving CWA-EQ-Catalog-2011.xml to CWA-EQ-Catalog-2011.xml Saving CWA-EQ-Catalog-2012.xml to CWA-EQ-Catalog-2012.xml Saving CWA-EQ-Catalog-2013.xml to CWA-EQ-Catalog-2013.xml Saving CWA-EQ-Catalog-2014.xml to CWA-EQ-Catalog-2014.xml Saving CWA-EQ-Catalog-2015.xml to CWA-EQ-Catalog-2015.xml Saving CWA-EQ-Catalog-2016.xml to CWA-EQ-Catalog-2016.xml Saving CWA-EQ-Catalog-2017.xml to CWA-EQ-Catalog-2017.xml Saving CWA-EQ-Catalog-2018.xml to CWA-EQ-Catalog-2018.xml Saving CWA-EQ-Catalog-2019.xml to CWA-EQ-Catalog-2019.xml Saving CWA-EQ-Catalog-2020.xml to CWA-EQ-Catalog-2020.xml Saving CWA-EQ-Catalog-2021.xml to CWA-EQ-Catalog-2021.xml Saving CWA-EQ-Catalog-2022.xml to CWA-EQ-Catalog-2022.xml Saving CWA-EQ-Catalog-2023.xml to CWA-EQ-Catalog-2023.xml Saving CWA-EQ-Catalog-2024.xml to CWA-EQ-Catalog-2024.xml
In [76]:
import pandas as pd
import xml.etree.ElementTree as ET
import os
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Point
import contextily as ctx
import numpy as np
years = np.arange(1990,2025)
# Define the region of interest (adjust these based on your region)
xlim = [119, 125] # Longitude range for Taiwan
ylim = [20, 27] # Latitude range for Taiwan
for year in years:
# Replace this with the path to your XML file
xml_file = f'CWA-EQ-Catalog-{year}.xml'
# Parse the XML file
tree = ET.parse(xml_file)
root = tree.getroot()
# Define the namespace to be used in the search
namespaces = {'ns': 'urn:cwa:gov:tw:cwacommon:0.1'}
# List to hold epicenter data
epicenters = []
# Iterate through all earthquake info in the XML
for earthquake_info in root.findall('.//ns:EarthquakeInfo', namespaces):
# Extract latitude and longitude of the epicenter
latitude = float(earthquake_info.find('ns:EpicenterLatitude', namespaces).text)
longitude = float(earthquake_info.find('ns:EpicenterLongitude', namespaces).text)
local_magnitude = float(earthquake_info.find('ns:LocalMagnitude', namespaces).text.strip())
# Append the data to the epicenters list
epicenters.append({'latitude': latitude, 'longitude': longitude, 'LocalMagnitude': local_magnitude})
# Create a DataFrame from the extracted epicenter data
df = pd.DataFrame(epicenters)
# Show the DataFrame
output_dir = 'earthquake_images'
# Categorize the LocalMagnitude into bins
bins = [3, 4, 5, 6, 7]
labels = ['3-4', '4-5', '5-6', '6-7']
df['MagnitudeCategory'] = pd.cut(df['LocalMagnitude'], bins=bins, labels=labels, right=False)
# Create a GeoDataFrame from the DataFrame
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
gdf['size'] = gdf['LocalMagnitude'] ** 3
# Check if CRS is None and set it to WGS84 if needed
if gdf.crs is None:
gdf = gdf.set_crs('EPSG:4326')
# Define a colormap for the categories
color_map = {
'3-4': 'blue',
'4-5': 'green',
'5-6': 'orange',
'6-7': 'red'
}
# Define alpha values for each category
alpha_map = {
'3-4': 0.1,
'4-5': 0.3,
'5-6': 0.5,
'6-7': 0.7
}
# Map the colors and alphas based on the category
gdf['color'] = gdf['MagnitudeCategory'].map(color_map)
gdf['alpha'] = gdf['MagnitudeCategory'].map(alpha_map)
# Remove rows with NaN colors (invalid categories)
gdf = gdf.dropna(subset=['color'])
# Plot the map
fig, ax = plt.subplots(figsize=(8,8))
gdf.plot(ax=ax, color=gdf['color'], alpha=gdf['alpha'], edgecolor='k', markersize=gdf['size'])
# Set xlim and ylim
ax.set_xlim(xlim)
ax.set_ylim(ylim)
# Add basemap (e.g., CartoDB Positron)
ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.CartoDB.Positron)
# Add a legend
import matplotlib.patches as mpatches
legend_labels = [mpatches.Patch(color=color_map[label], label=label) for label in labels]
legend = ax.legend(handles=legend_labels, title="2024\nLocal Magnitude", title_fontsize=14)
# Center align the legend title by accessing its text object
legend_title = legend.get_title()
legend_title.set_text(f"{year}\nLocal Magnitude")
legend_title.set_ha("center") # Align center horizontally
# Display the plot
ax.set_title("Earthquake Epicenters Categorized by Magnitude")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
output_path = os.path.join(output_dir, f"earthquake_{year}.png")
plt.tight_layout()
plt.savefig(output_path)
plt.show()
In [77]:
# Create a DataFrame from the extracted epicenter data
df = pd.DataFrame(epicenters)
# Create a GeoDataFrame from the DataFrame
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
# Check if CRS is None
if gdf.crs is None:
# Set the CRS, for example, to WGS84
gdf = gdf.set_crs('EPSG:4326')
# Your GeoDataFrame (example: 'gdf')
ax = gdf.plot(figsize=(10, 10), alpha=0.5, edgecolor='k')
# Add CartoDB Positron basemap
ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.CartoDB.Positron)
In [78]:
# Create a DataFrame from the extracted epicenter data
df = pd.DataFrame(epicenters)
# Create a GeoDataFrame from the DataFrame
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
# Check if CRS is None
if gdf.crs is None:
# Set the CRS, for example, to WGS84
gdf = gdf.set_crs('EPSG:4326')
# Scale marker size based on the LocalMagnitude (e.g., multiply by a factor)
gdf['size'] = gdf['LocalMagnitude'] ** 3 # Adjust multiplier as needed
# Plot the GeoDataFrame with marker sizes based on LocalMagnitude
fig, ax = plt.subplots(figsize=(10, 10))
gdf.plot(ax=ax, alpha=0.5, edgecolor='none', markersize=gdf['size'], color='red')
# Add CartoDB Positron basemap
ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.CartoDB.Positron)
# Labels and Title
ax.set_title('Epicenter Locations with Magnitude-based Marker Sizes')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
# Show the plot
plt.tight_layout()
plt.show()
In [79]:
output_dir = 'earthquake_images'
# Categorize the LocalMagnitude into bins
bins = [3, 4, 5, 6, 7]
labels = ['3-4', '4-5', '5-6', '6-7']
df['MagnitudeCategory'] = pd.cut(df['LocalMagnitude'], bins=bins, labels=labels, right=False)
# Create a GeoDataFrame from the DataFrame
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.longitude, df.latitude))
gdf['size'] = gdf['LocalMagnitude'] ** 3
# Check if CRS is None and set it to WGS84 if needed
if gdf.crs is None:
gdf = gdf.set_crs('EPSG:4326')
# Define a colormap for the categories
color_map = {
'3-4': 'blue',
'4-5': 'green',
'5-6': 'orange',
'6-7': 'red'
}
# Define alpha values for each category
alpha_map = {
'3-4': 0.1,
'4-5': 0.3,
'5-6': 0.5,
'6-7': 0.7
}
# Map the colors and alphas based on the category
gdf['color'] = gdf['MagnitudeCategory'].map(color_map)
gdf['alpha'] = gdf['MagnitudeCategory'].map(alpha_map)
# Plot the map
fig, ax = plt.subplots(figsize=(8,8))
gdf.plot(ax=ax, color=gdf['color'], alpha=gdf['alpha'], edgecolor='k', markersize=gdf['size'])
# Add basemap (e.g., CartoDB Positron)
ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.CartoDB.Positron)
# Add a legend
import matplotlib.patches as mpatches
legend_labels = [mpatches.Patch(color=color_map[label], label=label) for label in labels]
legend = ax.legend(handles=legend_labels, title="2024\nLocal Magnitude", title_fontsize=14)
# Center align the legend title by accessing its text object
legend_title = legend.get_title()
legend_title.set_text("2024\nLocal Magnitude")
legend_title.set_ha("center") # Align center horizontally
# Display the plot
ax.set_title("Earthquake Epicenters Categorized by Magnitude")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
output_path = os.path.join(output_dir, f"earthquake_{year}.png")
plt.tight_layout()
plt.savefig(output_path)
plt.show()