In [ ]:
from google.colab import files
uploaded = files.upload()  # This will prompt you to upload files
Upload widget is only available when the cell has been executed in the current browser session. Please rerun this cell to enable.
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()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
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)
No description has been provided for this image
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()
No description has been provided for this image
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()
No description has been provided for this image