task9.txt

Task: US Drought Analysis using Standardized Precipitation Index (SPI)

Produce a drought severity map for the contiguous United States for October 1, 2024, by calculating the Standardized Precipitation Index (SPI). Your objective is to use a long-term climate baseline to assess drought conditions and create a clear, publication-quality visualization and a brief descriptive report.

---
Data Download and Acquisition
---

CDS API Configuration:
- Use the exact dataset name: `"reanalysis-era5-land-monthly-means"`
- Variable: `["total_precipitation"]` (use list format)
- Domain: `"area": [50, -125, 25, -65]` (North, West, South, East format)
- Years: Include all years from 1950 to 2024 as strings: `["1950", "1951", ..., "2024"]`
- Months: Include all months as zero-padded strings: `["01", "02", ..., "12"]`
- Time: `["00:00"]`
- Product type: `["monthly_averaged_reanalysis"]`
- **Critical**: Set `"data_format": "grib"` and `"download_format": "unarchived"`

Example CDS request structure:
```python
dataset = "reanalysis-era5-land-monthly-means"
request = {
    "product_type": ["monthly_averaged_reanalysis"],
    "variable": ["total_precipitation"],
    "year": ["1950", "1951", ..., "2024"],
    "month": ["01", "02", ..., "12"],
    "time": ["00:00"],
    "data_format": "grib",
    "download_format": "unarchived",
    "area": [50, -125, 25, -65]
}
```

---
Data Loading Specifications
---

GRIB File Loading:
- **Critical Filter**: Use `filter_by_keys={'stepType': 'avgas'}` for monthly averaged data
- **Do NOT use**: `typeOfLevel` or `shortName` filters - these cause performance issues
- **Engine**: Use `engine="cfgrib"` with xarray
- **Variable Name**: The precipitation variable will be named `'tp'` in the loaded dataset

Correct loading pattern:
```python
ds = xr.open_dataset(
    "filename.grib", 
    engine="cfgrib", 
    filter_by_keys={'stepType': 'avgas'}
)
```

Domain Cropping:
- Use `.sel()` method with slice notation
- Latitude: `slice(50, 25)` (Note: ERA5 latitudes are ordered high to low)
- Longitude: `slice(-125, -65)`

```python
us_ds = ds.sel(
    latitude=slice(50, 25),
    longitude=slice(-125, -65)
)
```

---
SPI Calculation Implementation
---

Data Separation:
- Climatology: `slice('1995-01-01', '2023-12-31')` (29 years: 1995-2023)
- Target data: `slice('2024-10-01', '2024-10-01')` (October 2024 only)

Monthly Statistics Calculation:
```python
# Calculate climatology statistics grouped by month
monthly_mean = climatology.groupby('time.month').mean(dim='time')
monthly_std = climatology.groupby('time.month').std(dim='time')

# For October (month=10) analysis
target = target_year['tp'].sel(time=target_year['time.month'] == 10)
mean = monthly_mean['tp'].sel(month=10)
std = monthly_std['tp'].sel(month=10)
```

SPI Calculation with Error Handling:
```python
# Handle zero standard deviation
std = std.where(std > 0, np.nan)

# Calculate SPI with safe division
with np.errstate(divide='ignore', invalid='ignore'):
    spi = (target - mean) / std

# Handle infinite values
spi = spi.where(~np.isinf(spi), np.nan)
```

---
Visualization Requirements
---

Drought Categories and Colors:
- Define exact SPI thresholds: `[-2.5, -2, -1.5, -1, -0.5, 0]`
- Category labels: `['D4', 'D3', 'D2', 'D1', 'D0']`
- **Exact colors**: `['#730000', '#E60000', '#FFAA00', '#FCD37F', '#FFFF00', 'white']`

Plotting Strategy:
1. **Two-layer approach**: Plot non-drought areas (SPI ≥ 0) in white separately from drought areas
2. Use `ListedColormap` and `BoundaryNorm` for discrete categories
3. **Critical**: Set `add_colorbar=False` in plot calls, create custom colorbar

Map Configuration:
```python
# Create figure and axis
fig = plt.figure(figsize=(9, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# Plot non-drought areas first
non_drought = spi_data.where(spi_data >= 0)
non_drought.plot(ax=ax, transform=ccrs.PlateCarree(), 
                cmap=ListedColormap(['white']), add_colorbar=False)

# Plot drought areas with categories
drought = spi_data.where(spi_data < 0)
drought.plot(ax=ax, transform=ccrs.PlateCarree(),
            cmap=ListedColormap(colors[:-1]),  # Exclude white
            norm=BoundaryNorm(categories[:-1], len(colors)-1),
            add_colorbar=False)
```

Map Features and Styling:
- Projection: `ccrs.PlateCarree()`
- Required features: `COASTLINE`, `BORDERS`, `STATES`
- Extent: `[-125, -65, 25, 50]`
- Title format: `'US Drought Monitor Categories (SPI-based) for 2024-10-01\nClimatology: 1995-2023'`

Custom Colorbar:
```python
cbar = plt.colorbar(
    plot, ax=ax, label='Drought Category',
    ticks=[-2.25, -1.75, -1.25, -0.75, -0.25],  # Center of each category
    extend='neither', spacing='proportional'
)
cbar.set_ticklabels(['D4', 'D3', 'D2', 'D1', 'D0'])
```

---
Code Structure Requirements
---

Import Statements:
```python
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import ListedColormap, BoundaryNorm
import cdsapi
```

Error Handling:
- Include proper exception handling for file operations
- Handle division by zero in SPI calculations
- Use `np.errstate()` context manager for numerical warnings

Performance Considerations:
- Use the correct `stepType: 'avgas'` filter to avoid hanging
- Avoid unnecessary indexpath configurations
- Process data in logical chunks (climatology vs target period)

Output Requirements:
- Save plots with high DPI (300) and tight bounding box
- Use descriptive filenames with dates
- Include progress logging for long operations

---
Common Pitfalls to Avoid
---

1. **Wrong GRIB filters**: Do NOT use `'typeOfLevel': 'surface', 'shortName': 'tp'`
2. **Incorrect variable names**: ERA5-Land uses `'tp'`, not `'total_precipitation'`
3. **Indexpath issues**: Avoid setting `indexpath: None` - use empty string or omit
4. **Colorbar complexity**: Don't try to include non-drought areas in drought colorbar
5. **Date formatting**: Use proper datetime slicing with xarray
6. **Memory issues**: Separate climatology and target data processing

This detailed specification should prevent the performance and functionality issues encountered in the original implementation.