Selena’s Test!

Imports:

import fsspec
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import pyplot as plt
from matplotlib.dates import DateFormatter, AutoDateLocator,HourLocator,DayLocator
import numpy as np
import xarray as xr
from metpy.plots import USCOUNTIES
import pandas as pd
from datetime import datetime, timedelta
import matplotlib.colors as colors
from pyproj import Proj
import seaborn as sns
sns.set()

Load in hourly Precip Data and calculate total for storm duration.

zarr_url = f's3://mdmf/gdp/stageiv_combined.zarr/'

fs2 = fsspec.filesystem('s3', anon=True, endpoint_url='https://usgs.osn.mghpcc.org/')

ds = xr.open_dataset(fs2.get_mapper(zarr_url), engine='zarr',

                             backend_kwargs={'consolidated':True}, chunks={})

ds = ds.sortby('time')

# Extract the relevant time range
start_time = '2006-09-16 12:00:00'
end_time = '2006-09-16 12:00:00'

precipitation_data = ds['Total_precipitation_surface_1_Hour_Accumulation'].sel(time=slice(start_time, end_time))

# Sum the precipitation values over the specified time range
total_precipitation = precipitation_data.sum(dim='time')

# Convert total_precipitation to inches
#conversion_factor = 0.0393701
#total_precipitation = total_precipitation * conversion_factor

np.max(total_precipitation)

<xarray.DataArray 'Total_precipitation_surface_1_Hour_Accumulation' ()> Size: 4B
dask.array<_nanmax_skip-aggregate, shape=(), dtype=float32, chunksize=(), chunktype=numpy.ndarray>

xarray.DataArray

'Total_precipitation_surface_1_Hour_Accumulation'

• dask.array<chunksize=(), meta=np.ndarray>

  Array Chunk Bytes 4 B 4 B Shape () () Dask graph 1 chunks in 14 graph layers Data type float32 numpy.ndarray

• Coordinates: (0)
• Indexes: (0)
• Attributes: (0)

Plot total precip for entire CONUS.

ds

<xarray.Dataset> Size: 758GB
Dimensions:                                          (time: 191948, y: 881,
                                                      x: 1121, time_bounds_1: 2)
Coordinates:
    lat                                              (x, y) float32 4MB dask.array<chunksize=(102, 80), meta=np.ndarray>
    lon                                              (x, y) float32 4MB dask.array<chunksize=(102, 80), meta=np.ndarray>
  * time                                             (time) datetime64[ns] 2MB ...
Dimensions without coordinates: y, x, time_bounds_1
Data variables:
    Total_precipitation_surface_1_Hour_Accumulation  (time, y, x) float32 758GB dask.array<chunksize=(1597, 80, 102), meta=np.ndarray>
    crs                                              int32 4B ...
    time_bounds                                      (time, time_bounds_1) datetime64[ns] 3MB dask.array<chunksize=(47987, 1), meta=np.ndarray>
Attributes: (12/34)
    Conventions:                          CF-1.4
    Generating_process_or_model:          River Forecast Center Quantitative ...
    Metadata_Conventions:                 Unidata Dataset Discovery v1.0
    Originating_or_generating_Center:     US National Weather Service, Nation...
    Originating_or_generating_Subcenter:  Environmental Modeling Center
    _CoordSysBuilder:                     ucar.nc2.dataset.conv.CF1Convention
    ...                                   ...
    processing_level:                     Aggregated .grib archive from the N...
    project:                              Various
    summary:                              Mosaicked into a national product a...
    time_coverage_end:                    present
    time_coverage_start:                  2002-01-01T00:00:00Z
    title:                                United States Stage IV Quantitative...

xarray.Dataset

• Dimensions:
  • time: 191948
  • y: 881
  • x: 1121
  • time_bounds_1: 2
