From 42e7c5a025d1d1d482e89729adf563e707e0bc36 Mon Sep 17 00:00:00 2001
From: thopri <thopri@144-148.noc.soton.ac.uk>
Date: Wed, 11 Mar 2020 11:56:07 +0000
Subject: [PATCH] updated rotate script
---
unit_tests/rotate_pts.py | 33 +++++++++++++++++++++------------
1 file changed, 21 insertions(+), 12 deletions(-)
diff --git a/unit_tests/rotate_pts.py b/unit_tests/rotate_pts.py
index 494d989..3541b65 100644
--- a/unit_tests/rotate_pts.py
+++ b/unit_tests/rotate_pts.py
@@ -1,10 +1,11 @@
"""
Source: https://gist.github.com/LyleScott/d17e9d314fbe6fc29767d8c5c029c362
-Adapted Function to rotate NEMO grid for unit testing.
+Adapted Function to rotate NEMO grid for unit testing. rot variable defines how much rotation
+occurs, the results are then plotted in a figure showing original (blue points) and rotated (red) grids.
+The source coordinate grid is also plotted (green).
"""
-
from math import cos,sin,radians
import numpy as np
import xarray as xr
@@ -12,19 +13,23 @@ import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
+# TODO: add in auto max and min lat and lon
+# TODO: remove hard coded file names and directories.
+
def rotate_around_point(lat_in,lon_in, radians , origin=(0, 0)):
"""Rotate a point around a given point.
"""
+ # create empty lat and lon arrays that match input
new_lon = np.zeros(np.shape(lon_in))
new_lat = np.zeros(np.shape(lat_in))
-
+ # extract origin lat and lon as offset
offset_lat, offset_lon = origin
cos_rad = cos(radians)
sin_rad = sin(radians)
-
+ # cycle through lat and lon arrays and calcule new lat and lon values based on rotation
+ # and origin values
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
@@ -34,6 +39,7 @@ def rotate_around_point(lat_in,lon_in, radians , origin=(0, 0)):
def plot_grids(lat_in,lon_in,new_lat,new_lon,src_lat,src_lon):
+ # define lat and lon extents (define in future using input values?)
maxlat = 72
minlat = 32
maxlon = 18
@@ -55,38 +61,41 @@ def plot_grids(lat_in,lon_in,new_lat,new_lon,src_lat,src_lon):
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)
+ # tile 1D lat and lon to 2D arrays for plotting (src lat and lon only)
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)
+ # plot lat and lon for all grids
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='.')
+ # 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)
plt.show()
def _main():
+ #open source (parent) and destination (child) grids
ds = xr.open_dataset('unit_tests/test_data/dst_hgr_zps.nc')
src = xr.open_dataset('unit_tests/test_data/src_coordinates.nc')
+ # set rotation and origin point
rot = 45
theta = radians(rot)
origin = (54, -6)
+ # generate new lat and lon
new_lat,new_lon = rotate_around_point(ds.nav_lat.values[:],ds.nav_lon.values[:],theta,origin)
+ # plot orginal, rotatated and source lat and lon
plot_grids(ds.nav_lat.values[:],ds.nav_lon.values[:],new_lat,new_lon,src.latitude.values[:],src.longitude.values[:])
+ # save new lat lon to dataset
ds.nav_lat.values[:], ds.nav_lon.values[:] = new_lat,new_lon
+ # save dataset to netcdf
ds.to_netcdf('unit_tests/test_data/rot_'+str(rot)+'_dst_hgr_zps.nc')
+ # close data files
ds.close()
src.close()
-
-
-
-
-
-
-
if __name__ == '__main__':
_main()
\ No newline at end of file
--
GitLab