makegrid not working with npstere projection

[bward-mcgill]

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,yy

The 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 !!

[molinav]

@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!

[bward-mcgill]

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 !!