diff --git a/fmsgridtools/make_hgrid/gnomonic_grid.py b/fmsgridtools/make_hgrid/gnomonic_grid.py index 182f910..d9b320c 100644 --- a/fmsgridtools/make_hgrid/gnomonic_grid.py +++ b/fmsgridtools/make_hgrid/gnomonic_grid.py @@ -24,6 +24,7 @@ def make( output_length_angle: bool, do_schmidt: bool, do_cube_transform: bool, + transpose: bool, verbose: bool, ): @@ -183,5 +184,6 @@ def make( conformal=conformal, out_halo=out_halo, output_length_angle=output_length_angle, + transpose=transpose, verbose=verbose, ) diff --git a/fmsgridtools/make_hgrid/hgridobj.py b/fmsgridtools/make_hgrid/hgridobj.py index 8d3d78f..def2fbb 100644 --- a/fmsgridtools/make_hgrid/hgridobj.py +++ b/fmsgridtools/make_hgrid/hgridobj.py @@ -1,17 +1,18 @@ -import sys import ctypes +import sys + import numpy as np -from numpy.typing import NDArray import xarray as xr +from numpy.typing import NDArray from fmsgridtools.shared.gridtools_utils import get_provenance_attrs from fmsgridtools.shared.gridobj import GridObj def fill_cubic_grid_halo( - nx: int, - ny: int, - halo: int, - data: NDArray[np.float64], + nx: int, + ny: int, + halo: int, + data: NDArray[np.float64], data1_all: NDArray[np.float64], data2_all: NDArray[np.float64], tile: int, @@ -22,7 +23,7 @@ def fill_cubic_grid_halo( nyp = ny + joff nxph = nx + ioff + 2*halo nyph = ny + joff + 2*halo - + for i in range(nxph*nyph): data[i] = -9999. @@ -30,9 +31,9 @@ def fill_cubic_grid_halo( for j in range (1, nyp+1): for i in range(1, nxp+1): data[j*nxph+i] = data1_all[tile*nxp*nyp+(j-1)*nxp+(i-1)] - + ntiles = 6 - + if tile%2 == 1: lw = (tile+ntiles-1)%ntiles le = (tile+ntiles+2)%ntiles @@ -216,9 +217,10 @@ def write_out_hgrid( conformal: bool=True, out_halo: int=0, output_length_angle: bool=True, + transpose: bool=False, verbose: bool=False, ): - + var_dict={} pos_c = 0 pos_e = 0 @@ -267,7 +269,7 @@ def write_out_hgrid( ) var_dict['arcx'] = arcx - + """define dimension""" nx = self.nxl[n] ny = self.nyl[n] @@ -283,20 +285,20 @@ def write_out_hgrid( if n > 0: print(f"[INFO] XARRAY: n: {n} x[0]: {self.x[pos_c]} x[-1]: {self.x[pos_c-1]} x[-2]: {self.x[pos_c-2]} x[-3]: {self.x[pos_c-3]} x[-4]: {self.x[pos_c-4]} x[-5]: {self.x[pos_c-5]} x[-10]: {self.x[pos_c-10]}", file=sys.stderr) x = xr.DataArray( - data=self.x[pos_c:pos_c+nyp*nxp].reshape((nyp,nxp)), + data=self.x[pos_c:pos_c+nyp*nxp].reshape((nyp,nxp)).T if transpose else self.x[pos_c:pos_c+nyp*nxp].reshape((nyp,nxp)), dims=["nyp", "nxp"], attrs=dict( - units="degree_east", + units="degree_east", standard_name="geographic_longitude", ) ) var_dict['x'] = x y = xr.DataArray( - data=self.y[pos_c:pos_c+nyp*nxp].reshape((nyp, nxp)), + data=self.y[pos_c:pos_c+nyp*nxp].reshape((nyp, nxp)).T if transpose else self.y[pos_c:pos_c+nyp*nxp].reshape((nyp, nxp)), dims=["nyp", "nxp"], attrs=dict( - units="degree_north", + units="degree_north", standard_name="geographic_latitude", ) ) @@ -317,7 +319,7 @@ def write_out_hgrid( data=self.dx[pos_n:pos_n+nyp*nx].reshape((nyp, nx)), dims=["nyp", "nx"], attrs=dict( - units="meters", + units="meters", standard_name="grid_edge_x_distance", ) ) @@ -327,7 +329,7 @@ def write_out_hgrid( data=self.dy[pos_e:pos_e+ny*nxp].reshape((ny, nxp)), dims=["ny", "nxp"], attrs=dict( - units="meters", + units="meters", standard_name="grid_edge_y_distance", ) ) @@ -354,7 +356,7 @@ def write_out_hgrid( ) var_dict['angle_dy'] = angle_dy else: - if grid_type is not "gnomonic_ed": + if grid_type != "gnomonic_ed": raise RuntimeError("make_hgrid: out_halo > 0, only working for grid_type = 'gnomonic_ed'") if verbose: @@ -363,10 +365,10 @@ def write_out_hgrid( tmp_x = np.zeros(shape=(nxp+2*out_halo)*(nyp+2*out_halo), dtype=np.float64) fill_cubic_grid_halo(nx, ny, out_halo, tmp_x, self.x, self.x, n, 1, 1) x = xr.DataArray( - data=tmp_x.reshape((nyp+2*out_halo,nxp+2*out_halo)), + data=tmp_x.reshape((nyp+2*out_halo,nxp+2*out_halo)).T if transpose else tmp_x.reshape((nyp+2*out_halo,nxp+2*out_halo)), dims=["nyp", "nxp"], attrs=dict( - units="degree_east", + units="degree_east", standard_name="geographic_longitude", _FillValue=-9999., ) @@ -376,10 +378,10 @@ def write_out_hgrid( tmp_y = np.zeros(shape=(nxp+2*out_halo)*(nyp+2*out_halo), dtype=np.float64) fill_cubic_grid_halo(nx, ny, out_halo, tmp_y, self.y, self.y, n, 1, 1) y = xr.DataArray( - data=tmp_y.reshape((nyp+2*out_halo, nxp+2*out_halo)), + data=tmp_y.reshape((nyp+2*out_halo, nxp+2*out_halo)).T if transpose else tmp_y.reshape((nyp+2*out_halo, nxp+2*out_halo)), dims=["nyp", "nxp"], attrs=dict( - units="degree_north", + units="degree_north", standard_name="geographic_latitude", _FillValue = -9999., ) @@ -406,7 +408,7 @@ def write_out_hgrid( data=tmp_dx.reshape((nyp+2*out_halo, nx+2*out_halo)), dims=["nyp", "nx"], attrs=dict( - units="meters", + units="meters", standard_name="grid_edge_x_distance", _FillValue=-9999., ) @@ -419,7 +421,7 @@ def write_out_hgrid( data=tmp_dy.reshape((ny+2*out_halo, nxp+2*out_halo)), dims=["ny", "nxp"], attrs=dict( - units="meters", + units="meters", standard_name="grid_edge_y_distance", _FillValue=-9999., ) @@ -453,6 +455,17 @@ def write_out_hgrid( ) var_dict['angle_dy'] = angle_dy + if transpose: + ordering = xr.DataArray( + ["row_major"], + ) + else: + ordering = xr.DataArray( + ["column_major"], + ) + + var_dict['ordering'] = ordering + nx = self.nxl[n] ny = self.nyl[n] nxp = nx + 1 @@ -470,7 +483,7 @@ def write_out_hgrid( if verbose: print(f"About to close {outfile}") - prov_attrs = get_provenance_attrs(great_circle_algorithm=True) + prov_attrs = get_provenance_attrs(great_circle_algorithm=True) dataset = xr.Dataset( data_vars=var_dict @@ -539,5 +552,3 @@ def make_gridobj(self) -> "GridObj": else: dataset=self.dataset return GridObj(dataset=dataset) - - diff --git a/fmsgridtools/make_hgrid/make_hgrid.py b/fmsgridtools/make_hgrid/make_hgrid.py index 7c36927..76fb5a8 100644 --- a/fmsgridtools/make_hgrid/make_hgrid.py +++ b/fmsgridtools/make_hgrid/make_hgrid.py @@ -225,6 +225,12 @@ def lonlat( be set: --stretch_factor, --target_lon, and --target_lat. """, ) +@click.option( + "--transpose", + is_flag=True, + default=False, + help="This flag will transpose the x(lon) and y(lat) output arrays" +) @click.option( "--verbose", is_flag=True, @@ -250,6 +256,7 @@ def gnomonic( output_length_angle: bool, do_schmidt: bool, do_cube_transform: bool, + transpose: bool, verbose: bool ): gnomonic_grid.make( @@ -272,5 +279,6 @@ def gnomonic( output_length_angle=output_length_angle, do_schmidt=do_schmidt, do_cube_transform=do_cube_transform, + transpose=transpose, verbose=verbose, )