""" QPF (Manila) ======================================================== This example demonstrates how to perform operational deterministic QPF up to three hours using raingauge data from Manila and radar data from Subic. """ ########################################################### # Definitions # -------------------------------------------------------- # import os import numpy as np import pandas as pd from pyresample import utils import xarray import cartopy.feature as cfeature import cartopy.crs as ccrs import matplotlib import imageio import matplotlib.pyplot as plt from matplotlib.colors import BoundaryNorm from swirlspy.rad.uf_ph import read_uf_ph from swirlspy.qpe.utils import dbz2rr, rr2rf, locate_file, timestamps_ending from swirlspy.qpe.utils import multiple_acc from swirlspy.obs.rain import Rain from swirlspy.qpf.rover import rover from swirlspy.qpf.sla import sla from swirlspy.core.resample import grid_resample plt.switch_backend('agg') THIS_DIR = os.getcwd() os.chdir(THIS_DIR) start_time = pd.Timestamp.now() ############################################################# # Initialising # --------------------------------------------------- # # This section demonstrates extracting # radar reflectivity data. # # Step 1: Define a basetime. # # Supply basetime basetime = pd.Timestamp('20180811112000').floor('min') ############################################################################## # Step 2: Using basetime, generate timestamps of desired radar files # timestamps_ending() and locate files using locate_file(). # Obtain radar files dir = THIS_DIR + '/../tests/samples/uf_ph/sub/' located_files = [] radar_ts = timestamps_ending( basetime, duration=pd.Timedelta(minutes=60), interval=pd.Timedelta(minutes=10) ) for timestamp in radar_ts: located_files.append(locate_file(dir, timestamp)) ############################################################################# # Step 3: Define the target grid as a pyresample AreaDefinition. area_id = "epsg3123_240km" description = ("A 240 m resolution rectangular grid " "centred at Subic RADAR and extending to 240 km " "in each direction") proj_id = 'epsg3123' projection = ('+proj=tmerc +lat_0=0 ' '+lon_0=121 +k=0.99995 +x_0=500000 ' '+y_0=0 +ellps=clrk66 +towgs84=-127.62,-67.24,' '-47.04,-3.068,4.903,1.578,-1.06 +units=m ' '+no_defs') x_size = 500 y_size = 500 area_extent = (191376.04113, 1399386.68659, 671376.04113, 1879386.68659) area_def = utils.get_area_def( area_id, description, proj_id, projection, x_size, y_size, area_extent ) ############################################################################## # Step 4: Read data from radar files into xarray.DataArray # using read_uf_ph(). # reflectivity_list = [] # stores reflec from read_iris() for filename in located_files: reflec = read_uf_ph( filename, area_def=area_def, coord_label=['easting', 'northing'], indicator='cappi', elevation=2 ) reflectivity_list.append(reflec) ############################################################################ # Step 5: Assigning reflectivity xarrays at the last two timestamps to # variables for use during ROVER QPF. xarray1 = reflectivity_list[-2] xarray2 = reflectivity_list[-1] initialising_time = pd.Timestamp.now() ########################################################### # Running ROVER and Semi-Lagrangian Advection # ------------------------------------------- # # 1. Concatenate two reflectivity xarrays along time dimension. # 2. Run ROVER, with the concatenated xarray as the input. # 3. Perform Semi-Lagrangian Advection using the motion fields from rover. # Combining the two reflectivity DataArrays # the order of the coordinate keys is now ['y', 'x', 'time'] # as opposed to ['time', 'x', 'y'] reflec_concat = xarray.concat([xarray1, xarray2], dim='time') # Rover motion_u, motion_v = rover( reflec_concat ) rover_time = pd.Timestamp.now() # Semi Lagrangian Advection reflectivity = sla( reflec_concat, motion_u, motion_v, steps=30 ) sla_time = pd.Timestamp.now() ############################################################## # Concatenating observed and forecasted reflectivities # ------------------------------------------------------- # # 1. Add forecasted reflectivity to reflectivity_list. # 2. Concatenate observed and forecasted reflectivity # xarray.DataArrays along the time dimension. reflectivity_list.append(reflectivity[1:, ...]) reflectivity = xarray.concat(reflectivity_list, dim='time') concat_time = pd.Timestamp.now() ############################################# # Generating radar reflectivity maps # ----------------------------------- # # #. Define the colour scale and the format of the # reflectivity plots. # #. Initialise the figure and the cartopy GeoAxes. Use the # to_cartopy_crs() method of the AreaDefinition # to obtain the projection system of the axes. # #. Add coastlines using the add_feature() method # of cartopy GeoAxes. # #. Plot data using xarray.plot() and make desired adjustments # to labels. # #. Plot quiver using the quiver() method of cartopy GeoAxes. # # In this example, only hourly images will be plotted. # # Defining colour scale and format levels = [ -32768, 10, 15, 20, 24, 28, 32, 34, 38, 41, 44, 47, 50, 53, 56, 58, 60, 62 ] cmap = matplotlib.colors.ListedColormap([ '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433', '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200', '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) # generate timesteps for plotting timesteps = [basetime+pd.Timedelta(minutes=60*i) for i in range(5)] # Defining quiver parameters qx = motion_u.coords['easting'].values[::5] qy = motion_u.coords['northing'].values[::5] qu = motion_u.values[::5, ::5] qv = motion_v.values[::5, ::5] # Defining coastlines hires = cfeature.GSHHSFeature( scale='h', levels=[1], edgecolor='black', facecolor='none' ) # Plotting images = [] for time in timesteps: # Figure plt.figure(figsize=(28, 21)) plt.axis("equal") # Axes crs = area_def.to_cartopy_crs() ax = plt.axes(projection=crs) # Coastlines ax.add_feature(hires) # Plot reflectivity data quadmesh = reflectivity.sel(time=time).plot( cmap=cmap, norm=norm, extend='neither' ) # Customising labels cbar = quadmesh.colorbar cbar.ax.set_ylabel( reflectivity.attrs['long_name']+'['+reflectivity.attrs['units']+']', fontsize=28 ) cbar.ax.tick_params(labelsize=24) ax.xaxis.set_visible(True) ax.yaxis.set_visible(True) ax.xaxis.set_tick_params(labelsize=24) ax.yaxis.set_tick_params(labelsize=24) ax.xaxis.label.set_size(28) ax.yaxis.label.set_size(28) # Quiver ax.quiver( qx, qy, qu, qv, pivot='mid', regrid_shape=20 ) # Title plt.title( "Quantitative Precipitation Forecast with Reflectivity Fields\n" f"Base time: {basetime}H\n" f"Valid time {pd.Timestamp(time)}H", fontsize=32 ) savepath = os.path.join( THIS_DIR, ("../tests/outputs/rover-output-map-" f"{pd.Timestamp(time).strftime('%Y%m%d%H%M')}.png") ) plt.savefig(savepath) images.append(imageio.imread(savepath)) imageio.mimsave( THIS_DIR + "/../tests/outputs/rover-output-map.gif", images, duration=1 ) radar_image_time = pd.Timestamp.now() ########################################################################## # Accumulating hourly rainfall for 3-hour forecast # ------------------------------------------------ # # Hourly accumulated rainfall is calculated every 30 minutes, the first # endtime is the basetime i.e. T+0min. # # #. Convert reflectivity in dBZ to rainrates in mm/h with dbz2rr(). # #. Changing time coordinates of xarray from start time to endtime. # #. Convert rainrates to rainfalls in 10 mins with rr2rf(). # #. Accumulate hourly rainfall every 30 minutes using multiple_acc(). # Convert reflectivity to rainrates rainrates = dbz2rr(reflectivity, a=300, b=1.4) # Converting the coordinates of xarray from start to endtime rainrates_endtime = rainrates.copy() rainrates_endtime.coords['time'] = \ [ pd.Timestamp(t) + pd.Timedelta(minutes=10) for t in rainrates_endtime.coords['time'].values ] # Convert rainrates to accumulated rainfalls every 10 minutes with rr2rf(). rainfalls = rr2rf(rainrates_endtime, scan_duration=10) # Accumulate hourly rainfall every 30 minutes acc_rf = multiple_acc( rainfalls, basetime, basetime+pd.Timedelta(hours=3) ) acc_time = pd.Timestamp.now() ####################################################################### # Plotting rainfall maps # --------------------------------------- # # 1. Define the colour scheme. # 2. Obtaining the crs of the projection system using to_cartopy_crs() # method of AreaDefinition. # 3. Defining the zoom or area extent. Tuple order is (x0, x1, y0, y1) # as opposed to pyresample (x0, y0, x1, y1). # 4. Initialise figure. # 5. Initialise cartopy GeoAxes. # 6. Adding area extent defined in (3) to axes. # 7. Adding coastlines (this examples uses GSHHS). # 8. Plot xarray using xarray.plot() and make desired adjustments to labels. # 9. Adding title (if required). # Defining the colour scheme levels = [ 0, 0.5, 2, 5, 10, 20, 30, 40, 50, 70, 100, 150, 200, 300, 400, 500, 600, 700 ] cmap = matplotlib.colors.ListedColormap([ '#ffffff', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646', '#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063', '#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd' ]) norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True) # Defining projection crs = area_def.to_cartopy_crs() # Defining zoom extent r = 64000 proj_site = xarray1.proj_site zoom = ( proj_site[0]-r, proj_site[0]+r, proj_site[1]-r, proj_site[1]+r ) images = [] for time in acc_rf.coords['time'].values: # Define figure plt.figure(figsize=(28, 21)) plt.axis('equal') # Plotting axes ax = plt.axes(projection=crs) # setting extent ax.set_extent(zoom, crs=crs) # Adding coastlines ax.add_feature(hires) # Plot data quadmesh = acc_rf.sel(time=time).plot( cmap=cmap, norm=norm, extend='neither' ) # Customising labels cbar = quadmesh.colorbar cbar.ax.set_ylabel( acc_rf.attrs['long_name']+'['+acc_rf.attrs['units']+']', fontsize=28 ) cbar.ax.tick_params(labelsize=24) ax.xaxis.set_visible(True) ax.yaxis.set_visible(True) ax.xaxis.set_tick_params(labelsize=24) ax.yaxis.set_tick_params(labelsize=24) ax.xaxis.label.set_size(28) ax.yaxis.label.set_size(28) # Add title t_minus = pd.Timestamp(time) - \ basetime-pd.Timedelta(minutes=60) t_minus = round(t_minus.total_seconds()) // 60 t_plus = pd.Timestamp(time) - basetime t_plus = round(t_plus.total_seconds()) // 60 plt.title( "Quantitative Precipitation Forecast with Hourly" " Accumulated Rainfall.\n" f"Basetime: {basetime} Valid time: {pd.Timestamp(time)}\n" "Start time: t" f"{t_minus:+} minutes End time: t{t_plus:+} minutes", fontsize=32 ) savepath = THIS_DIR + \ ("/../tests/outputs/forecast-" + f"{pd.Timestamp(time).strftime('%Y%m%d%H%M')}.png") plt.savefig(savepath) images.append(imageio.imread(savepath)) imageio.mimsave( THIS_DIR + "/../tests/outputs/forecast.gif", images, duration=1 ) rf_image_time = pd.Timestamp.now() ############################################################################## # Extract the rainfall values at a specified location # ------------------------------------------------------------------ # # In this example, the rainfall values at the location is assumed # to be the same as the nearest gridpoint. # # 1. Read information regarding the radar stations into a pandas.DataFrame. # 2. Extract the rainfall values at the nearest gridpoint to location # for given timesteps (in this example, 30 minute intervals). # 3. Store rainfall values over time in a pandas.DataFrame. # 4. Plot the time series of rainfall at different stations. # Getting radar station coordinates df = pd.read_csv( os.path.join(THIS_DIR, "../tests/samples/manila_rg_list.csv"), delim_whitespace=True, usecols=[0, 3, 4] ) # Extract rainfall values at gridpoint closest to the # location specified for given timesteps and storing it # in pandas.DataFrame. rf_time = [] for time in acc_rf.coords['time'].values: rf = [] for index, row in df.iterrows(): rf.append(acc_rf.sel( time=time, northing=row[2], easting=row[1], method='nearest' ).values) rf_time.append(rf) rf_time = np.array(rf_time) station_rf = pd.DataFrame( data=rf_time, columns=df.iloc[:, 0], index=pd.Index( acc_rf.coords['time'].values, name='time' ) ) print(station_rf) loc_stn = \ ['BAA', 'BUM', 'PAF', 'QUL', 'ZAP', 'ZAA'] loc_stn_drop = [ stn for stn in station_rf.columns.to_list() if stn not in loc_stn ] df_loc = station_rf.drop(loc_stn_drop, axis=1) print(df_loc) # Plotting time series graph for selected stations ax = df_loc.plot(title="Time Series of Hourly Accumulated Rainfall", grid=True) ax.set_ylabel("Hourly Accumulated Rainfall (mm)") plt.savefig(THIS_DIR+"/../tests/outputs/qpf_time_series.png") extract_time = pd.Timestamp.now() ###################################################################### # Checking run time of each component # -------------------------------------------------------------------- # print(f"Start time: {start_time}") print(f"Initialising time: {initialising_time}") print(f"Rover time: {rover_time}") print(f"SLA time: {sla_time}") print(f"Concatenating time: {concat_time}") print(f"Plotting radar image time: {radar_image_time}") print(f"Accumulating rainfall time: {acc_time}") print(f"Plotting rainfall map time: {rf_image_time}") print(f"Extracting and plotting time series time: {extract_time}") print(f"Time to initialise: {initialising_time-start_time}") print(f"Time to run rover: {rover_time-initialising_time}") print(f"Time to perform SLA: {sla_time-rover_time}") print(f"Time to concatenate xarrays: {concat_time - sla_time}") print(f"Time to plot radar image: {radar_image_time - concat_time}") print(f"Time to accumulate rainfall: {acc_time - radar_image_time}") print(f"Time to plot rainfall maps: {rf_image_time-acc_time}") print(f"Time to extract and plot time series: {extract_time-rf_image_time}")