How to calculate precipitation rate (mm/hr) from WRF model output variable?

[Ankan]

Dear all,
I want to calculate the precipitation rate (mm/hr) from the WRF model output and validate it with hourly satellite observations (e.g., GPM-IMERG). After searching online, I learned from this conversation (Calculate precipitation from WRF model output) that I need to just add 'RAINC' and 'RAINNC' to get the precipitation (i.e., RAINC+RAINNC). Now, I did this to compute the precipitation. However, after checking the hourly progression of these two variables, , i.e., 'RAINC' and 'RAINNC', I realized that the WRF model is not giving the precipitation value for every hour; rather, it is adding every hour precipitation value. For example, whatever precipitation values are there in the 1st hour, in the 2nd hour, it is also considering the value of the 1st hour as well as the the value for the 2nd hour, and so on. But, I want to extract the exact precipitation value for every hour so that I can see the hourly variation of precipitation rate.
So, can anyone please guide me on how to calculate the precipitation rate (mm/hr) from these WRF model output variables? Any help will be really appreciated. Thank you for your time.
With regards,
Ankan

 

[Ming Chen]

Ankan,
You are right that RAINC and RAINNC are precipitation accumulated from the initiation of WRF run.
You can turn on the option prec_acc_dt = 60 to namelist.input (&physics) and you will have hourly output of prec_acc_c, prec_acc_nc and snow_acc_nc.
Please let me know whether this option works.

 

[kcmonaka]

Hello. I am struggling with handling precipitation. I am just interested displaying precipitation soo i do not know how to handle RAINC and RAINNC. I am interested in just calculating total precipitation for a day. I am running a nested domain (d01 at15km and d02 at 5km) now i have the wrfout files and every calculation i tried using python it shows strange figures. If you have any other code or ncl script please assist. apologies for hijacking the post


python script for adding RAINC and RAINNC
import os
import xarray as xr

# Function to process a single date and domain
def process_date_and_domain(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 variables (RAINNC and RAINC)
rainnc = wrf_data['RAINNC']
rainc = wrf_data['RAINC']

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

# Sum the RAINNC and RAINC to get total precipitation for this time step
total_precip = rainnc + rainc

# 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({
'TOTAL_PRECIP': total_precip,
'XLAT': latitude,
'XLONG': longitude
})
else:
# If data for the day already exists, sum the new total_precip with the existing daily total
day_data[date_part]['TOTAL_PRECIP'] += total_precip

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

# Save the aggregated daily total precipitation 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}_daily_total.nc')
data.to_netcdf(output_file_path)
print(f'Aggregated daily total precipitation data for {date_part} saved to {output_file_path}')

# Main interactive loop
def interactive_process():
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
process_date_and_domain(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()


python script for Visualizing
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:
total_precip = wrf_file.variables['TOTAL_PRECIP'][:] # Read TOTAL_PRECIP 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()

# Replace the file paths with your actual file paths
wrf_file_path = '/home/wrf/Build_WRF/WRF/run/wrfout_d01_2023-12-26_daily_total.nc'
shapefile_path = '/home/wrf/Build_WRF/WRF/run/SHP/Bw_disricts.shp'
wrf_file_date = '2023-12-26'

# 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)

 

[William.Hatheway]

Hello. I am struggling with handling precipitation. I am just interested displaying precipitation soo i do not know how to handle RAINC and RAINNC. I am interested in just calculating total precipitation for a day. I am running a nested domain (d01 at15km and d02 at 5km) now i have the wrfout files and every calculation i tried using python it shows strange figures. If you have any other code or ncl script please assist. apologies for hijacking the post


python script for adding RAINC and RAINNC
import os
import xarray as xr

# Function to process a single date and domain
def process_date_and_domain(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 variables (RAINNC and RAINC)
rainnc = wrf_data['RAINNC']
rainc = wrf_data['RAINC']

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

# Sum the RAINNC and RAINC to get total precipitation for this time step
total_precip = rainnc + rainc

# 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({
'TOTAL_PRECIP': total_precip,
'XLAT': latitude,
'XLONG': longitude
})
else:
# If data for the day already exists, sum the new total_precip with the existing daily total
day_data[date_part]['TOTAL_PRECIP'] += total_precip

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

# Save the aggregated daily total precipitation 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}_daily_total.nc')
data.to_netcdf(output_file_path)
print(f'Aggregated daily total precipitation data for {date_part} saved to {output_file_path}')

# Main interactive loop
def interactive_process():
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
process_date_and_domain(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()


python script for Visualizing
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:
total_precip = wrf_file.variables['TOTAL_PRECIP'][:] # Read TOTAL_PRECIP 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()

# Replace the file paths with your actual file paths
wrf_file_path = '/home/wrf/Build_WRF/WRF/run/wrfout_d01_2023-12-26_daily_total.nc'
shapefile_path = '/home/wrf/Build_WRF/WRF/run/SHP/Bw_disricts.shp'
wrf_file_date = '2023-12-26'

# 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)

can you open this as a new issue as it is different from the original post? The admins like to keep the issues separated to help future users.

Also this might be helpful too:

wrfpython-talk - Google Groups

groups.google.com