• Coordinates: (3)
  • lat
    (x, y)
    float32
    dask.array<chunksize=(102, 80), meta=np.ndarray>

    units :

    degrees_north

    Array Chunk Bytes 3.77 MiB 31.88 kiB Shape (1121, 881) (102, 80) Dask graph 132 chunks in 2 graph layers Data type float32 numpy.ndarray

  • lon
    (x, y)
    float32
    dask.array<chunksize=(102, 80), meta=np.ndarray>

    units :

    degrees_east

    Array Chunk Bytes 3.77 MiB 31.88 kiB Shape (1121, 881) (102, 80) Dask graph 132 chunks in 2 graph layers Data type float32 numpy.ndarray

  • time
    (time)
    datetime64[ns]
    2002-01-01 ... 2023-11-30T23:00:00

    bounds :

    time_bounds

    long_name :

    GRIB forecast or observation time

    standard_name :

    time

    array(['2002-01-01T00:00:00.000000000', '2002-01-01T01:00:00.000000000',
           '2002-01-01T02:00:00.000000000', ..., '2023-11-30T21:00:00.000000000',
           '2023-11-30T22:00:00.000000000', '2023-11-30T23:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (3)
  • Total_precipitation_surface_1_Hour_Accumulation
    (time, y, x)
    float32
    dask.array<chunksize=(1597, 80, 102), meta=np.ndarray>

    cell_methods :

    time: sum (interval: 1 hr)

    description :

    Total precipitation

    grid_mapping :

    crs: lat lon

    long_name :

    Total precipitation (1_Hour Accumulation) @ Ground or water surface

    standard_name :

    precipitation_amount

    units :

    kg m^-2

    Array Chunk Bytes 706.20 GiB 49.71 MiB Shape (191948, 881, 1121) (1597, 80, 102) Dask graph 29964 chunks in 3 graph layers Data type float32 numpy.ndarray

  • crs
    ()
    int32
    ...

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    geographic_crs_name :

    WGS 84

    grid_mapping_name :

    latitude_longitude

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    inverse_flattening :

    298.257223563

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    reference_ellipsoid_name :

    WGS 84

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    [1 values with dtype=int32]

  • time_bounds
    (time, time_bounds_1)
    datetime64[ns]
    dask.array<chunksize=(47987, 1), meta=np.ndarray>

    long_name :

    bounds for time

    Array Chunk Bytes 2.93 MiB 374.90 kiB Shape (191948, 2) (47987, 1) Dask graph 8 chunks in 3 graph layers Data type datetime64[ns] numpy.ndarray

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2002-01-01 00:00:00', '2002-01-01 01:00:00',
                   '2002-01-01 02:00:00', '2002-01-01 03:00:00',
                   '2002-01-01 04:00:00', '2002-01-01 05:00:00',
                   '2002-01-01 06:00:00', '2002-01-01 07:00:00',
                   '2002-01-01 08:00:00', '2002-01-01 09:00:00',
                   ...
                   '2023-11-30 14:00:00', '2023-11-30 15:00:00',
                   '2023-11-30 16:00:00', '2023-11-30 17:00:00',
                   '2023-11-30 18:00:00', '2023-11-30 19:00:00',
                   '2023-11-30 20:00:00', '2023-11-30 21:00:00',
                   '2023-11-30 22:00:00', '2023-11-30 23:00:00'],
                  dtype='datetime64[ns]', name='time', length=191948, freq=None))

