Creating vertical cross-section along specific lat/lon path

[Michael Pletcher]

Hello all,

I'm working on creating a vertical cross-section along a flight path using lat/lon data collected during the flight but I'm running into an error. 

Here's the code I have so far:

u_wind_ms = getvar(ncfile, 'ua', units = 'm/s')

p_hpa         = getvar(ncfile, 'p',  units = 'hPa')

dawn_data = Dataset("file", 'r')
dawn_lats  = np.array(dawn_data.variables["Profile_Latitude"])
dawn_lons = np.array(dawn_data.variables["Profile_Longitude"])

zonal_cs    = []

for i, j, k, l in zip(dawn_lats, dawn_lons, dawn_lats[1:], dawn_lons[1:]):

    u_zonal_cross = vertcross(u_wind_ms, p_hpa, levels = np.arange(1000, 199.9, -5.), start_point = CoordPair(lat = i, lon = j), end_point = CoordPair(lat = k, lon = l), latlon = True, wrfin = ncfile, meta = True)

    zonal_cs.append(u_zonal_cross)

but I run into the following error:

Traceback (most recent call last):
  File "wrf_crosssec.py", line 99, in <module>
    u_zonal_cross         = vertcross(u_wind_ms, p_hpa, levels = np.arange(1000, 199.9, -5.), start_point = CoordPair(lat = i, lon = j), end_point = CoordPair(lat = k, lon = l), latlon = True, wrfin = ncfile, meta = True)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/metadecorators.py", line 1732, in func_wrapper
    return _set_cross_meta(wrapped, instance, args, kwargs)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/metadecorators.py", line 991, in _set_cross_meta
    levels, autolevels)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/interputils.py", line 236, in get_xy_z_params
    var2dz = _interp2dxy(z, xy)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/decorators.py", line 511, in func_wrapper
    return wrapped(*args, **kwargs)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/decorators.py", line 137, in func_wrapper
    return wrapped(*args, **kwargs)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/decorators.py", line 299, in func_wrapper
    result = wrapped(*new_args, **new_kargs)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/decorators.py", line 392, in func_wrapper
    result = wrapped(*new_args, **new_kargs)
  File "/uufs/chpc.utah.edu/common/home/u1325975/software/pkg/miniconda3/lib/python3.7/site-packages/wrf/extension.py", line 139, in _interp2dxy
    xy)
ValueError: unexpected array size: new_size=4, got array with arr_size=2

Does anyone have an idea what this error means and what could be causing it? What I'm trying to do with the above program is create many vertical cross-sections using lat/lon data collected during a flight and combine all of those cross-sections into one large cross-section.

Any help is greatly appreciated.

Thanks,

Michael

[Michael Pletcher]

Hello all,

I realized that I was using the wrong domain when generating the cross-sections which gave me the above error. However, I'm now running into a different issue.

The intervals between each lat/lon point in the recorded flight data are not uniform, which causes vertcross to output cross-sections that aren't all the same shape. For example, when I print

print(np.array(u_zonal_cross.shape))

in the loop, it outputs these array shapes for example

[161   6]

[161   5]

[161   3]

[161   2]

Which I cannot append to a list, since they have different shapes. Does anyone have an idea on how to get around this? 

Thanks,

Michael

[Michael Pletcher]

I forgot to mention the error I receive when I attempt to print the shape of the appended array

ValueError: could not broadcast input array from shape (161,6) into shape (161,)

[Marco Miani]

Greetings, Michael.

WRF-python package is build on xarray which, in turn, offers splendid features and functionalities when it comes to merge data. Eg.:

import numpy as np
import xarray as xr

da1 = xr.DataArray(np.arange(6).reshape(2, 3), [("latlon", ["lat/lon1", "lat/lon2"]), ("lev", [10, 20, 30])])
da2 = xr.DataArray(np.arange(9).reshape(3, 3), [("latlon", ["lat/lon3", "lat/lon4","lat/lon5"]), ("lev", [10, 20, 30])])
da3 = xr.DataArray(np.arange(12).reshape(4, 3), [("latlon", ["lat/lon6", "lat/lon7","lat/lon8","lat/lon9"]), ("lev", [10, 20, 30])])
print(da1)
print("\n\n")
print(da2)
print("\n\n")
print(da3)

#merge them all along the same direction
total = xr.concat([da1, da2, da3], dim='latlon') # the vertical coordinate is always consistent, so use that to "chain" and concat along varying dimension...

#print
print(total)

print("\nNew dimensions (aggregated) are: ",total.shape)

Hope this helped.
Bye

Marco

