Precipitation

[KC181818]

Please please assist with scripts to process precipitation. I have installed NCL and python but am still strugling. my main issue is that when i generate scripts to extract RAINC (not cummulative time step one "RAINNC") the resulting map is that of accumulated rain. but i just want to get RAINC for a specific day. see code below

please see code to extract RAINC

import os
import xarray as xr

Function to process a single date and domain for RAINC only

def process_date_and_domain_rainc(base_path, date_part, domain):
# Search for files matching the selected date and domain
wrf_files = [file_name for file_name in os.listdir(base_path)
if file_name.startswith(f'wrfout_{domain}_{date_part}')]

if not wrf_files:
    print(f'No files found for the date: {date_part} and domain: {domain}')
    return

# Dictionary to store cumulative data for the selected day
day_data = {}

for file_name in wrf_files:
    # Path to the input WRF NetCDF file
    input_file_path = os.path.join(base_path, file_name)

    # Open the WRF NetCDF file using xarray
    wrf_data = xr.open_dataset(input_file_path)

    try:
        # Extract the precipitation variable (RAINC)
        rainc = wrf_data['RAINC']

        # Extract latitude (XLAT) and longitude (XLONG) variables
        latitude = wrf_data['XLAT']
        longitude = wrf_data['XLONG']

        # If this is the first file for the day, initialize the data for that day
        if date_part not in day_data:
            day_data[date_part] = xr.Dataset({
                'RAINC_TOTAL': rainc,
                'XLAT': latitude,
                'XLONG': longitude
            })
        else:
            # If data for the day already exists, sum the new rainc with the existing daily total
            day_data[date_part]['RAINC_TOTAL'] += rainc

    except KeyError:
        print(f'Error: RAINC, XLAT, or XLONG not found in file: {file_name}')
    finally:
        # Close the opened NetCDF file
        wrf_data.close()

# Save the aggregated daily total RAINC data for the day to a NetCDF file
for date_part, data in day_data.items():
    output_file_path = os.path.join(base_path, f'wrfout_{domain}_{date_part}_RAINC_daily_total.nc')
    data.to_netcdf(output_file_path)
    print(f'Aggregated daily total RAINC data for {date_part} saved to {output_file_path}')

Main interactive loop

def interactive_process_rainc():
base_path = '/data/run/' # Directory path for input and output files

while True:
    # Prompt for the date
    date_part = input("Enter the date (YYYY-MM-DD) you want to process, or 'quit' to exit: ")
    if date_part.lower() == 'quit':
        print("Exiting the script.")
        break

    # Prompt for the domain (e.g., d01 or d02)
    domain = input("Enter the domain (e.g., d01 or d02): ").strip()

    # Process the selected date and domain for RAINC
    process_date_and_domain_rainc(base_path, date_part, domain)

    # Ask if the user wants to process another date
    continue_choice = input("Do you want to process another date? (yes/no): ").strip().lower()
    if continue_choice != 'yes':
        print("Exiting the script.")
        break

Run the interactive process

interactive_process_rainc()

please see code to visualize

import netCDF4 as nc
import matplotlib.pyplot as plt
import geopandas as gpd
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from datetime import datetime
from matplotlib.colors import ListedColormap
import matplotlib.ticker as mticker

Function to read the total precipitation data from the aggregated WRF output file

def read_total_precipitation(wrf_file_path):
with nc.Dataset(wrf_file_path, 'r') as wrf_file:
# Adjusted variable name from 'TOTAL_PRECIP' to 'RAINC_TOTAL'
total_precip = wrf_file.variables['RAINC_TOTAL'][:] # Read RAINC_TOTAL directly
lat = wrf_file.variables['XLAT'][:]
lon = wrf_file.variables['XLONG'][:]

return total_precip, lat, lon

Function to plot precipitation on a map using cartopy

def plot_precipitation_on_map(precipitation_data, lat, lon, shapefile_path, wrf_file_date, legend_colors,
title='WRF Total Precipitation_D01_15km',
colorbar_label='Total Precipitation (mm)'):
# Debugging: Print shapes of the arrays
print("Precipitation Data Shape:", precipitation_data.shape)
print("Latitude Shape:", lat.shape)
print("Longitude Shape:", lon.shape)

# Ensure matching dimensions (remove the time dimension if needed)
precipitation_data = precipitation_data.squeeze()
lat = lat.squeeze()
lon = lon.squeeze()

# Ensure dimensions match after squeezing
if lat.shape != precipitation_data.shape or lon.shape != precipitation_data.shape:
    raise ValueError("Dimension mismatch between latitude, longitude, and precipitation data.")

# Update the colormap and ensure vmax is set to 300
cmap = ListedColormap(legend_colors)

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='-')
ax.add_feature(cfeature.LAND, edgecolor='black')

# Plot the precipitation data with vmax set to 300
im = ax.pcolormesh(lon, lat, precipitation_data,
                   cmap=cmap, transform=ccrs.PlateCarree(), vmax=300)

# Add the shapefile overlay
shapefile = gpd.read_file(shapefile_path)
shapefile.plot(ax=ax, facecolor='none', edgecolor='black', linewidth=0.5)

# Set extent for visualization (you can customize this based on your region)
#ax.set_extent([10, 40, -30, -16], crs=ccrs.PlateCarree())


# Set extent for visualization (you can customize this based on your region)
ax.set_extent([5, 60, -45, -5], crs=ccrs.PlateCarree())

# Add gridlines with labels
gl = ax.gridlines(draw_labels=True, linewidth=0)
gl.top_labels = False
gl.right_labels = False
gl.left_labels = True
gl.bottom_labels = True

# Set gridline intervals to show whole numbers for latitude and longitude
gl.ylocator = mticker.MaxNLocator(integer=True)
gl.xlocator = mticker.MaxNLocator(integer=True)

# Add colorbar with the updated range and label
cbar = plt.colorbar(im, orientation='horizontal')
cbar.set_label(colorbar_label)

# Add title and labels
date_str = datetime.strptime(wrf_file_date, '%Y-%m-%d').strftime('%Y-%m-%d')
plt.title(f'{title} ({date_str})', fontsize=16, loc='center')
plt.xlabel('Longitude', fontsize=12)
plt.ylabel('Latitude', fontsize=12)

# Display the plot
plt.show()

Adjusted file path to match the output naming convention from the previous steps

wrf_file_path = '/home/wrf/Build_WRF/WRF/run/wrfout_d02_2023-12-27_RAINC_daily_total.nc'
shapefile_path = '/home/wrf/Build_WRF/WRF/run/SHP/Bw_disricts.shp'
wrf_file_date = '2024-03-15'

Read total precipitation data from the NetCDF file

total_precip_data, lat, lon = read_total_precipitation(wrf_file_path)

Define a color scheme for precipitation

precipitation_colors = [(1.0, 1.0, 1.0), (0.6, 0.6, 1.0), (0.0, 0.0, 1.0),
(0.0, 0.8, 0.0), (0.8, 1.0, 0.8), (1.0, 1.0, 0.0),
(1.0, 0.8, 0.0), (1.0, 0.4, 0.0), (1.0, 0.0, 0.0)]

Plot precipitation on map with shapefile overlay

plot_precipitation_on_map(total_precip_data, lat, lon, shapefile_path, wrf_file_date, precipitation_colors)