Wind vector rotation troubles (rotate_vector)

[JacobCarley-NOAA]

I am currently working with a gridded data set that has its wind vectors already in an earth-relative frame. When looking to rotate the wind vectors appropriately for plotting on a Basemap projection I've run into a strange problem where the vectors, quite simply, look rather odd. I ran a few tests with simple test points to make sure the rotate_vector routine appeared to be working correctly, and it seemed fine. But when I run a whole model output grid through the rotate_vector routine, it seems to produce wind fields which don't look quite right (see below).

Here are some stats associated with the point that has the red star on the plot with the title 'Rotated vectors'

Lat and Lon: 39.11, -70.0144

Original U and V:     1.30,  7.14
Rotated U and V:    -1.33 , 7.13

Now note how the rotated wind at this point is aligned along the -70 meridian, which would suggest that in earth relative terms we might expect the wind components to be more along the lines of u=0 and v=7. However, judging from above that's certainly not the case.

Now, if I rotate that same vector, by itself, I get the following rotated wind:

Rotated U and V: -0.03 7.26

And see the second plot for a visual:

It would seem desirable to have the same behavior for the single point that we have for the gridded set of winds.

The only thing special about the gridded data, that I have noticed, is that the latitudes can vary in a non-standard way with increasing x dimension, west to east (e.g. 35.5N, 34.2N, 33.1N, 36.2N, etc.). It wasn't clear if this was okay within rotate_vector. In some simple tests it didn't seem to be a problem.

It would be nice to be able to use this routine for not only a set of gridded data, but observations as well. It's worth noting that observations wouldn't necessarily be ordered in a nice, regular way given the nature that observations tend to be irregularly spaced (e.g. surface stations co-located with airports).

Unfortunately I don't have any suggestions for a solution, but am hopeful others may have an idea.

Below is the snippet of code used to generate the first plot above.

Thanks! Jacob

import nemsio
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# Get the input file
f='restart_file'
nio=nemsio.nemsfile(f)

skip=25
u10=u[::skip,::skip,0]
v10=v[::skip,::skip,0]
lats=nio.lats[::skip,::skip]
lons=nio.lons[::skip,::skip]

# Create the figure
fig=plt.figure(figsize=(18, 6)) 

# Domain covers NE CONUS
llcrnrlon=-84.0
llcrnrlat=35.0
urcrnrlon=-60.0
urcrnrlat=49.0
res='l'

m = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,\
          rsphere=(6378137.00,6356752.3142),\
          resolution=res,projection='lcc',\
          lat_1=25.0,lon_0=-95.0)

#The vector rotation to the Basemap projection just specified
u10_rot, v10_rot, x, y = m.rotate_vector(u10, v10, lons, lats, returnxy=True)

parallels = np.arange(-80.,90,5.)
meridians = np.arange(0.,360.,5.)

# - First sublot is without rotation
ax = fig.add_subplot(121)
ax.set_title('Without rotation')

m.drawmapboundary(fill_color='aqua')
m.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
m.drawcoastlines(color = '0.15')

m.drawparallels(parallels)
m.drawmeridians(meridians)

m.barbs(x, y, u10, v10, pivot='middle', barbcolor='black',zorder=10)

# - Second subplot is with rotation

ax = fig.add_subplot(122)
ax.set_title('Rotated vectors')

m.drawmapboundary(fill_color='aqua')
m.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
m.drawcoastlines(color = '0.15')
m.drawparallels(parallels)
m.drawmeridians(meridians)
m.barbs(x, y, u10_rot, v10_rot, 
    pivot='middle', barbcolor='black',zorder=10)

m.scatter(-70.0144,39.11,s=175,color='red',marker='*',latlon=True)

[micahcochran]

How would someone go about installing the nemsio package?

[JacobCarley-NOAA]

Hi Micah,

It's a local package put together for reading files associated with a model I work with. Unfortunately, it's not in a form which is community friendly (and the file associated with it is large, ~28G).

However, I can provide a numpy.savez_compressed file that contains only the necessary fields along with the associated python to open and read it. Would this be suitable? The file comes out to about 36M in size. The plotted winds are from a different date, but issue is still present.

-Jacob

[guziy]

Yes @jrcarley, it would be great )) I think there is lots of people who would like to look at it.

Cheers

[JacobCarley-NOAA]

Alrighty - thanks for looking into this!

GitHub seems to prefer file sizes <10MB so I've tossed it onto DropBox. To get the file, use the following link:

https://www.dropbox.com/s/nyr5gta8bqoyash/rotate_vector_inputs.npz?dl=0

I've updated the code I pasted in my previous post below so that the data at the link above may be read and plotted:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# Read the npz file
dat=np.load('rotate_vector_inputs.npz')
u=dat['u']
v=dat['v']
nio_lons=dat['lons']
nio_lats=dat['lats']

