import warnings
warnings.filterwarnings("ignore")
# This suppresses warning messages to keep our output clean.

# Let us install GeoPandas if it is not already here.
import sys
try:
    import geopandas
except ImportError:
    !{sys.executable} -m pip install geopandas
    

# Let us import the main tools.
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt

# Data setup
# Let us download a US States shapefile from a simple URL.
# Source: Census Bureau. This is public data.
url = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_500k.zip"
import os
shpzip = "cb_2018_us_state_500k.zip"
if not os.path.exists(shpzip):
    import urllib.request
    urllib.request.urlretrieve(url, shpzip)
gdf = gpd.read_file("zip://" + shpzip)
print("Shape:", gdf.shape)
gdf.head()

# The GeoDataFrame is like a spreadsheet, but one column is 'geometry'.
print("Columns:", list(gdf.columns))
print("Geometry column sample:")
print(gdf['geometry'].head(2))

# Let us plot the map using gdf.plot().
gdf.plot(figsize=(10,6), edgecolor="black", color="lightblue")
plt.title("US States Map")
plt.show()

# Show some western US states only.
western = ["CA", "WA", "OR", "NV", "AZ", "NM"]
gdf_west = gdf[gdf['STUSPS'].isin(western)]
gdf_west.plot(figsize=(8,5), edgecolor="black", color="orange")
plt.title("Western US States")
plt.show()

# Let us label the states on the map.
ax = gdf_west.plot(figsize=(8,5), color="wheat", edgecolor="black")
for x, y, label in zip(gdf_west.geometry.centroid.x, gdf_west.geometry.centroid.y, gdf_west['NAME']):
    ax.text(x, y, label, fontsize=9)
plt.title("Western States with Labels")
plt.show()

# Practice: Try entering your own abbreviation.
abbr = input("Enter a US state abbreviation (e.g. CA): ")
if abbr in list(gdf['STUSPS']):
    state_row = gdf[gdf['STUSPS'] == abbr]
    state_row.plot(figsize=(5,5), color="lightgreen", edgecolor="black")
    plt.title(f"Map of {abbr}")
    plt.show()
else:
    print("Sorry, that state is not found.")
    

# Making new columns: let us compute the state area in square kilometers.
gdf['area_km2'] = gdf.to_crs(epsg=3395).geometry.area / 1e6
gdf[['STUSPS', 'NAME', 'area_km2']].head()

# Plot the top 5 biggest states by area.
largest5 = gdf.nlargest(5, 'area_km2')
ax = gdf.plot(color='lightgrey', edgecolor='grey', figsize=(10,6))
largest5.plot(ax=ax, color='darkblue', edgecolor='white')
for x, y, abbr in zip(largest5.geometry.centroid.x, largest5.geometry.centroid.y, largest5['STUSPS']):
    ax.text(x, y, abbr, fontsize=11, color='white')
plt.title('Top 5 Largest US States')
plt.show()

# Save the filtered western states as a new shapefile.
gdf_west.to_file("western_states.shp")
print("Shapefile saved as western_states.shp")

# Download USGS earthquake data (past week, all M1+).
quakes_url = "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.csv"
quakes = pd.read_csv(quakes_url)
print("Quake count:", quakes.shape[0])
quakes[['place','mag','latitude','longitude']].head()

# Convert earthquake locations to GeoDataFrame points.
from shapely.geometry import Point
quake_gdf = gpd.GeoDataFrame(
    quakes, geometry=[Point(xy) for xy in zip(quakes.longitude, quakes.latitude)], crs='EPSG:4326')
# Plot earthquakes and state boundaries together.
ax = gdf.plot(figsize=(10,6), color='white', edgecolor='black')
quake_gdf.plot(ax=ax, markersize=quake_gdf['mag']*3, alpha=0.5, color='red')
plt.title('US Map: Recent Earthquakes')
plt.show()

# How many earthquakes in California this week?
from shapely.geometry import Polygon
ca_shape = gdf[gdf['STUSPS']=='CA'].geometry.values[0]
quakes_in_ca = quake_gdf[quake_gdf.within(ca_shape)]
print(f"Number of earthquakes in CA: {len(quakes_in_ca)}")

# Clean up files created earlier (optional).
import glob
for ext in ["shp","shx","dbf","prj","cpg"]:
    for fname in glob.glob(f"western_states.{ext}"):
        os.remove(fname)