makegrid not working with npstere projection
bward-mcgill opened this issue · 2 comments
Hello,
I'm trying to reproduce the transform_vector() function. Basically, I want to plot my wind vector rotated and interpolated in a north polar stereographic projection. Here's the source code :
def transform_vector(self,uin,vin,lons,lats,nx,ny,returnxy=False,checkbounds=False,order=1,masked=False):
"""
Rotate and interpolate a vector field (``uin,vin``) from a
lat/lon grid with longitudes = ``lons`` and latitudes = ``lats``
to a ``ny`` by ``nx`` map projection grid.
...
"""
delon = lons[1:]-lons[0:-1]
delat = lats[1:]-lats[0:-1]
if min(delon) < 0. or min(delat) < 0.:
raise ValueError('lons and lats must be increasing!')
# check that lons in -180,180 for non-cylindrical projections.
if self.projection not in _cylproj:
lonsa = np.array(lons)
count = np.sum(lonsa < -180.00001) + np.sum(lonsa > 180.00001)
if count > 1:
raise ValueError('grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)')
# allow for wraparound point to be outside.
elif count == 1 and math.fabs(lons[-1]-lons[0]-360.) > 1.e-4:
raise ValueError('grid must be shifted so that lons are monotonically increasing and fit in range -180,+180 (see shiftgrid function)')
lonsout, latsout, x, y = self.makegrid(nx,ny,returnxy=True)
# interpolate to map projection coordinates.
uin = interp(uin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
vin = interp(vin,lons,lats,lonsout,latsout,checkbounds=checkbounds,order=order,masked=masked)
# rotate from geographic to map coordinates.
return self.rotate_vector(uin,vin,lonsout,latsout,returnxy=returnxy)So my understanding is that the function is quite simple. There is some verifications at the beginning in order to make sure that the longitude goes to -180 to 180 and that both longitude and latitude are well ordered (I guess that this is important for interp()). After that it calls makegrid() to create a nice equally spaced grid in the projection. After that it interpolates with interp() . Finally it rotate the vectors in the projection with rotate_vector. I was able to use the function with a regular rectilinear grid. In such case, the a grid can be fully represented only by two 1d vectors (the latitude is constant when we are moving on a longitude line) so I can specify 1d vector as input for transform_vector() and it works. Here's what it looks like :
However, I am working with a rotated pole curvilinear grid : my grid need to be represented by two 2d matrix since the latitude is changing when we are moving on a constant longitude line. Therefore, I cannot use transform_vector() anymore. From what I saw, the only problem comes in interp() function (it can only take a 1d array for lons and lats). The rest should still work with 2d matix (rotate_vector() still work with 2d matrix and it should not make any difference for makegrid()). Therefore, I implemented a function that is really similar to transform_vector() but it uses a different method for the interpolation.
def create_wind_projection(path_a, file_a, path_g, file_g, basem, wind_var, nlon, nlat):
import pandas as pd
import xarray as xr
import os
#Read file wind a create a new netCDF file with only winds in it
file_data=path_a+"/"+file_a
name=file_a[:-3]
file_wind=path_a+"/"+name+"_wnd.nc"
for wnd in wind_var:
if wnd=='uwnd' or wnd=='wndewd':
uwind=read_data(path_a, [file_a], wnd, nlon, nlat)*1.9434
elif wnd=='vwnd' or wnd=='wndnwd':
vwind=read_data(path_a, [file_a], wnd, nlon, nlat)*1.9434
df = pd.DataFrame()
dsWind=df.to_xarray()
dsWind['uwnd']=(("latitude","longitude"),uwind)
dsWind['vwnd']=(("latitude","longitude"),vwind)
dsWind.to_netcdf(file_wind)
#Create the grid, basem is the projection that I'm using.
lonsout, latsout=basem.makegrid(41,41)
#Create instruction to interpolate with CDO
new_grid_info=path_a+"/remap_griddes.info"
np.savetxt(path_a+"/lat_remap.dat", latsout[:,1], fmt='%.2f')
np.savetxt(path_a+"/lon_remap.dat", lonsout[1], fmt='%.2f')
os.system("echo 'gridtype = lonlat' > "+new_grid_info)
os.system("echo 'gridsize = '"+str(latsout.shape[0]*latsout.shape[1])+" >> "+new_grid_info)
os.system("echo 'xsize = '"+str(latsout.shape[0])+" >> "+new_grid_info)
os.system("echo 'ysize = '"+str(latsout.shape[1])+" >> "+new_grid_info)
os.system("echo 'xvals =' >> "+new_grid_info+" && cat "+path_a+"/lon_remap.dat >>"+new_grid_info)
os.system("echo 'yvals =' >> "+new_grid_info+" && cat "+path_a+"/lat_remap.dat >>"+new_grid_info)
#Interpolate
os.system('/opt/cdo/bin/cdo remapbil,'+new_grid_info+" "+file_wind+" "+file_wind+"_interp>/dev/null $
os.system('rm -f '+file_wind)
os.system('mv '+file_wind+"_interp"+" "+file_wind)
dsWind_new=xr.open_dataset(file_wind)
dsWind_new=dsWind_new[wind_var]
#Rotate in the projection
for wnd in wind_var:
if wnd=='uwnd' or wnd=='wndewd':
uwind_new=read_data(path_a, [name+"_wnd.nc"], wnd, 41, 41)
elif wnd=='vwnd' or wnd=='wndnwd':
vwind_new=read_data(path_a, [name+"_wnd.nc"], wnd, 41, 41)
lat_new=np.nan_to_num(dsWind_new[['lat']]['lat'].values.squeeze())
lon_new=np.nan_to_num(dsWind_new[['lon']]['lon'].values.squeeze())
urot,vrot,xx,yy = basem.rotate_vector(uwind_new,vwind_new,lon_new,lat_new,returnxy=True)
return urot,vrot,xx,yyThe function is working except when I'm using the a stereographic projection for basem. For example, when I use :
basem = Basemap(projection='npstere', boundinglat=35,lon_0=270, resolution='l')
Here's what I get :
Here's what it gives me for latsout and lonsout :
lonsout = -135.00 -133.45 -131.82 -130.10 -128.29 -126.38 -124.38 ...-43.53
latsout = 18.64 20.06 21.47 22.86 24.23 25.57 26.88 28.14 20.06 ... 18.64
In don't understand why because I'm calling makegrid() in the same way that it is called in transform_vector() but for some reason it doesn't work for me. What confuse me even more is that when I use the mercator projection the makegrid() works perfectly. For example when I use basem = Basemap(projection='merc',llcrnrlat=-80,urcrnrlat=80, llcrnrlon=-180,urcrnrlon=180,lat_ts=20,resolution='c')
I now get :
I hope that my problem is clear. If anyone has any idea on why the makegrid() function doesn't work when I'm calling it in my create_wind_projection() function for a npstere projection, please let me know !!
@bward-mcgill Did you manage to solve the problem yourself? If not, I would like to keep the issue open to remember that it still needs a fix. I try to do my best but I cannot find always enough time to manage all the issues. Sorry!
Hello @molinav, I figured that the problem didn't come from the makegrid() function, but from an other part of the function that I implemented. I assumed that the makegrid was returning a rectilinear grid, but I was wrong since an equally spaced grid in a stereographic projection is in fact curvilinear. Therefore, I add to change a little bit the way I was interpolating my field and now I fixed it.
That being said, I still think that I could be nice to add the possibility to use the transform_vector with 2-D longitudes and latitudes so it could work for people that are working with non-rectilinear grid.
Thanks !!