diff --git a/unit_tests/coord_gen.py b/unit_tests/coord_gen.py
index 2408da9baa522d63fdb982d346f7be69081820b9..c05e3bf082d6bce6e4e57a7dcd4c8b3fc45fa8b9 100644
--- a/unit_tests/coord_gen.py
+++ b/unit_tests/coord_gen.py
@@ -10,7 +10,10 @@ example grid set up
import numpy as np
from netCDF4 import Dataset
-def set_hgrid(dx, dy, jpi, jpj, zoffx, zoffy):
+def set_hgrid(dx, dy, jpi, jpj, zoffx=0, zoffy=0,sf=1):
+ if sf > 1:
+ jpi = int(jpj - (jpi/sf))+1
+ jpj = int(jpj -(jpj/sf))+1
# Set grid positions [km]
latt = np.zeros((jpi, jpj))
lont = np.zeros((jpi, jpj))
@@ -21,6 +24,10 @@ def set_hgrid(dx, dy, jpi, jpj, zoffx, zoffy):
lonf = np.zeros((jpi, jpj))
latf = np.zeros((jpi, jpj))
+
+
+
+
for i in range(0, jpi):
lont[i, :] = zoffx * dx * 1.e-3 + dx * 1.e-3 * np.float(i)
lonu[i, :] = zoffx * dx * 1.e-3 + dx * 1.e-3 * (np.float(i) + 0.5)
diff --git a/unit_tests/rotate_pts.py b/unit_tests/rotate_pts.py
index c673f5f7d49afae72b58dd622a4aacad91e15a28..c2a529336b4b5774da13d543d9d9cb21c6efdd25 100644
--- a/unit_tests/rotate_pts.py
+++ b/unit_tests/rotate_pts.py
@@ -40,7 +40,7 @@ def rotate_around_point(lat_in,lon_in, radians , origin=(0, 0)):
return new_lat, new_lon
-def plot_grids(lat_in,lon_in,new_lat,new_lon,src_lat,src_lon):
+def plot_grids(lat_in,lon_in,new_lat,new_lon,off_lat,off_lon,src_lat,src_lon):
# define lat and lon extents (define in future using input values?)
#maxlat = 72
@@ -48,16 +48,16 @@ def plot_grids(lat_in,lon_in,new_lat,new_lon,src_lat,src_lon):
#maxlon = 18
#minlon = -28
- plt.figure(figsize=[18, 8]) # a new figure window
+ plt.figure(figsize=[18, 18]) # a new figure window
- ax = plt.subplot(121)#, projection=ccrs.PlateCarree()) # specify (nrows, ncols, axnum)
+ ax = plt.subplot(221)#, 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 = plt.subplot(222)#, 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
@@ -69,11 +69,16 @@ def plot_grids(lat_in,lon_in,new_lat,new_lon,src_lat,src_lon):
#src_lat = np.tile(src_lat, (np.shape(src_lon)[1], 1))
#src_lat = np.rot90(src_lat)
+ ax2 = plt.subplot(223)
+ ax3 = plt.subplot(224)
+
# plot lat and lon for all grids
ax.plot(lon_in, lat_in, color='blue', marker='.', linestyle="")
ax.plot(src_lon,src_lat, color='green', marker='o', linestyle="")
ax1.plot(new_lon,new_lat, color='red', marker='.', linestyle="")
ax1.plot(src_lon,src_lat, color='green',marker='o',linestyle="")
+ ax2.plot(off_lon,off_lat, color='yellow',marker='.',linestyle="")
+ ax2.plot(src_lon, src_lat, color='green', marker='o', linestyle="")
# tweak margins of subplots as tight layout doesn't work
plt.subplots_adjust(left=0.01, right=1, top=0.9, bottom=0.05,wspace=0.01)
@@ -86,15 +91,13 @@ def _main():
dy = 1000 # units in Km
jpi = 15
jpj = 15
- zoffx = 0
- zoffy = 0
jpk = 10
max_dep = 100
min_dep = 10
- z_end_dim = 10.0
+ z_end_dim = 1
h_fname = 'unit_tests/test_data/test_src__hgr_zps.nc'
z_fname = 'unit_tests/test_data/test_src_zgr_zps.nc'
- grid_h1 = cg.set_hgrid(dx,dy,jpi,jpj,zoffx,zoffy)
+ grid_h1 = cg.set_hgrid(dx,dy,jpi,jpj)
grid_z1 = cg.set_zgrid(grid_h1,jpk,max_dep,min_dep,z_end_dim)
write_coord_H = cg.write_coord_H(h_fname,grid_h1)
write_coord_Z = cg.write_coord_Z(z_fname,grid_h1,grid_z1)
@@ -111,10 +114,11 @@ def _main():
jpk = 10
max_dep = 100
min_dep = 10
- z_end_dim = 10.0
+ z_end_dim = 1
+ sf = 10
h_fname = 'unit_tests/test_data/test_dst_hgr_zps.nc'
z_fname = 'unit_tests/test_data/test_dst_zgr_zps.nc'
- grid_h2 = cg.set_hgrid(dx,dy,jpi,jpj,zoffx,zoffy)
+ grid_h2 = cg.set_hgrid(dx,dy,jpi,jpj,zoffx,zoffy,sf)
grid_z2 = cg.set_zgrid(grid_h2,jpk,max_dep,min_dep,z_end_dim)
write_coord_H = cg.write_coord_H(h_fname,grid_h2)
write_coord_Z = cg.write_coord_Z(z_fname,grid_h2,grid_z2)
@@ -127,9 +131,33 @@ def _main():
origin = (7,7)
- new_lat,new_lon = rotate_around_point(grid_h2['latt'],grid_h2['lont'],theta,origin)
+ rot_lat,rot_lon = rotate_around_point(grid_h2['latt'],grid_h2['lont'],theta,origin)
+
+ #offset coords
+ dx = 100 # units in km
+ dy = 100 # units in Km
+ jpi = 100
+ jpj = 100
+ zoffx = 25
+ zoffy = 25
+ jpk = 10
+ max_dep = 100
+ min_dep = 10
+ z_end_dim = 1
+ sf = 10
+ h_fname = 'unit_tests/test_data/test_dst_offset_hgr_zps.nc'
+ z_fname = 'unit_tests/test_data/test_dst_offset_zgr_zps.nc'
+ grid_h3 = cg.set_hgrid(dx,dy,jpi,jpj,zoffx,zoffy,sf)
+ grid_z3 = cg.set_zgrid(grid_h2,jpk,max_dep,min_dep,z_end_dim)
+ write_coord_H = cg.write_coord_H(h_fname,grid_h3)
+ write_coord_Z = cg.write_coord_Z(z_fname,grid_h3,grid_z3)
+ if write_coord_H + write_coord_Z == 0:
+ print("Success!")
+
+
# plot orginal, rotatated and source lat and lon
- plot_grids(grid_h2['latt'],grid_h2['lont'],new_lat,new_lon,grid_h1['latt'],grid_h1['lont'])
+ plot_grids(grid_h2['latt'],grid_h2['lont'],rot_lat,rot_lon,grid_h3['latt'], \
+ grid_h3['lont'],grid_h1['latt'],grid_h1['lont'])
# save new lat lon to dataset
#ds.nav_lat.values[:], ds.nav_lon.values[:] = new_lat,new_lon
# save dataset to netcdf