# Skip factor so we do not overwhelm the plot (this is ~ a 3 km grid)
skip=25
u10=u[::skip,::skip,0]
v10=v[::skip,::skip,0]
lats=nio_lats[::skip,::skip]
lons=nio_lons[::skip,::skip]

# Create the figure
fig=plt.figure(figsize=(18, 6)) 

# Domain covers NE CONUS
llcrnrlon=-84.0
llcrnrlat=35.0
urcrnrlon=-60.0
urcrnrlat=49.0
res='l'

m = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,\
          rsphere=(6378137.00,6356752.3142),\
          resolution=res,projection='lcc',\
          lat_1=25.0,lon_0=-95.0)

#The vector rotation to the Basemap projection just specified
u10_rot, v10_rot, x, y = m.rotate_vector(u10, v10, lons, lats, returnxy=True)

parallels = np.arange(-80.,90,5.)
meridians = np.arange(0.,360.,5.)

# - First sublot is without rotation
ax = fig.add_subplot(121)
ax.set_title('Without rotation')

m.drawmapboundary(fill_color='aqua')
m.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
m.drawcoastlines(color = '0.15')

m.drawparallels(parallels)
m.drawmeridians(meridians)

m.barbs(x, y, u10, v10, 
    pivot='middle', barbcolor='black',zorder=10)

# - Second subplot is with rotation

ax = fig.add_subplot(122)
ax.set_title('Rotated vectors')

m.drawmapboundary(fill_color='aqua')
m.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
m.drawcoastlines(color = '0.15')
m.drawparallels(parallels)
m.drawmeridians(meridians)
m.barbs(x, y, u10_rot, v10_rot, 
    pivot='middle', barbcolor='black',zorder=10)

m.scatter(-70.0144,39.11,s=175,color='red',marker='*',latlon=True)
plt.show()

[efiring]

Have you tried appending .copy() to each array after subsampling, so that you will have C-contiguous arrays? There is a bug that is fixed in master but not in 1.07 involving arrays that are not in C order.

[guziy]

@jrcarley: It seems like the plots with the data from dropbox are different? Or is it because my version of basemap is different (I am using the latest master)? Could you please, also specify the indices corresponding to the star (i, j)?

https://github.com/guziy/PyNotebooks/blob/master/basemap_issue_%23269.ipynb

Cheers

[guziy]

I've figured out the indices, and the plot for 1 point looks very similar to the field plots...

https://github.com/guziy/PyNotebooks/blob/master/basemap_issue_%23269.ipynb

[JacobCarley-NOAA]

@guziy, Sorry I had meant to mention that the date corresponding to the data I had provided is different from the date in the first example. What is important, though, is that the problem remains. My apologies for not specifying the i,j location. On the thinned/subsampled array the points are: i=64, j=40. For the example I attached, at the point of the star we have:

I and J of the thinned/subsampled grid:  (i=64, j=40)
Lat and Lon: 39.11, -70.0144

Original U and V:     7.39,  4.12 
Rotated U and V:     6.17,  5.80

Though the problematic barbs are eastward of the red star for the example data provided here.

@efiring, I did attempt .copy() after your suggestion and unfortunately did not see a change in the behavior. I have also tried passing the full set of winds through the rotation prior to subsampling and the issue remained.

To update: I've a routine that is separate from Basemap which will do the rotation correctly - but it's not general and depends upon the projection being Lambert Conic Conformal or Polar Stereographic. Below is an example with comparisons of un-rotated (leftmost), Basemap rotated (center), and the alternate rotation method (far right):

Note the differences in the barbs east of the red star. The Basemap rotation indicates the winds are mostly southerly while the the alternate approach suggests more a southwesterly wind. The simplified, working code for the alternate rotation is as follows:

def alt_rotate(true_lat,lov_lon,earth_lons,uin,vin,proj,inverse=False):
  #  Rotate winds from LCC relative to earth relative (or vice-versa if inverse==true)
  #
  # Input args:
  #  true_lat = True latitidue for LCC projection (single value in degrees)
  #  lov_lon  = The LOV value (e.g. - -95.0) (single value in degrees)
  #                              "Lov = orientation of the grid; i.e. the east longitude value of
  #                              the meridian which is parallel to the Y-axis (or columns of the grid)
  #                              along which latitude increases as the Y-coordinate increases (the 
  #                              orientation longitude may or may not appear on a particular grid).
  #             
  #  earth_lons = Earth relative longitudes (in degrees)
  #  uin, vin     = Input winds to rotate
  #  
  # Returns:
  #  uout, vout = Output, rotated winds
  #-----------------------------------------------------------------------------------------------------

  if lov_lon > 0.: lov_lon=lov_lon-360.  
  dtr=np.pi/180.0             # Degrees to radians

  # Compute rotation constant which is also
  # known as the Lambert cone constant.  In the case 
  # of a polar stereographic projection, this is one.

  if proj.lower()=='lcc':
    rotcon_p=np.sin(true_lat*dtr)
  elif proj.lower() in ['stere','spstere', 'npstere']:
    rotcon_p=1.0

  angles = rotcon_p*(earth_lons-lov_lon)*dtr 
  sinx2 = np.sin(angles)
  cosx2 = np.cos(angles)

  # Steps below are element-wise products
  if inverse==False:
    # Return the earth relative winds
    uout = cosx2*uin+sinx2*vin
    vout =-sinx2*uin+cosx2*vin
  elif inverse==True:
    # Return the grid relative winds
    uout = cosx2*uin-sinx2*vin
    vout = sinx2*uin+cosx2*vin

  return uout,vout

