Questions about wrf.vertcross and having values below sea-level.

[Bas Cam]

Hi all,

Currently I am working on a research project about a specific wind near the coast of the Caspian Sea, Kazakhstan. In order to study this wind I also make vertical cross sections of an area (wrf.vertcross). Everything works fine, except that the plotted wind speed (or any other 3d variable) only plots a selected range in the atmosphere until sea level. But.. the Caspian Sea is -28 below sea level, so the vertical cross section does not show the lowest 28m of the figure. Does anyone know how to implement this part or extrapolate it below sea level? Or does it have to do that the WRF-output file does not go lower than that and thus the initial WRF settings have to be adjusted in order to get the values below sea level?  

Below the code can be found. In addition I also added the figure in the attachments. Many thanks in advance!

from netCDF4 import Dataset

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.cm import get_cmap

import cartopy.crs as ccrs

import cartopy

from cartopy.feature import NaturalEarthFeature

from wrf import (getvar, to_np, vertcross, smooth2d, CoordPair, GeoBounds,

                 get_cartopy, latlon_coords, interpline, cartopy_xlim, cartopy_ylim)

# Open the NetCDF file

WRF_File1 = Dataset("wrfout_bautino.kz.default_370m_20210127_2100.nc")

# Set the start point and end point for the cross section

start_point = CoordPair(lat=44.627243, lon=50.217994)

end_point = CoordPair(lat=44.55243, lon=50.460756)

# Get the WRF variables

ht = getvar(WRF_File1, "z", timeidx=-1)

ter = getvar(WRF_File1, "ter", timeidx=-1)

tc = getvar(WRF_File1, "tc", timeidx=-1)

wspd =  getvar(WRF_File1, "wspd_wdir", units="m s-1",timeidx=-1)[0,:]

pres = getvar(WRF_File1, "pressure", timeidx=-1)

# Compute the vertical cross-section interpolation.  Also, include the

# lat/lon points along the cross-section in the metadata by setting latlon

# to True.

wspd_cross = vertcross(wspd, ht, wrfin=WRF_File1, start_point=start_point,

                       end_point=end_point, latlon=True, meta=True)

# To remove the slight gap between the wspd contours and terrain due to the

# contouring of gridded data, a new vertical grid spacing, and model grid

# staggering, fill in the lower grid cells with the first non-missing value

# for each column.

# Make a copy of the z cross data. Let's use regular numpy arrays for this.

wspd_cross_filled = np.ma.copy(wspd_cross)

# For each cross section column, find the first index with non-missing

# values and copy these to the missing elements below.

for i in range(wspd_cross_filled.shape[-1]):

    column_vals = wspd_cross_filled[:,i]

    # Let's find the lowest index that isn't filled. The nonzero function

    # finds all unmasked values greater than 0. Since 0 is a valid value

    # for wspd, let's change that threshold to be -200 wspd instead.

    first_idx = int(np.transpose((column_vals > -200).nonzero())[0])

    wspd_cross_filled[0:first_idx, i] = wspd_cross_filled[first_idx, i]

# Get the terrain heights along the cross section line

ter_line = interpline(ter, wrfin=WRF_File1, start_point=start_point,

                      end_point= end_point)

# Get the lat/lon points

lats, lons = latlon_coords(tc)

# Get the cartopy projection object

cart_proj = get_cartopy(tc)

# Create a figure that will have 3 subplots

fig = plt.figure(figsize=(12,9))

ax_wspd = fig.add_subplot(1,1,1)

#ax_dbz = fig.add_subplot(1,1,1)

# Make the contour plot for wind speed

levels = np.linspace(6, 31, num=26)

#levels = np.linspace(-10, 12, num=23)

xs = np.arange(0, wspd_cross.shape[-1], 1)

ys = to_np(wspd_cross.coords["vertical"])

wspd_contours = ax_wspd.contourf(xs,

                                 ys,

                                 to_np(wspd_cross_filled),20,levels=levels,cmap=get_cmap("jet"))

# Add the color bar

cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd)

cb_wspd.ax.tick_params(labelsize=10)

cb_wspd.set_label('m s-1', rotation=0)

# Fill in the mountain area

ht_fill = ax_wspd.fill_between(xs, -30, ter_line,

                                facecolor="saddlebrown")

# Set the x-ticks to use latitude and longitude labels

coord_pairs = to_np(wspd_cross.coords["xy_loc"])

x_ticks = np.arange(coord_pairs.shape[0])

x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]

# Set the desired number of x ticks below

num_ticks = 5

thin = int((len(x_ticks) / num_ticks) + .5)

coord_pairs = wspd_cross.coords["xy_loc"]

x_ticks = np.arange(coord_pairs.shape[0])

x_labels = [pair.latlon_str() for pair in to_np(coord_pairs)]

ax_wspd.set_xticks(x_ticks[0:150:15])

ax_wspd.set_xticklabels(x_labels[0:150:15], rotation=22.5, fontsize=9)

# Set the x-axis and  y-axis labels

ax_wspd.set_xlabel("Latitude, Longitude", fontsize=12)

ax_wspd.set_ylabel("Height (m)", fontsize=12)

ax_wspd.set_ylim([-30,700])

# Add a title

ax_wspd.set_title("Cross-Section of wind speed over Bautino 27th of January 2100", {"fontsize" : 14})

[Thomas Owen Mazzetti]

define the levels of your cross section in the vertcross call by defining levels as an array including values of your vertical coordinate "height/Z" below zero.

--
You received this message because you are subscribed to the Google Groups "wrfpython-talk" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wrfpython-tal...@ucar.edu.
To view this discussion on the web visit https://groups.google.com/a/ucar.edu/d/msgid/wrfpython-talk/17ee2bd2-fe7e-4c27-b2a6-2769147dc433n%40ucar.edu.

[Bas Cam]

Thank you Thomas, that works! Luckily a simple solution..

Op donderdag 18 maart 2021 om 16:35:47 UTC+1 schreef thomasowe...@gmail.com:

[Thomas Owen Mazzetti]

you are welcome. 

Just a tip, before going to ask someone a question,

 look at the documentation for the functions you are using to see what options are available and to see if you are using it the way it was intended. 

Most of the times the answers can be pretty obvious, and only when things get way complex like using numerous functions, data formats,...

might an issue be pretty difficult to solve for someone who is unfamiliar with a particular library or something.