• Attributes: (34)

  Conventions :

  CF-1.4

  Generating_process_or_model :

  River Forecast Center Quantitative Precipitation estimate mosaic generated by NCEP

  Metadata_Conventions :

  Unidata Dataset Discovery v1.0

  Originating_or_generating_Center :

  US National Weather Service, National Centres for Environmental Prediction (NCEP)

  Originating_or_generating_Subcenter :

  Environmental Modeling Center

  _CoordSysBuilder :

  ucar.nc2.dataset.conv.CF1Convention

  acknowledgement :

  Data originated at the National Center for Environmental Prediction: http://www.emc.ncep.noaa.gov/mmb/ylin/pcpanl/stage4/

  cdm_data_type :

  Grid

  contributor_name :

  Center for Integrated Data Analytics

  contributor_role :

  Aggregation and Serving via THREDDS

  date_created :

  ongoing

  featureType :

  GRID

  geospatial_lat_max :

  53

  geospatial_lat_min :

  24

  geospatial_lat_resolution :

  4

  geospatial_lat_units :

  km

  geospatial_lon_max :

  -66

  geospatial_lon_min :

  -125

  geospatial_lon_resolution :

  4

  geospatial_lon_units :

  km

  geospatial_vertical_units :

  km

  history :

  Aggregated from NCEP Stage IV Analysis by dblodgett@usgs.gov 2016-01-01

  id :

  stageiv

  institution :

  National Center for Environmental Prediction

  keywords :

  Atmosphere > Precipitation > Precipitation Amount Atmosphere > Precipitation > Rain

  keywords_vocabulary :

  GCMD Science Keywords

  license :

  Freely available

  naming_authority :

  cida.usgs.gov

  processing_level :

  Aggregated .grib archive from the NCEP. The original Hydrologic Rainfall Analysis Projection grid cell locations redefined by unprojecting grid cell locations using spherical earth projection algorithms. The cell lat/lon locations should be interpreted as being on a modern geotetic datum. This transformation is needed because the original grid cell location's lat/lon values were on a geodetic datum but were assumed to be on a spherical datum. See: https://pubs.usgs.gov/fs/2013/3035/pdf/fs2013-3035.pdf for more information. Note that this publication reverences the National Precipitation Verification Unit. The data represented here are aggregated from the same source as the National Precipitation Verification Unit used prior to being shut down. Data processing notes can be found here: http://www.emc.ncep.noaa.gov/mmb/ylin/pcpanl/

  project :

  Various

  summary :

  Mosaicked into a national product at NCEP, from the regional hourly multi-sensor (radar+gauges) precipitation analyses produced by the 12 River Forecast Centers over CONUS. Some manual QC done at the RFCs. The Stage II/IV job is run at 33min past the top of each hour. Hourly Stage IV is re-made hourly (if there is new input after valid time for the next 23 hours, then again at 1/3/5/7 days after valid time. Source: http://www.emc.ncep.noaa.gov/mmb/ylin/pcpanl/stage4/ Note: This dataset is not truely national. Some RFCs do not produce hourly precip estimates and others may not be present in the dataset. Please review the NCEP documentation for more information.

  time_coverage_end :

  present

  time_coverage_start :

  2002-01-01T00:00:00Z

  title :

  United States Stage IV Quantitative Precipitation Archive

lon = ds['lon'].values
lat = ds['lat'].values

# Transpose the total_precipitation array if needed
total_precipitation = total_precipitation.transpose()

# Create a larger figure with Cartopy
fig, ax = plt.subplots(figsize=(12, 8), subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_extent([-128, -66.5, 24, 50], crs=ccrs.PlateCarree())

# Add map features
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.STATES, linestyle=':')
ax.add_feature(cfeature.LAND, edgecolor='black')

# Plot the total precipitation on the map
c = ax.contourf(lon, lat, total_precipitation, cmap='magma_r', levels=range(1, 11), extend='max', transform=ccrs.PlateCarree())
colorbar = plt.colorbar(c, ax=ax, label='Total Precipitation (inches)', shrink=0.5)  # Adjust the shrink parameter as needed

# Add title and labels
plt.title('Total Precipitation from {} to {}'.format(start_time, end_time))
plt.xlabel('Longitude')
plt.ylabel('Latitude')

# Show the plot
plt.show()

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_physical/ne_50m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_cultural/ne_50m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/50m_cultural/ne_50m_admin_1_states_provinces_lakes.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Zoom in on NE CONUS.

fig, ax = plt.subplots(figsize=(12, 12), subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_extent([-82, -66.5, 35, 48], crs=ccrs.PlateCarree())
ax.set_extent([-76, -68.5, 39, 43], crs=ccrs.PlateCarree())

# Add map features
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.STATES, linestyle=':', linewidth=2)
ax.add_feature(cfeature.LAND, edgecolor='black')
ax.add_feature(USCOUNTIES,edgecolor='grey', linewidth=1)

# Plot the total precipitation on the map   (levels=range(1,10) OR WHAT I HAVE NOW
c = ax.contourf(lon, lat, total_precipitation, cmap='gist_ncar', levels=np.arange(0.1, 100, 1), extend='max', transform=ccrs.PlateCarree())
colorbar = plt.colorbar(c, ax=ax, label='Total Precipitation (inches)', shrink=0.9)  # Adjust the shrink parameter as needed

# Add title and labels
plt.title('Total Precipitation from {} to {}'.format(start_time, end_time))
plt.xlabel('Longitude')
plt.ylabel('Latitude')

# Show the plot
plt.show()

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.11/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_1_states_provinces_lakes.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Downloading file 'us_counties_20m.cpg' from 'https://github.com/Unidata/MetPy/raw/v1.6.2/staticdata/us_counties_20m.cpg' to '/home/runner/.cache/metpy/v1.6.2'.

Downloading file 'us_counties_20m.dbf' from 'https://github.com/Unidata/MetPy/raw/v1.6.2/staticdata/us_counties_20m.dbf' to '/home/runner/.cache/metpy/v1.6.2'.

Downloading file 'us_counties_20m.prj' from 'https://github.com/Unidata/MetPy/raw/v1.6.2/staticdata/us_counties_20m.prj' to '/home/runner/.cache/metpy/v1.6.2'.

Downloading file 'us_counties_20m.shx' from 'https://github.com/Unidata/MetPy/raw/v1.6.2/staticdata/us_counties_20m.shx' to '/home/runner/.cache/metpy/v1.6.2'.

Downloading file 'us_counties_20m.shp' from 'https://github.com/Unidata/MetPy/raw/v1.6.2/staticdata/us_counties_20m.shp' to '/home/runner/.cache/metpy/v1.6.2'.