updated unit test grid generation

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)

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