Air Quality Data Analysis¶

# Load Data

# Assuming 'air_quality.csv' is in the same directory as the notebook
try:
    df = pd.read_csv('air_quality.csv')
    print("CSV loaded successfully.")
    # Display basic info and first few rows
    print(df.info())
    display(df.head())
except FileNotFoundError:
    print("Error: 'air_quality.csv' not found. Please ensure the file is in the correct directory.")
except Exception as e:
    print(f"An error occurred while loading the CSV: {e}")

# Data Cleaning and Preparation

if 'df' in locals(): # Check if df was loaded successfully
    # Clean column names (remove leading/trailing spaces, replace spaces/hyphens with underscores)
    df.columns = df.columns.str.strip().str.replace(' ', '_').str.replace('-', '_')

    # Convert Start_Date to datetime
    df['Start_Date'] = pd.to_datetime(df['Start_Date'])

    # Extract Year based on the Time_Period column for annual averages
    # For Winter YYYY-YY, use the latter year. For Summer YYYY, use YYYY. For Annual YYYY, use YYYY.
    def get_analysis_year(row):
        period = str(row['Time_Period'])
        start_year = row['Start_Date'].year
        if 'Winter' in period and '-' in period:
            try:
                return int(period.split('-')[1]) + 2000 # Assumes YY format like 14-15
            except:
                 return start_year + 1 # Fallback if format is unexpected
        elif 'Annual Average' in period:
             # If start date is Dec, it's likely for the *next* year's average
             if row['Start_Date'].month == 12:
                 return start_year + 1
             else:
                 return start_year
        else: # Summer or other cases where start date year is usually correct
            return start_year

    df['Analysis_Year'] = df.apply(get_analysis_year, axis=1)

    # Extract Season from Time Period
    def get_season(time_period):
        period_str = str(time_period).lower()
        if 'winter' in period_str:
            return 'Winter'
        elif 'summer' in period_str:
            return 'Summer'
        elif 'annual' in period_str:
            return 'Annual'
        else:
            return 'Other'

    df['Season'] = df['Time_Period'].apply(get_season)

    # Convert Data_Value to numeric, coercing errors
    df['Data_Value'] = pd.to_numeric(df['Data_Value'], errors='coerce')
    # Drop rows where Data_Value could not be converted
    original_rows = len(df)
    df.dropna(subset=['Data_Value'], inplace=True)
    if original_rows > len(df):
        print(f"Dropped {original_rows - len(df)} rows with non-numeric Data_Value.")

    # Display cleaned info and sample
    print("\nCleaned DataFrame Info:")
    print(df.info())
    display(df[['Start_Date', 'Analysis_Year', 'Season', 'Time_Period', 'Data_Value']].head())
else:
    print("DataFrame 'df' not available for cleaning.")

# Overview of Indicators

if 'df' in locals():
    indicator_counts = df['Name'].value_counts()
    print("Indicators Measured:")
    print(indicator_counts)

    plt.figure(figsize=(10, max(6, len(indicator_counts) * 0.5))) # Adjust height based on number of indicators
    sns.barplot(y=indicator_counts.index, x=indicator_counts.values, palette='viridis', orient='h')
    plt.title('Number of Records per Indicator')
    plt.xlabel('Number of Records')
    plt.ylabel('Indicator Name')
    plt.tight_layout()
    plt.show()
else:
    print("DataFrame 'df' not available for indicator overview.")

# Filter for NO2 and PM2.5

if 'df' in locals():
    df_pollutants = df[df['Name'].isin(['Nitrogen dioxide (NO2)', 'Fine particles (PM 2.5)'])].copy()
    print(f"Filtered DataFrame shape for NO2 and PM2.5: {df_pollutants.shape}")
    display(df_pollutants.head())
else:
     print("DataFrame 'df' not available for filtering.")

# Temporal Trends (Annual Average)