The alternate method/routine is just a point of reference and not a suggested modification to Basemap, especially since it's not terribly general. To use the routine in the example I've provided, simply call as follows (and update the other fig.add_subplot() settings from 12x to 13x:

  u10_rot_alt, v10_rot_alt = alt_rotate(lat_1,lon_0,lons,u10,v10,'lcc',inverse=True)
  ax = fig.add_subplot(133)
  ax.set_title('Alternate Rotated vectors')

  m.drawmapboundary(fill_color='aqua')
  m.fillcontinents(color='#cc9955', lake_color='aqua', zorder = 0)
  m.drawcoastlines(color = '0.15')
  m.drawparallels(parallels)
  m.drawmeridians(meridians)

  m.barbs(x, y, u10_rot_alt, v10_rot_alt, 
      pivot='middle', barbcolor='black',zorder=10)
  m.scatter(-70.0144,39.11,s=175,color='red',marker='*',latlon=True)
  plt.show()

[Chun-ChihWang]

Hi all,

I don’t know if I can reopen this issue.

I was recently doing some tests with Basemap’s rotate_vector utility on my WRF simulation data. I plan to plot the 10-m wind barbs on a map that has the same projection as the model grid. In my understanding, if someone is plotting wind vectors/barbs from a WRF on a map with the same projection as the model itself, then no rotation needs to be done to the model u and v wind components since they are already grid relative (Please correct me if I am wrong on this).

Out of curiosity, I converted my model u and v from grid-relative to Earth-relative using the following formula according to this website: https://www-k12.atmos.washington.edu/~ovens/wrfwinds.html Uearth = Ucosalpha - Vsinalpha Vearth = Vcosalpha + Usinalpha , where cosalpha and sinalpha are given in the WRF output used for converting grid-relative winds to Earth-relative. Then I tried to rotate the Earth-relative u and v back to model-grid-relative using Basemap’s rotate_vector, after creating a map object with the same projection properties as the model itself. Theoretically, I should arrive back at the same grid-relative u and v components. However, when I plotted the grid-relative barbs rotated back from the Earth-relative u and v using rotate_vector on my map, I found some major differences between these barbs and the original model grid-relative barbs (see attached Fig. 1. Darkred barbs=original grid-relative; Green barbs=grid-relative rotated from Earth-relative using rotate_vector-“batch”-mode). Interestingly, the grid-relative barbs rotated back from Earth-relative winds behave similarly to the strange-looking barbs that Jacob got when rotating the barbs using rotate_vector. At some points, the winds are inconveniently from the due east or due south, where the u and v printout from the model output clearly suggests otherwise.

After finding this error, I tried the alt_rotate vector code provided by Jacob, with a few modifications to allow two reference latitudes. The modifications were based on the formula on WRF-Python’s website: https://wrf-python.readthedocs.io/en/latest/user_api/generated/wrf.uvmet.html?highlight=cone#wrf.uvmet The modified alt_rotate code is also attached.

The grid-relative winds calculated from the Earth-relative winds using alt_rotate (by setting inverse=True) appear to closely match the original grid-relative winds from the model, with the largest difference only being ~1 cm/s. The plotted barbs also match well with the original grid-relative winds on the map (see attached Fig. 2. Darkred barbs=original grid-relative; Blue barbs=grid-relative rotated from Earth-relative using alt_rotate. The blue barbs cover the darkred barbs.).

I also conducted another test similar to the one that Jacob performed before where the rotation from Earth-relative back to grid-relative using rotate_vector was done element-wise for each point using nested for loops. Using this method, the derived grid-relative winds appear to be correct and match the original grid-relative winds well (see attached Fig. 3. Darkred barbs=original grid-relative; Green barbs=grid-relative rotated from Earth-relative using rotate_vector-element-wise. The green barbs cover the darkred barbs). The largest difference between the grid-relative winds rotated back from the Earth-relative winds and the original grid-relative winds is also ~1 cm/s.

Can someone familiar with Basemap provide an explanation as to why there is a difference between the resulting u and v components when rotate_vector is called in “batch” mode (supplying 2D inputs) versus called in element-wise mode?

To reproduce the above, I have attached my script and test data to this thread. I was using Basemap version 1.2.1.

Any help would be greatly appreciated, David

test_data.zip wrf_windbarb_rotation_TEST.txt alt_rotateVec.txt

[JacobCarley-NOAA]

Hi David,

In my understanding, if someone is plotting wind vectors/barbs from a WRF on a map with the same projection as the model itself, then no rotation needs to be done to the model u and v wind components since they are already grid relative

This is correct.

As for why the basemap vector rotation routine behaves strangely - unfortunately I never arrived at a satisfactory conclusion. The alternative rotation you've tried, that I also use(d), does work. For what it's worth I think most folks have moved on to cartopy https://scitools.org.uk/cartopy/docs/latest/

-Jacob

[Chun-ChihWang]

Hi Jacob,

OK. No problem. Thank you for your response.

Since all of my plots were done in Basemap, I will just stick with Basemap for now but use the alt_rotate routine whenever I need to convert Earth-relative winds back to WRF (or map) grid-relative winds.

Cheers, David

[molinav]

Thanks for the feedback, @Chun-ChihWang! The issue may stay here open until an appropriate solution is found. I do not use the barbs plotting regularly and I do not have the knowledge, but please feel free to provide any necessary feedback and I will try to take a look when I have some time slot free.

[DavidWang505]

Hi, all. I just have some updates after working on this for a while. Another quick fix to the strange-looking vectors when rotating Earth-relative winds back to grid-relative winds using Basemap's rotate_vector is to recast the input Earth-relative U and V components from float32 to float64 before passing them to the rotate_vector function. The resulting U and V are very similar (within 1 cm/s) to the native grid-relative U and V components given by the WRF model. I don't know why there is a sensitivity to the input array data type (i.e., float32 vs float64). Maybe the Basemap developers will have a better idea of why that is.

I was making a transition from Basemap to Cartopy. However, it appears that Cartopy has its own issues with vector rotation, as described in David Oven's plotting WRF model winds plotting in Cartopy example (https://www-k12.atmos.washington.edu/~ovens/wrfwinds.html) and this post in wrf-python forum here (https://groups.google.com/a/ucar.edu/g/wrfpython-talk/c/5Kcz-qFQXAE/m/l0hX0JpLGwAJ). I also checked WRF-Python's plotting example of 500 mb heights and barbs using either Cartopy or Basemap (https://wrf-python.readthedocs.io/en/latest/plot.html) and the barbs do look off in the Cartopy plot. The resulting figures should look the same between Cartopy vs. Basemap, but obviously, they don't.

There is a code to correct the vector rotation issue with Cartopy described in David Oven's example and it seems to fix the problem, if I use Cartopy to plot my barbs. Plotting the resultant U and V from the correction function in Basemap still didn't give me the same result, so it appeared that the problem was specific to Cartopy.

It seems in Cartopy (I use v0.23.0), if I want to plot the barbs in the same map projection as the WRF model itself, the code below is the way to do it, # To get the Cartopy projection instance using WRF-Python map_proj = wrf.get_cartopy(wrfin=input_wrf) # create a Cartopy projection instance, which will be the same as the WRF simulation

# A bunch of codes to read in the lon, lat, and native grid-relative U and V. I will call the WRF grid-relative U as Ua and grid-relative V as Va.

fig = plt.figure() ax = fig.add_subplot(projection=map_proj)

# First, transform the lon and lat coordinates into map projection coordinates x, y, _ = map_proj.transform_points(ccrs.PlateCarree(), lon, lat).T

# This will transpose the resultant x, y, so transpose them back to the same shape as the original lon and lat fields x=x.T y=y.T

# Plot barbs in Cartopy ax.barbs(x, y, Ua, Va, transform=map_proj) # my way

Note this is different from what other people have done. They did the below instead:

ax.barbs(lon, lat, Ua, Va, transform=ccrs.PlateCarree()) # other people's way

which I don't think is correct because when transform=ccrs.PlateCarree() is specified, I think it is telling Cartopy that the input U and V are Earth-relative, but they are not. Again, I might be wrong on this.

The plotted barbs using "my way" in Cartopy produced the same result as plotting them using Basemap in the same map projection. That's why I think "my way" is correct unless I plotted the wind barbs in Basemap incorrectly.

I will also consult with the folks in the WRF-Python forum and see what they have to say about this. At this point, I am confused about what approach I should take to plot the barbs correctly. It seems like nobody has a definitive answer yet.