Skip to content

Commit 96699ab

Browse files
committed
Merge branch 'feature/MosaicObj' of https://github.com/hderry/FMSgridtools into feature/MosaicObj
2 parents 45a7f9e + 7113843 commit 96699ab

8 files changed

Lines changed: 530 additions & 0 deletions

File tree

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
import click
2+
import solo_mosaic
3+
import regional_mosaic
4+
5+
6+
ATMOS_MOSAIC_HELP = "specify the atmosphere mosaic information \
7+
This file contains list of tile files which specify \
8+
the grid information for each tile. Each grid is \
9+
required to be logically rectangular grid. The \
10+
file name cannot be 'mosaic.nc'"
11+
12+
TILE_FILE_HELP = "Grid file name of all tiles in the mosaic. \
13+
The file name should be the \
14+
relative path (exclude the absolute file path) \
15+
The absolute file path will be dir/tile_file. \
16+
the default for tile_file will be'horizontal_grid.tile#.nc'"
17+
18+
19+
20+
mosaic_name = click.option('--mosaic_name',
21+
default='mosaic',
22+
help="mosaic name; The output file will be mosaic_name.nc.")
23+
24+
sea_level = click.option('--sea_level', type=int, default=0)
25+
26+
ocean_topog = click.option('--ocean_topog',
27+
type=click.Path(exists=True))
28+
29+
30+
31+
@click.group()
32+
def make_mosaic():
33+
'''MAKE MOSAIC CLI'''
34+
35+
@make_mosaic.command()
36+
@mosaic_name
37+
@click.option("--num_tiles",
38+
type=int,
39+
required=True,
40+
help="Number of tiles in the mosaic.")
41+
42+
@click.option('--dir_name',
43+
type=click.Path(exists=True, file_okay=False),
44+
help="the directory that contains all the tile grid files." )
45+
46+
@click.option('--tile_file', '-f',
47+
multiple=True, type=click.Path(exists=True, file_okay=True),
48+
help=TILE_FILE_HELP)
49+
50+
@click.option('--periodx',
51+
default=0,
52+
help = "Specify the period in x-direction of mosaic. \
53+
Default value is 0 (not periodic)")
54+
55+
@click.option('--periody', default=0,
56+
help="Specify the period in y-direction of mosaic. \
57+
Default value is 0 (not periodic)")
58+
59+
def solo(num_tiles,
60+
dir_name,
61+
mosaic_name,
62+
tile_file,
63+
periodx,
64+
periody):
65+
66+
solo_mosaic.make(num_tiles,
67+
dir_name,
68+
mosaic_name,
69+
tile_file,
70+
periodx,
71+
periody)
72+
73+
74+
@make_mosaic.command()
75+
@click.option('--global_mosaic',
76+
type=click.Path(exists=True),
77+
help="global_mosaic Specify the mosaic file for the global grid.")
78+
79+
@click.option('--regional_file',
80+
type=click.Path(exists=True),
81+
help="regional_file Specify the regional model output file.")
82+
83+
def regional(global_mosaic,
84+
regional_file):
85+
86+
regional_mosaic.make(global_mosaic,
87+
regional_file)
88+
89+
@make_mosaic.command()
90+
@mosaic_name
91+
@sea_level
92+
@ocean_topog
93+
@click.option('--input_mosaic',
94+
type=click.Path(exists=True))
95+
96+
def quick(input_mosaic,
97+
mosaic_name,
98+
ocean_topog,
99+
sea_level,
100+
land_frac_file,
101+
land_frac_field):
102+
pass
103+
104+
@make_mosaic.command()
105+
@mosaic_name
106+
@sea_level
107+
@ocean_topog
108+
@click.option('--atmos_mosaic',
109+
type=click.Path(exists=True),
110+
help = ATMOS_MOSAIC_HELP)
111+
112+
@click.option('--ocean_mosaic',
113+
type=click.Path(exists=True))
114+
115+
@click.option('--land_mosaic', type=click.Path(exists=True))
116+
117+
@click.option('--wave_mosaic', type=click.Path(exists=True))
118+
119+
@click.option('--interp_order', default=2)
120+
121+
@click.option('--area_ratio_thresh', type=float, default=1e-6)
122+
123+
@click.option('--check')
124+
125+
@click.option('--rotate_poly')
126+
127+
def coupler(atmos_mosaic,
128+
ocean_mosaic,
129+
ocean_topog,
130+
mosaic_name,
131+
sea_level,
132+
land_mosaic,
133+
wave_mosaic,
134+
interp_order,
135+
area_ratio_thresh,
136+
check,
137+
rotate_poly):
138+
pass
139+
140+
#coupler.make(atmos_mosaic,
141+
# ocean_mosaic,
142+
# ocean_topog,
143+
# mosaic_name,
144+
# sea_level,
145+
# land_mosaic,
146+
# wave_mosaic,
147+
# interp_order,
148+
# area_ratio_thresh,
149+
# check
150+
# rotate_poly)
151+
152+
153+
154+
if __name__ == '__main__':
155+
make_mosaic()
Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
import sys
2+
import xarray as xr
3+
import numpy as np
4+
from FMSgridtools.shared.mosaicobj import MosaicObj
5+
6+
def make(global_mosaic,
7+
regional_file):
8+
""" Generates a horizontal grid and solo mosaic for a regional output.
9+
The created grid and solo mosaic could be used to regrid regional
10+
output data onto regular lat-lon grid."""
11+
12+
#get tile number from regional file
13+
try:
14+
tile = int(list(filter(str.isdigit, regional_file))[0])
15+
except IndexError:
16+
sys.exit("make_regional_mosaic: tile number not found in the regional_file name")
17+
18+
ds = xr.open_dataset(regional_file)
19+
nx = ds.sizes['grid_xt_sub01']
20+
ny = ds.sizes['grid_yt_sub01']
21+
22+
indx = ds.grid_xt_sub01.values
23+
indy = ds.grid_yt_sub01.values
24+
25+
i_min = np.min(indx)
26+
i_max = np.max(indx)
27+
j_min = np.min(indy)
28+
j_max = np.max(indy)
29+
30+
if i_max-i_min+1 != nx:
31+
print("Error: make_regional_mosaic: i_max-i_min+1 != nx")
32+
if j_max-j_min+1 != ny:
33+
print("Error: make_regional_mosaic: j_max-j_min+1 != ny")
34+
35+
grid = MosaicObj(global_mosaic).griddict()
36+
try:
37+
x,y = grid[f'tile{tile}'].x, grid[f'tile{tile}'].y
38+
except KeyError:
39+
sys.exit("make_regional_mosaic: gridtile with matching tile number not found")
40+
xarr = x[round(2*j_min - 2):round(2*j_min - 2 + 2*ny+1), round(2*i_min - 2):round(2*i_min - 2 + 2*nx+1)]
41+
yarr = y[round(2*j_min - 2):round(2*j_min - 2 + 2*ny+1), round(2*i_min - 2):round(2*i_min - 2 + 2*nx+1)]
42+
43+
outfile = f"regional_grid.tile{tile}.nc"
44+
write_out_regional_grid(tile, xarr,
45+
yarr, outfile)
46+
47+
mosaic = "regional_mosaic.nc"
48+
regional_mosaic = MosaicObj(mosaic_name="regional_mosaic",
49+
gridfiles=np.array([outfile]),
50+
gridtiles=np.array([f"tile{tile}"]))
51+
regional_mosaic.write_out_regional_mosaic(mosaic)
52+
53+
print("\nCongratulations: You have successfully run regional mosaic")
54+
55+
56+
def write_out_regional_grid(tile, xarr, yarr, outfile):
57+
58+
tile = xr.DataArray(
59+
data=f'tile{tile}',
60+
attrs=dict(standard_name="grid_tile_spec")).astype('|S255')
61+
62+
x = xr.DataArray(
63+
data=xarr,
64+
dims=["nyp", "nxp"],
65+
attrs=dict(
66+
standard_name="geographic_longitude", units="degree_east"))
67+
68+
y = xr.DataArray(
69+
data=yarr,
70+
dims=["nyp", "nxp"],
71+
attrs=dict(
72+
standard_name="geographic_latitude", units="degree_north"))
73+
74+
ds = xr.Dataset(
75+
data_vars={"tile": tile,
76+
"x": x,
77+
"y": y})
78+
79+
encoding = {'x': {'_FillValue': None},
80+
'y': {'_FillValue': None}}
81+
82+
ds.to_netcdf(outfile, encoding=encoding)
83+
84+
Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
import sys
2+
import numpy as np
3+
from FMSgridtools.shared.gridobj import GridObj
4+
from FMSgridtools.shared.mosaicobj import MosaicObj
5+
import pyfrenctools
6+
7+
def make(num_tiles,
8+
dir_name,
9+
mosaic_name,
10+
tile_file,
11+
periodx,
12+
periody) -> None:
13+
"""Generates mosaic information between tiles. The mosaic information includes:
14+
list of tile files, list of contact region,
15+
specified by index, contact type."""
16+
17+
contacts = []
18+
contact_index = []
19+
tilefiles = list(tile_file)
20+
nfiles = len(tilefiles)
21+
22+
if nfiles == 0:
23+
if num_tiles == 1:
24+
tilefiles.append("horizontal_grid.tile1.nc")
25+
else:
26+
for n in range(num_tiles):
27+
tilefiles.append(f"horizontal_grid.tile{n}.nc")
28+
else:
29+
if nfiles != num_tiles:
30+
sys.exit("Error: number tiles provided through --ntiles"
31+
"does not match number of grid files passed through")
32+
33+
gridtiles = [f'tile{i}' for i in range(1,nfiles+1)]
34+
35+
36+
grid = MosaicObj(ntiles=num_tiles, gridfiles = tilefiles,
37+
gridtiles = gridtiles).griddict()
38+
39+
grid_data = {}
40+
for tile in gridtiles:
41+
grid_data[tile] = {}
42+
grid_data[tile]['x'] = grid[tile].x
43+
grid_data[tile]['y'] = grid[tile].y
44+
grid_data[tile]['nxp'] = grid[tile].nxp
45+
grid_data[tile]['nyp'] = grid[tile].nyp
46+
47+
ncontact = 0
48+
#FIND CONTACT REGIONS
49+
for n in range(1,num_tiles):
50+
for m in range(n,num_tiles+1):
51+
contact = Contact(n, m, grid_data[f'tile{n}']['nxp'],
52+
grid_data[f'tile{m}']['nxp'],
53+
grid_data[f'tile{n}']['nyp'],
54+
grid_data[f'tile{m}']['nyp'],
55+
grid_data[f'tile{n}']['x'],
56+
grid_data[f'tile{m}']['x'],
57+
grid_data[f'tile{n}']['y'],
58+
grid_data[f'tile{m}']['y'],
59+
periodx, periody)
60+
count, istart1, iend1, jstart1, jend1, istart2, iend2, jstart2, jend2 = contact.align_contact()
61+
62+
if count > 0:
63+
contacts.append(f"{mosaic_name}:tile{n}::{mosaic_name}:tile{m}")
64+
ncontact+=1
65+
66+
67+
contact_index.append(f"{istart1}:{iend1},{jstart1}:{jend1}::{istart2}:{iend2},{jstart2}:{jend2}")
68+
69+
print("\nCongratulations: You have successfully run solo mosaic")
70+
print(f"NOTE: There are {ncontact} contacts\n")
71+
72+
if ncontact > 0:
73+
mosaic = MosaicObj(mosaic_name=mosaic_name,
74+
gridlocation=dir_name,
75+
gridfiles=tilefiles,
76+
gridtiles=gridtiles,
77+
contacts=contacts,
78+
contact_index=contact_index)
79+
mosaic.write_out_mosaic(f'{mosaic_name}.nc')

FMSgridtools/shared/mosaicobj.py

100755100644
File mode changed.

0 commit comments

Comments
 (0)