if 'df_pollutants' in locals():
    annual_avg = df_pollutants[df_pollutants['Season'] == 'Annual']

    if not annual_avg.empty:
        # Calculate mean and std dev across locations for each year/pollutant
        annual_trend = annual_avg.groupby(['Analysis_Year', 'Name'])['Data_Value'].agg(['mean', 'std']).reset_index()

        plt.figure(figsize=(14, 7))
        for pollutant in annual_trend['Name'].unique():
            pollutant_data = annual_trend[annual_trend['Name'] == pollutant]
            plt.errorbar(pollutant_data['Analysis_Year'], pollutant_data['mean'], yerr=pollutant_data['std'], 
                         label=pollutant, marker='o', capsize=5)

        plt.title('NYC Annual Average Pollutant Levels Over Time (Mean ± SD across locations)')
        plt.xlabel('Year')
        plt.ylabel('Concentration (Units vary by pollutant)')
        plt.legend(title='Pollutant')
        plt.xticks(rotation=45)
        # Ensure x-axis shows integer years if possible
        plt.gca().xaxis.set_major_locator(plt.MaxNLocator(integer=True))
        plt.tight_layout()
        plt.show()
        print("Insight: General decreasing trend observed for both pollutants.")
    else:
        print("No 'Annual' season data found for NO2 or PM2.5.")
else:
    print("Filtered DataFrame 'df_pollutants' not available.")

# Seasonal Variation

if 'df_pollutants' in locals():
    seasonal_data = df_pollutants[df_pollutants['Season'].isin(['Winter', 'Summer'])]

    if not seasonal_data.empty:
        plt.figure(figsize=(12, 7))
        sns.boxplot(data=seasonal_data, x='Name', y='Data_Value', hue='Season', palette='coolwarm')
        plt.title('Seasonal Variation of NO2 and PM2.5 Levels')
        plt.xlabel('Pollutant')
        plt.ylabel('Concentration (Units vary)')
        # Consider log scale if distributions are highly skewed
        # plt.yscale('log') 
        plt.legend(title='Season')
        plt.tight_layout()
        plt.show()
        print("Insight: Concentrations are generally higher in Winter than in Summer.")
    else:
        print("No 'Winter' or 'Summer' season data found for NO2 or PM2.5.")
else:
    print("Filtered DataFrame 'df_pollutants' not available.")

# Geographical Variation (Community Districts - Annual Average for a recent year)

if 'df_pollutants' in locals() and 'annual_avg' in locals() and not annual_avg.empty:
    # Find the most recent year with 'Annual' data
    recent_annual_year = annual_avg['Analysis_Year'].max()
    geo_data_recent = annual_avg[(annual_avg['Analysis_Year'] == recent_annual_year) & (annual_avg['Geo_Type_Name'] == 'CD')]

    if not geo_data_recent.empty:
        # Separate plots for NO2 and PM2.5 due to different scales
        pollutants_to_plot = geo_data_recent['Name'].unique()
        num_pollutants = len(pollutants_to_plot)
        
        fig, axes = plt.subplots(num_pollutants, 1, figsize=(15, 8 * num_pollutants), sharex=False)
        if num_pollutants == 1:
            axes = [axes] # Make it iterable if only one pollutant
            
        fig.suptitle(f'Average Annual Pollutant Levels by Community District ({int(recent_annual_year)})', fontsize=16, y=1.02)

        for i, pollutant in enumerate(pollutants_to_plot):
            ax = axes[i]
            data_subset = geo_data_recent[geo_data_recent['Name'] == pollutant].sort_values('Data_Value', ascending=False)
            
            sns.barplot(data=data_subset, y='Geo_Place_Name', x='Data_Value', ax=ax, palette='coolwarm', orient='h')
            ax.set_title(pollutant)
            ax.set_xlabel(f'Concentration ({data_subset["Measure_Info"].iloc[0] if not data_subset.empty else ""})')
            ax.set_ylabel('Community District')
            
        plt.tight_layout()
        plt.show()
        print("Insight: Shows variation in pollutant levels between different Community Districts.")
    else:
        print(f"No annual Community District (CD) data found for the most recent year ({int(recent_annual_year)}).")
elif 'df_pollutants' not in locals():
     print("Filtered DataFrame 'df_pollutants' not available.")
else:
    print("No annual average data available for geographical analysis.")