[Michael Pletcher]

Hi Marco,

Thanks so much for the input.

I've been trying to understand the data structure you use in your variables. Here's what I think each part means in the context of my data:

1. The "np.arange(6).reshape(2, 3)" is the reshaped vertcross output data. In my case, if the output data's shape was (161,2) then I would do np.arange(322).reshape(u_zonal_cross)

2. The ("latlon", ["lat/lon1", "lat/lon2"]) are the two lat/lon points generated from vertcross. In my case, this would be u_zonal_cross[0,:] with a length of 2 in this example.

3, The "lev", [10, 20, 30] are the vertical coordinates generated from vertcross. In my case, this would be u_zonal_cross.shape[:,0] with a length of 161.

Can you correct me if I'm wrong? When I run the following code

u_zonal_cross = vertcross(u_wind_ms, p_hpa, levels = np.arange(1000, 199.9, -5.), start_point = CoordPair(lat = i, lon = j), end_point = CoordPair(lat = k, lon = l), latlon = True, wrfin = ncfile, meta = True)

cross_5            = xr.DataArray(np.arange(805).reshape(u_zonal_cross.T), [("latlon", u_zonal_cross[0,:]), ("lev", u_zonal_cross[:,0])])

I get the warning UserWarning: 'latlon' is set to True, but 'field3d' is not of type xarray.DataArray and contains no coordinate information warnings.warn("'latlon' is set to True, but 'field3d' is "

and the program stops with the error

Traceback (most recent call last):

  File "wrf_crosssec.py", line 136, in <module>
    cross_5  = xr.DataArray(np.arange(805).reshape(u_zonal_cross.T), [("latlon", u_zonal_cross[0,:]), ("lev", u_zonal_cross[:,0])])
TypeError: 'DataArray' object cannot be interpreted as an integer

Let me know if you need anything else from me.

Thank you,

Michael

[Marco Miani]

Hi Michael.

I see my last message only served generating more confusion, instead of actually helping you in your problem solution. Let me clarify.

1. The code I had originally posted, was a mere example to demonstrate one of the (many) possible ways of solving this, using xarray. 
2. You could have used that as inspiration in shaping a possible solution to your problem. However, I may deduce you are not confident with xarray, so I decided to code it for you (a very easy and simple version). I have attached 2 figures to this message: please refer to them, when reading the following text.
   
3. Please be my guest and follow me in this little trip across the Italian boot where, on board of our drone, we are overflying 3 different regions (NORTH, CENTER and SOUTH)
4. These 3 flight paths (measuring events) took place at 3 different times (e.g., Mon - Wed - Friday), so the time for each flight set is monotonically increasing, but we may assume that flight duration was smaller then one single output time-step of WRF model. This assumption is important.
   
5. These 3 flight paths were different in duration, as the battery (or fuel on board) was different for each flight. This means that, at constant cruising speed, the corresponding vertical cross section has stretches over different lengths (i.e. same number of vertical levels, but larger amount of horizontal coordinates)
6. When extracting from WRF, we have to set a different time (see point 4).
7. When extracting using vertcross function a new coordinate variable is assigned to the returned vertical cross section: dim = 'cross_line_idx'. We will be using that new dimension to concatenate (i.e. stitch all datasets together), and join them.
   
8. The 3 arrays you had posted in your first message ("u_zonal_cross") already are xarray type of variables. This means that, with the right syntax, you could persuade python to concat them. My advice for you is not "downgrading" to numpy (np(u_zonal_cross) etc.... ) but, instead, remaining in the xarray realm. This would make everything easier when manipulating data. xarray is a very powerful tool and was specifically designed and tailored to handle multidimensional, labelled, structured data like these. Also, xarray allows you to include metadata, which are crucial when manipulating large data and combining
9. I now take the 3 flight paths and concatenate them together, after stating the dimension along which I wan to perform this operation: 
   
   total = xr.concat([cross1, cross2,cross3], dim = 'cross_line_idx')     - you see how easy it is, within xarray, to combine data.
   
10. You will see that the extracted cross sections has lat lons (xaxis) of ALL the points shown on the map.
    

The result is the cross section I have enclosed to this email. There is also a map. Please note that mapping and plotting data is a functionality that is natively embedded into xarray (using either basemap, or cartopy - I strongly suggest the second option). You just manipulate your variable (concat, slice, mean etc), and then with a one-liner, you plot it in the projection you most like. No nasty coding or coordinate/projection transformations. 

I hope this helps you. If you need more info or have doubts, please feel free to respond to this email. Figures follow.

Happy modelling

Bye bye

Marco