{ "cells": [ { "cell_type": "markdown", "id": "3664680a-c7f2-404e-99b6-50705a75c69e", "metadata": {}, "source": [ "# MAIAC AOD Extraction — Step by Step\n", "\n", "**Required libraries & environment** \n", "- Python 3.x \n", "- `pandas` (table I/O), `numpy` (numeric ops) \n", "- `pyhdf` (read HDF4; provides `from pyhdf.SD import SD, SDC`) \n", "- stdlib: `os`, `glob`, `datetime`, `time`, `warnings`\n", "\n", "Install (Conda recommended on Derecho):\n", "```bash\n", "conda install -c conda-forge pandas numpy pyhdf\n", "```\n", "\n", "**Data prerequisites** \n", "- Precomputed **site→grid index map** with MAIAC indices per site: \n", " `/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/epa_amtic/benzene_sites_satellite_indices.csv`\n", "- Benzene dataset: \n", " `/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/epa_amtic/benzene_hcho_amtic_regional.csv`\n", "- MAIAC AOD (MCD19A2CMG) daily HDFs stored as `/maiac_aod/YYYY/DDD/MCD19A2CMG.AYYYYDDD*.hdf`\n", "\n", "---\n", "\n", "## Extraction workflow\n", "\n", "1. **Load MAIAC indices** \n", " Build a dictionary keyed by `SITEID` with `MAIAC_LAT_INDEX`, `MAIAC_LON_INDEX`, `MAIAC_DISTANCE_KM`.\n", "\n", "2. **Load benzene dataset** \n", " Add an empty `AOD` column, parse `DATE`, and collect unique dates.\n", "\n", "3. **Find MAIAC files** \n", " For each date, compute `YYYY` and `DDD` (Julian day) and look under `/maiac_aod/YYYY/DDD/` for a valid `MCD19A2CMG.AYYYYDDD*.hdf` file.\n", "\n", "4. **Read the HDF file** \n", " Using `pyhdf.SD`, open the file and read the `AOD_055` dataset. Also read attributes: `scale_factor`, `add_offset`, `_FillValue`, and `valid_range`.\n", "\n", "5. **Extract site-level AOD** \n", " For each site on that date: \n", " - Use precomputed indices to slice `aod_data[lat_idx, lon_idx]`. \n", " - QC filter: ignore fill values or values outside `valid_range`. \n", " - Scale: `scaled = raw * scale_factor + add_offset`. \n", " - Keep only values in [0, 5]; else set `NaN`.\n", "\n", "6. **Update dataframe** \n", " Assign extracted values into the `AOD` column of benzene records.\n", "\n", "7. **Progress reporting** \n", " Print updates every ~50 dates: records/sec, dates/min, files found, successful extractions, ETA.\n", "\n", "8. **Finalize** \n", " Save enriched dataset to: \n", " `/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/epa_amtic/benzene_hcho_amtic_regional_with_aod.csv` \n", " Print summary statistics and sample results.\n", "\n", "---\n", "\n", "## Why this approach is efficient\n", "- **Index-based lookup** avoids per-record geolocation searches. \n", "- **Date-wise batching** ensures each daily HDF file is opened only once. \n", "- **Robust QC & scaling** ensures only valid physical AOD values are included. \n", "- **Graceful handling of missing files** allows the pipeline to skip dates without crashing.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d479f7dc-ca9a-41fa-b51d-18128f209c35", "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import os\n", "import glob\n", "from datetime import datetime\n", "import time\n", "import warnings\n", "\n", "warnings.filterwarnings('ignore')\n", "\n", "def load_maiac_indices():\n", " \"\"\"Load MAIAC indices from mapping file\"\"\"\n", " \n", " indices_file = '/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/epa_amtic/benzene_sites_satellite_indices.csv'\n", " \n", " print(f\"Loading indices from: {indices_file}\")\n", " indices_df = pd.read_csv(indices_file)\n", " print(f\"Loaded indices for {len(indices_df)} sites\")\n", " \n", " # Convert to dictionary for fast lookup\n", " indices_dict = {}\n", " for _, row in indices_df.iterrows():\n", " indices_dict[row['SITEID']] = {\n", " 'maiac_lat_idx': row['MAIAC_LAT_INDEX'],\n", " 'maiac_lon_idx': row['MAIAC_LON_INDEX'],\n", " 'maiac_distance': row['MAIAC_DISTANCE_KM']\n", " }\n", " \n", " return indices_dict\n", "\n", "def find_valid_maiac_file(maiac_base, target_date):\n", " \"\"\"Find valid MAIAC AOD file for specific date\"\"\"\n", " \n", " date_obj = datetime.strptime(target_date, '%Y-%m-%d')\n", " year = date_obj.year\n", " day_of_year = date_obj.timetuple().tm_yday\n", " \n", " year_dir = os.path.join(maiac_base, str(year))\n", " day_dir = os.path.join(year_dir, f\"{day_of_year:03d}\")\n", " \n", " if not os.path.exists(day_dir):\n", " return None\n", " \n", " # Find valid HDF files\n", " all_files = os.listdir(day_dir)\n", " date_code = f\"A{year}{day_of_year:03d}\"\n", " \n", " for filename in all_files:\n", " if (filename.startswith(f'MCD19A2CMG.{date_code}') and \n", " filename.endswith('.hdf') and \n", " '?details' not in filename):\n", " \n", " full_path = os.path.join(day_dir, filename)\n", " if os.path.getsize(full_path) > 1000:\n", " return full_path\n", " \n", " return None\n", "\n", "def extract_aod_for_date(aod_file, site_records, indices_dict):\n", " \"\"\"Extract AOD for all sites on a specific date\"\"\"\n", " \n", " try:\n", " from pyhdf.SD import SD, SDC\n", " \n", " hdf = SD(aod_file, SDC.READ)\n", " sds = hdf.select('AOD_055')\n", " \n", " # Get data and scaling parameters\n", " aod_data = sds.get()\n", " attrs = sds.attributes()\n", " \n", " scale_factor = attrs.get('scale_factor', 0.001)\n", " add_offset = attrs.get('add_offset', 0.0)\n", " fill_value = attrs.get('_FillValue', -28672)\n", " valid_range = attrs.get('valid_range', [0, 6000])\n", " \n", " results = {}\n", " \n", " # Process each site record for this date\n", " for idx, row in site_records.iterrows():\n", " site_id = row['SITEID']\n", " \n", " if site_id not in indices_dict:\n", " results[idx] = np.nan\n", " continue\n", " \n", " lat_idx = indices_dict[site_id]['maiac_lat_idx']\n", " lon_idx = indices_dict[site_id]['maiac_lon_idx']\n", " \n", " if lat_idx < 0 or lon_idx < 0:\n", " results[idx] = np.nan\n", " continue\n", " \n", " # Extract and scale AOD value\n", " raw_value = aod_data[lat_idx, lon_idx]\n", " \n", " if (raw_value == fill_value or \n", " raw_value < valid_range[0] or \n", " raw_value > valid_range[1]):\n", " results[idx] = np.nan\n", " else:\n", " scaled_aod = (raw_value * scale_factor) + add_offset\n", " \n", " # Final validity check\n", " if scaled_aod < 0 or scaled_aod > 5:\n", " results[idx] = np.nan\n", " else:\n", " results[idx] = float(scaled_aod)\n", " \n", " sds.endaccess()\n", " hdf.end()\n", " \n", " return results\n", " \n", " except Exception as e:\n", " print(f\" Error extracting AOD from {os.path.basename(aod_file)}: {e}\")\n", " return {}\n", "\n", "def main():\n", " \"\"\"Main extraction function\"\"\"\n", " \n", " print(\"=\"*60)\n", " print(\"MAIAC AOD EXTRACTION FOR BENZENE DATASET\")\n", " print(\"=\"*60)\n", " print(f\"Start time: {datetime.now()}\")\n", " print()\n", " \n", " # File paths\n", " maiac_base = '/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/satellite/maiac_aod/'\n", " benzene_file = '/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/epa_amtic/benzene_hcho_amtic_regional.csv'\n", " output_file = '/glade/derecho/scratch/shahriar/research_projects/benzene_prediction/01_raw_data/epa_amtic/benzene_hcho_amtic_regional_with_aod.csv'\n", " \n", " # Load data\n", " print(\"LOADING DATA\")\n", " print(\"-\" * 20)\n", " indices_dict = load_maiac_indices()\n", " if not indices_dict:\n", " print(\"ERROR: Could not load MAIAC indices\")\n", " return None\n", " \n", " print(f\"Loading benzene data from: {benzene_file}\")\n", " df_benzene = pd.read_csv(benzene_file)\n", " print(f\"Loaded {len(df_benzene):,} benzene measurements\")\n", " print()\n", " \n", " # Initialize AOD column\n", " df_benzene['AOD'] = np.nan\n", " \n", " # Process by unique dates\n", " df_benzene['DATE_PARSED'] = pd.to_datetime(df_benzene['DATE'])\n", " unique_dates = df_benzene['DATE_PARSED'].dt.strftime('%Y-%m-%d').unique()\n", " unique_dates = sorted(unique_dates)\n", " \n", " print(\"PROCESSING DATES\")\n", " print(\"-\" * 20)\n", " print(f\"Total unique dates: {len(unique_dates)}\")\n", " print(f\"Date range: {unique_dates[0]} to {unique_dates[-1]}\")\n", " print()\n", " \n", " processed_records = 0\n", " successful_extractions = 0\n", " files_found = 0\n", " dates_processed = 0\n", " \n", " start_time = time.time()\n", " \n", " for i, date_str in enumerate(unique_dates):\n", " # Skip dates before MODIS Terra launch\n", " if date_str < '2000-02-24':\n", " continue\n", " \n", " dates_processed += 1\n", " \n", " # Progress reporting every 50 dates\n", " if dates_processed % 50 == 0 or dates_processed == 1:\n", " elapsed = time.time() - start_time\n", " if dates_processed > 1:\n", " rate = processed_records / elapsed\n", " dates_per_min = (dates_processed / elapsed) * 60\n", " remaining_dates = len(unique_dates) - i - 1\n", " eta = remaining_dates / dates_per_min if dates_per_min > 0 else 0\n", " \n", " print(f\"PROGRESS UPDATE:\")\n", " print(f\" Date {dates_processed}: {date_str}\")\n", " print(f\" Processed {i+1:,}/{len(unique_dates):,} dates ({(i+1)/len(unique_dates)*100:.1f}%)\")\n", " print(f\" Records/sec: {rate:.1f} | Dates/min: {dates_per_min:.1f}\")\n", " print(f\" Files found: {files_found} | Successful extractions: {successful_extractions:,}\")\n", " print(f\" ETA: {eta:.1f} minutes\")\n", " print()\n", " \n", " # Get records for this date\n", " date_mask = df_benzene['DATE_PARSED'].dt.strftime('%Y-%m-%d') == date_str\n", " date_records = df_benzene[date_mask]\n", " \n", " # Find MAIAC file for this date\n", " aod_file = find_valid_maiac_file(maiac_base, date_str)\n", " \n", " if not aod_file:\n", " processed_records += len(date_records)\n", " continue\n", " \n", " files_found += 1\n", " \n", " # Extract AOD for all sites on this date\n", " aod_results = extract_aod_for_date(aod_file, date_records, indices_dict)\n", " \n", " # Store results\n", " for idx, aod_value in aod_results.items():\n", " df_benzene.loc[idx, 'AOD'] = aod_value\n", " if not np.isnan(aod_value):\n", " successful_extractions += 1\n", " \n", " processed_records += len(date_records)\n", " \n", " # Final timing and statistics\n", " total_time = time.time() - start_time\n", " \n", " print(\"=\"*60)\n", " print(\"EXTRACTION COMPLETE!\")\n", " print(\"=\"*60)\n", " print(f\"End time: {datetime.now()}\")\n", " print(f\"Total runtime: {total_time/60:.1f} minutes ({total_time/3600:.2f} hours)\")\n", " print(f\"Processing rate: {processed_records/total_time:.1f} records/second\")\n", " print()\n", " \n", " # Results summary\n", " print(\"EXTRACTION SUMMARY:\")\n", " print(\"-\" * 20)\n", " print(f\"Total records processed: {processed_records:,}\")\n", " print(f\"Dates with MAIAC files: {files_found:,}\")\n", " print(f\"Successful AOD extractions: {successful_extractions:,}\")\n", " \n", " success_rate = successful_extractions / processed_records * 100 if processed_records > 0 else 0\n", " print(f\"Overall success rate: {success_rate:.1f}%\")\n", " print()\n", " \n", " # AOD statistics\n", " valid_aod = df_benzene['AOD'].dropna()\n", " if len(valid_aod) > 0:\n", " print(\"AOD STATISTICS:\")\n", " print(\"-\" * 15)\n", " print(f\"Valid AOD values: {len(valid_aod):,}\")\n", " print(f\"AOD range: {valid_aod.min():.4f} to {valid_aod.max():.4f}\")\n", " print(f\"Mean AOD: {valid_aod.mean():.4f}\")\n", " print(f\"Median AOD: {valid_aod.median():.4f}\")\n", " print(f\"Standard deviation: {valid_aod.std():.4f}\")\n", " print()\n", " \n", " # Clean up temporary column\n", " df_benzene = df_benzene.drop('DATE_PARSED', axis=1)\n", " \n", " # Save results\n", " print(\"SAVING RESULTS:\")\n", " print(\"-\" * 15)\n", " print(f\"Output file: {output_file}\")\n", " df_benzene.to_csv(output_file, index=False)\n", " print(f\"Results saved successfully!\")\n", " \n", " # Verify saved file\n", " file_size_mb = os.path.getsize(output_file) / (1024 * 1024)\n", " print(f\"Output file size: {file_size_mb:.1f} MB\")\n", " print()\n", " \n", " # Show sample results\n", " print(\"SAMPLE RESULTS:\")\n", " print(\"-\" * 15)\n", " sample_data = df_benzene[df_benzene['AOD'].notna()].head(10)\n", " print(sample_data[['SITEID', 'STATE', 'DATE', 'BENZENE', 'HCHO', 'AOD']].to_string(index=False))\n", " print()\n", " \n", " print(\"=\"*60)\n", " print(\"EXTRACTION COMPLETED SUCCESSFULLY!\")\n", " print(\"=\"*60)\n", " \n", " return df_benzene\n", "\n", "if __name__ == \"__main__\":\n", " result = main()" ] }, { "cell_type": "code", "execution_count": null, "id": "8f6c33e7-112c-4358-920f-e952e69e80d4", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.13" } }, "nbformat": 4, "nbformat_minor": 5 }