diff --git a/pynemo_37.yml b/pynemo_37.yml index 8d24b7fa6bd0fd15d2968f5c1de609c1e0129d2d..21ab2e875aa28671e7e74823bf9fbd2dfc2067aa 100644 --- a/pynemo_37.yml +++ b/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 diff --git a/test_scripts/bdy_var_plot.py b/test_scripts/bdy_var_plot.py index 763688c73db1f5fe52acdac2b8a4f1f2f4b0714e..be1e8136ca3fe4fc2f585a20e90b59b8f6c4ebcd 100644 --- a/test_scripts/bdy_var_plot.py +++ b/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() diff --git a/unit_tests/namelist_unit_test.bdy b/unit_tests/namelist_unit_test.bdy index 48746cf4ec6c186c82b1b5bb34e90e11daaf64cc..b5368956db59adc18f0f74e8e048296b5d227df3 100644 --- a/unit_tests/namelist_unit_test.bdy +++ b/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 diff --git a/unit_tests/resample_netcdf.py b/unit_tests/resample_netcdf.py new file mode 100644 index 0000000000000000000000000000000000000000..8c352ba3c78d69dc8b2599a0f39c1547ad5241c1 --- /dev/null +++ b/unit_tests/resample_netcdf.py @@ -0,0 +1,38 @@ +# -*- 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'}) + diff --git a/unit_tests/rotate_pts.py b/unit_tests/rotate_pts.py new file mode 100644 index 0000000000000000000000000000000000000000..494d989e92984e5690c0c7fbb8a6222cac784d76 --- /dev/null +++ b/unit_tests/rotate_pts.py @@ -0,0 +1,92 @@ +""" +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