Visualize One Chunk of a Zarr Variable

In [ ]:
import numpy as np
import cartopy
import cartopy.crs as ccrs
from datetime import datetime,timedelta
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.colors as mcolors
import xarray as xr
import zarr
import s3fs
import boto3
import numcodecs as ncd
import struct
In [ ]:
data = pathToMetadata
x = xr.open_dataset(data)
lat = x.latitude.data
lon = x.longitude.data
In [ ]:
# Input Lat/Lon you are interested in to find chunk

def findChunk(latPoint,lonPoint):
    absLat = np.abs(lat - latPoint)
    absLon = np.abs(lon - lonPoint)
    c = np.maximum(absLon,absLat)
    index = np.where(c == c.min())

    x = int(index[0])
    y = int(index[1])
    
    # Divide by chunk size - ignoring remainder
    chunkX = (x//150)
    chunkY = (y//150)
    
    return chunkX, chunkY


def getLatLons(chunkX,chunkY):
    x0 = chunkX*150
    y0 = chunkY*150
    xRange = x0+150
    yRange = y0+150
    
    chunkLats = lat[x0:xRange,y0:yRange]
    chunkLons = lon[x0:xRange,y0:yRange]
    
    return chunkLats, chunkLons

def nearestGridPoint(ptLat,ptLon,latGrid,lonGrid):
    abslat = np.abs(ptLat - latGrid)
    abslon = np.abs(ptLon - lonGrid)
    c = np.maximum(abslon,abslat)
    ind = np.where(c == c.min())

    i = int(ind[0])
    j = int(ind[1])
    
    return i, j


def perdelta(start, end, delta):
    current = start
    while current < end:
        yield current
        current +=delta
In [62]:
ptLat = 40.2
ptLon = -111.335

x, y = findChunk(ptLat,ptLon)

chunkID = '0.{}.{}'.format(x,y)

newLats, newLons = getLatLons(x,y)

i, j = nearestGridPoint(ptLat,ptLon,newLats,newLons)
In [ ]:
# Set up hourly time stamps
start = datetime(2020, 10, 18, 0, 0)
end = datetime(2020, 10, 20, 15, 0)

dates=[]
for times in perdelta(start, end, timedelta(hours=1)):
    dates.append(datetime.strftime(times, '%m/%d %HZ'))
In [63]:
tempFilePath = yourTempFilePath

KEY = amazonKey
SECRET = amazonSecret
REGION = 'us-west-1'
BUCKET = 'hrrrzarr'
DAY = '20201021'

levelVar = '/surface/PRES/surface/PRES/'


s3 = boto3.resource(
    service_name='s3',
    region_name=REGION,
    aws_access_key_id=KEY,
    aws_secret_access_key=SECRET
    )
    

for hr in range(0,1):
    fileName = DAY+'_{:02d}z_fcst.zarr'.format(hr)
    fileKey = 'sfc/'+DAY+'/'+fileName+levelVar+chunkID
    newFile = tempFilePath+'{:02d}_chunk'.format(hr)
    s3.meta.client.download_file(BUCKET, fileKey, newFile)
In [42]:
chunk = pathToChunk
In [ ]:
# Decompress the chunk you want
data = np.zeros((36,24))

for h in range(0,24):
    chunk = tempFilePath+'{:02d}_chunk'.format(h)
    buf = ncd.blosc.decompress(open(chunk, 'rb').read())

    chunk = np.frombuffer(buf, dtype='<f4')

    if len(chunk) == 810000:
        n = 36
    elif len(chunk) == 405000:
        n = 18
    elif len(chunk) == 765000:
        n = 34
    elif len(chunk) == 382500:
        n = 17

    hrrrVar = np.reshape(chunk,((n,150,150)))
    data[0:n,h] = hrrrVar[:,i,j]

Plot a Timeseries

In [ ]:
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,
                               AutoMinorLocator)

fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(111)

plt.scatter(dates[1:37],data[:,0],s=60,edgecolor='k',zorder=99)
plt.plot(dates[1:37],data[:,0],linestyle = '--',label='00 UTC',zorder=98)
plt.scatter(dates[4:22],data[:18,3],s=60,edgecolor='k',zorder=99)
plt.plot(dates[4:22],data[:18,3],linestyle = '--',label='03 UTC',zorder=98)
plt.scatter(dates[7:43],data[:,6],s=60,edgecolor='k',zorder=99)
plt.plot(dates[7:43],data[:,6],linestyle = '--',label='06 UTC',zorder=98)
plt.scatter(dates[10:28],data[:18,9],s=60,edgecolor='k',zorder=99)
plt.plot(dates[10:28],data[:18,9],linestyle = '--',label='09 UTC',zorder=98)
plt.scatter(dates[13:49],data[:,12],s=60,edgecolor='k',zorder=99)
plt.plot(dates[13:49],data[:,12],linestyle = '--',label='12 UTC',zorder=98)

ax.xaxis.set_major_locator(MultipleLocator(6))
plt.legend()
plt.ylabel('Gust (m/s)',fontsize=14)
plt.xlabel('Time',fontsize=14)

Plot a spatial map

In [ ]:
import cartopy.io.img_tiles as cimgt


fig = plt.figure(figsize=(15,12))
#ax = plt.axes(projection=ccrs.Mercator())
stamen_terrain = cimgt.Stamen('terrain-background')
ax = plt.axes(projection=stamen_terrain.crs)

ax.coastlines()
#ax.set_extent([lon.min()-1,lon.max()+1,lat.min()-1,lat.max()-2])
ax.set_extent([-125,-102,28,50])
states = cartopy.feature.NaturalEarthFeature(
            category='cultural', scale='50m',facecolor='none',
            name='admin_1_states_provinces_shp')

ax.add_image(stamen_terrain,8)
ax.add_feature(states,edgecolor='black',zorder=97,linewidth=1)
ax.add_feature(cartopy.feature.BORDERS,edgecolor='black',zorder=97,linewidth=1)

ax.scatter(iGridLon,iGridLat,s=8,c='black',transform=ccrs.PlateCarree(),zorder=98)
ax.scatter(jGridLon,jGridLat,s=8,c='black',transform=ccrs.PlateCarree(),zorder=98)