updated unit tests with rotate grid script

pynemo_37.yml

@@ -10,6 +10,7 @@ dependencies:
   - pip=20.0.2
   - pandas=1.0.1
   - pytest=5.3.5
+  - xarray=0.15.0
   - pip:
     - idna==2.9
     - lxml==4.5.0

test_scripts/bdy_var_plot.py

@@ -161,7 +161,7 @@ def plot_bdy(fname, bdy_ind, bdy_dst, bdy_brk, varnam, t, rw):
                 try:
                     gdep = np.squeeze(rootgrp.variables['depthv'][:, :, :])
                 except KeyError:
-                    print 'depth variable not found'
+                    print ('depth variable not found')
 
     rootgrp.close()
 

unit_tests/namelist_unit_test.bdy

@@ -30,12 +30,12 @@
 !------------------------------------------------------------------------------
 !  grid information 
 !------------------------------------------------------------------------------
-   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/inputs/subset_coordinates.nc'
-   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/inputs/subset_coordinates.nc'
-   sn_dst_hgr = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_C/mesh_hgr_zps.nc'
-   sn_dst_zgr = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_C/mesh_zgr_zps.nc'
-   sn_src_msk = '/Users/thopri/Projects/PyNEMO/inputs/subset_bathy.nc'
-   sn_bathy   = 'http://opendap4gws.jasmin.ac.uk/thredds/noc_msm/dodsC/pynemo_grid_C/NNA_R12_bathy_meter_bench.nc'
+   sn_src_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/src_coordinates.nc'
+   sn_src_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/src_coordinates.nc'
+   sn_dst_hgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/dst_hgr_zps.nc'
+   sn_dst_zgr = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/dst_zgr_zps.nc'
+   sn_src_msk = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/src_bathy.nc'
+   sn_bathy   = '/Users/thopri/Projects/PyNEMO/unit_tests/test_data/dst_bathy.nc'
 
 !------------------------------------------------------------------------------
 !  I/O 

unit_tests/resample_netcdf.py

+# -*- coding: utf-8 -*-
+"""
+Resample netcdf code, used for making small low res testing datasets.
+"""
+
+import xarray as xr
+import numpy as np
+
+# subset X Y dims netcdf
+ds = xr.open_dataset('unit_tests/test_data/NNA_R12_bathy_meter_bench.nc')
+
+new_y = np.linspace(ds.y[0],ds.y[-1],100)
+
+new_x = np.linspace(ds.x[0],ds.x[-1],85)
+
+dsi = ds.interp(x=new_x,y=new_y)
+
+dsi.to_netcdf('unit_tests/test_data/dst_bathy.nc')
+
+# subset lat lon dims netcdf
+
+ds = xr.open_dataset('unit_tests/test_data/subset_bathy.nc')
+
+new_lon = np.linspace(ds.longitude[0],ds.longitude[-1],45)
+
+new_lat = np.linspace(ds.latitude[0],ds.latitude[-1],30)
+
+dsi = ds.interp(latitude=new_lat,longitude=new_lon)
+
+dsi.to_netcdf('unit_tests/test_data/resample_bathy.nc')
+
+
+# drop unneeded variables
+dsd = xr.Dataset.drop_vars(ds,'zos')
+
+# remove unwanted time
+ds_sel = dsd.sel({'time':'01-01-2017'})
+

unit_tests/rotate_pts.py

