import pandas as pd

# 1. Load the dataset
df = pd.read_csv("The_Tax_Burden_on_Tobacco__1970-2019.csv")

# 2. Filter for the desired submeasure and years 1980–1989
filtered = df[
    (df["SubMeasureDesc"] == "Cigarette Consumption (Pack Sales Per Capita)") &
    (df["Year"].between(1980, 1989))
].copy()

# 3. Convert 'Data_Value' to numeric
filtered["Data_Value"] = pd.to_numeric(filtered["Data_Value"], errors="coerce")

# 4. Compute average per state (abbreviation and name)
avg_1980s = (
    filtered.groupby(["LocationAbbr", "LocationDesc"])["Data_Value"]
    .mean()
    .reset_index()
    .rename(columns={"Data_Value": "Avg_Pack_Sales_Per_Capita_1980s"})
)

# 5. Save to a new CSV
avg_1980s.to_csv("Cigarette_Pack_Sales_Per_Capita_1980s.csv", index=False)

# 6. (Optional) View the result
print(avg_1980s.head())

  LocationAbbr LocationDesc  Avg_Pack_Sales_Per_Capita_1980s
0           AK       Alaska                           124.68
1           AL      Alabama                           115.37
2           AR     Arkansas                           125.29
3           AZ      Arizona                           109.24
4           CA   California                           104.23

import pygmt
import geopandas as gpd
import pandas as pd
import os

# --- Get current working directory ---
cwd = os.getcwd()

# --- Load cigarette dataset ---
data = pd.read_csv("Cigarette_Pack_Sales_Per_Capita_1980s.csv")

# --- Load US states shapefile ---
shapefile_path = os.path.join(cwd, "Shapefiles", "cb_2018_us_state_500k.shp")
states_gdf = gpd.read_file(shapefile_path)

# --- Merge shapefile with your data ---
states_gdf = states_gdf.merge(data, left_on="STUSPS", right_on="LocationAbbr", how="left")

# --- Create PyGMT figure ---
fig = pygmt.Figure()

region = [-125, -66.5, 24, 49.5]  # continental USA

fig.basemap(region=region, projection="M12c",
            frame=["f"])
fig.coast(land="gray90", water="lightblue",
          shorelines="1/0.5p,black", borders=["1/0.5p,black"])

# --- Create colormap ---
pygmt.makecpt(
    cmap="amp",
    series=[data["Avg_Pack_Sales_Per_Capita_1980s"].min(),
            data["Avg_Pack_Sales_Per_Capita_1980s"].max()],
    continuous=True
)

# --- Plot each state colored by your data ---
fig.plot(
    data=states_gdf,
    pen="0.3p,gray10",
    fill="+z",  # fill polygons by z-value
    cmap=True,
    aspatial="Z=Avg_Pack_Sales_Per_Capita_1980s"
)

# --- Add colorbar ---
fig.colorbar(frame=["x+lAverage Cigarette Pack Sales per Capita (1980s) (Lower 48)"])

# --- Show figure ---
#fig.show()

#fig.savefig("cigarette_map.png")

makecpt [WARNING]: Without inc in -T option, -Z has no effect (ignored)