metpy

how to make a cross section of frontogenesis by metpy


I want to make a cross section of frontogenesis my code is

  domain_1 = os.path.abspath(filenames_in)
info = os.path.join(domain_1, filename)
data = xr.open_dataset(info)
data = data.metpy.parse_cf().squeeze()
print(data)
p1 = data['level'].values[:]* units.hPa
start = (56, -9)
end = (56.5, -6)
cross = cross_section(data, start, end).set_coords(('latitude', 'longitude'))
r,le,t,w =xr.broadcast(cross['r'],cross['level'],cross['t'],cross['w'])
# print('r',r.shape)
# print(le)
# print('w',w.shape)
# print(w)
theta=mpcalc.potential_temperature(level,t)
cross['Potential_temperature'] = xr.DataArray(theta, coords=t.coords, dims=t.dims)
ptemp=cross['Potential_temperature']
cross['u'] = cross['u']
u = cross['u']
cross['v'] = cross['v']
v = cross['v']
# Set subset slice for the geographic extent of data to limit download
lon_slice = slice(-40, 20)
lat_slice = slice(70, 40)
# Grab lat/lon values (GFS will be 1D)
lats = data.latitude.sel(latitude=lat_slice).values
lons = data.longitude.sel(longitude=lon_slice).values
# Compute dx and dy spacing for use in vorticity calculation
dx, dy = mpcalc.lat_lon_grid_deltas(lons,lats)
front = mpcalc.frontogenesis(ptemp, u, v, dx[None,:, :], dy[None,:, :], x_dim=-1, y_dim=-2)
tmpk = data.t.metpy.sel(
    latitude=lat_slice, longitude=lon_slice).metpy.unit_array.squeeze()
uwnd = data['u'].metpy.sel(
    latitude=lat_slice, longitude=lon_slice).metpy.unit_array.squeeze()
vwnd = data['v'].metpy.sel(
    latitude=lat_slice, longitude=lon_slice).metpy.unit_array.squeeze()
potent_temp=mpcalc.potential_temperature(p1[:,None,None], tmpk)
fronto = mpcalc.frontogenesis(potent_temp, uwnd, vwnd, dx[None, :, :], dy[None, :, :], x_dim=-1, y_dim=-2)
# f=cross_section(fronto, start, end).set_coords(('latitude', 'longitude'))

there are always some bugs. how to fix it I do not know if i need to calculate the frontogenesis first and then use cross-function or use cross function first and then calculate frontogenesis. I tries both of them but I cannot get the cross section.

ValueError: operands could not be broadcast together with shapes (0,120,241) (19,100) 

if i calculate the frontogenesis first and then use cross-function

lon_slice = slice(-40, 20) 
lat_slice = slice(70, 40) 
subset = data.metpy.sel(latitude=lat_slice, longitude=lon_slice)
subset['potential_temperature'] = mpcalc.potential_temperature(
subset['level'],
subset['t'])
subset['frontogenesis'] = mpcalc.frontogensis(
subset['potential_temperature'],
subset['u'],
subset['v'])
start = (56, -9) end = (56.5, -6)
cross = cross_section(subset, start, end).set_coords(('latitude', 'longitude'))

i always meet this error

ValueError: Data missing required coordinate information. Verify that your data have been parsed by MetPy with proper x and y dimension coordinates and added crs coordinate of the correct projection for each variable.

Solution

  • Given that MetPy's frontogensis calculation operates on a 2D horizontal grid, you will indeed need to calculate frontogensis prior to taking the cross section. Similarly, subsetting should be done prior to calculations, as MetPy does not yet support delayed computations. Additionally, instead of manually handling grid deltas and unit arrays, you will typically run into less errors when you let MetPy handle those for you, like so:

    domain_1 = os.path.abspath(filenames_in)
    info = os.path.join(domain_1, filename)
    data = xr.open_dataset(info)
    data = data.metpy.parse_cf().squeeze()
    print(data)
    
    lon_slice = slice(-40, 20)
    lat_slice = slice(70, 40)
    subset = data.metpy.sel(latitude=lat_slice, longitude=lon_slice)
    
    subset['potential_temperature'] = mpcalc.potential_temperature(
        subset['level'],
        subset['t']
    )
    subset['frontogenesis'] = mpcalc.frontogensis(
        subset['potential_temperature'],
        subset['u'],
        subset['v']
    )
    
    start = (56, -9)
    end = (56.5, -6)
    cross = cross_section(subset, start, end).set_coords(('latitude', 'longitude'))
    

    (This assumes that your data have proper 'units' attributes and aligned grids in horizontal and vertical. If this is not the case, some extra data cleaning will be needed after opening.)