+"""
+Source: https://gist.github.com/LyleScott/d17e9d314fbe6fc29767d8c5c029c362
+
+Adapted Function to rotate NEMO grid for unit testing.
+
+"""
+
+from math import cos,sin,radians
+import numpy as np
+import xarray as xr
+import matplotlib.pyplot as plt
+import cartopy
+import cartopy.crs as ccrs
+
+def rotate_around_point(lat_in,lon_in, radians , origin=(0, 0)):
+    """Rotate a point around a given point.
+    """
+    new_lon = np.zeros(np.shape(lon_in))
+    new_lat = np.zeros(np.shape(lat_in))
+
+    offset_lat, offset_lon = origin
+    cos_rad = cos(radians)
+    sin_rad = sin(radians)
+
+    for indx in range(np.shape(lat_in)[0]):
+        for indy in range(np.shape(lon_in)[1]):
+
+            adjusted_lat = (lat_in[indx,indy] - offset_lat)
+            adjusted_lon = (lon_in[indx,indy] - offset_lon)
+            new_lat[indx,indy] = offset_lat + cos_rad * adjusted_lat + -sin_rad * adjusted_lon
+            new_lon[indx,indy] = offset_lon + sin_rad * adjusted_lat + cos_rad * adjusted_lon
+
+    return new_lat, new_lon
+
+def plot_grids(lat_in,lon_in,new_lat,new_lon,src_lat,src_lon):
+
+    maxlat = 72
+    minlat = 32
+    maxlon = 18
+    minlon = -28
+
+    plt.figure(figsize=[18, 8])  # a new figure window
+
+    ax = plt.subplot(121, projection=ccrs.PlateCarree())  # specify (nrows, ncols, axnum)
+    ax.set_extent([minlon,maxlon,minlat,maxlat],crs=ccrs.PlateCarree()) #set extent
+    ax.set_title('Original Grid', fontsize=20,y=1.05) #set title
+    ax.add_feature(cartopy.feature.LAND, zorder=0)  # add land polygon
+    ax.add_feature(cartopy.feature.COASTLINE, zorder=10)  # add coastline polyline
+    gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=2, color='black', alpha=0.5, linestyle='--', draw_labels=True)
+
+    ax1 = plt.subplot(122, projection=ccrs.PlateCarree())  # specify (nrows, ncols, axnum)
+    ax1.set_extent([minlon,maxlon,minlat,maxlat],crs=ccrs.PlateCarree()) #set extent
+    ax1.set_title('Rotated Grid', fontsize=20,y=1.05) #set title
+    ax1.add_feature(cartopy.feature.LAND, zorder=0)  # add land polygon
+    ax1.add_feature(cartopy.feature.COASTLINE, zorder=10)  # add coastline polyline
+    gl1 = ax1.gridlines(crs=ccrs.PlateCarree(), linewidth=2, color='black', alpha=0.5, linestyle='--', draw_labels=True)
+
+    src_lon = np.tile(src_lon, (np.shape(src_lat)[0], 1))
+    src_lat = np.tile(src_lat, (np.shape(src_lon)[1], 1))
+    src_lat = np.rot90(src_lat)
+
+    ax.plot(lon_in, lat_in, color='blue', marker='.')
+    ax.plot(src_lon,src_lat, color='green', marker='.')
+    ax1.plot(new_lon,new_lat, color='red', marker='.')
+    ax1.plot(src_lon,src_lat, color='green',marker='.')
+    plt.subplots_adjust(left=0.01, right=1, top=0.9, bottom=0.05,wspace=0.01)
+
+    plt.show()
+
+
+def _main():
+    ds = xr.open_dataset('unit_tests/test_data/dst_hgr_zps.nc')
+    src = xr.open_dataset('unit_tests/test_data/src_coordinates.nc')
+    rot = 45
+    theta = radians(rot)
+    origin = (54, -6)
+    new_lat,new_lon = rotate_around_point(ds.nav_lat.values[:],ds.nav_lon.values[:],theta,origin)
+    plot_grids(ds.nav_lat.values[:],ds.nav_lon.values[:],new_lat,new_lon,src.latitude.values[:],src.longitude.values[:])
+    ds.nav_lat.values[:], ds.nav_lon.values[:] = new_lat,new_lon
+    ds.to_netcdf('unit_tests/test_data/rot_'+str(rot)+'_dst_hgr_zps.nc')
+    ds.close()
+    src.close()
+
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    _main()
\ No newline at end of file