Skip to content

Commit

Permalink
First attempt to resolve issue #22
Browse files Browse the repository at this point in the history
 - add tests to ensure issue #22 stays resolved
 - when no land points are in the ocean grid, skip writing atmos/land exchange file
jr3cermak committed May 12, 2022

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
1 parent 80f8ddc commit 81d771a
Showing 7 changed files with 503 additions and 7 deletions.
5 changes: 4 additions & 1 deletion docs/development/CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -4,7 +4,10 @@

Updates

- Issue #22: fixing up creation of mosaics and exchange grids depending on
- Resolve #21: makeSoloMosaic() will create the appropriate exchange grids
depending if land points exist within an ocean grid.
- Tests are now included to make sure issue #22 stays operational.
- Issue #21: fixing up creation of mosaics and exchange grids depending on
what is seen in the land mask. Fix for ocean grid `tile` variable to
make sure it is of type `NC_CHAR`.

10 changes: 5 additions & 5 deletions docs/development/TODO.md
Original file line number Diff line number Diff line change
@@ -3,12 +3,12 @@
## Milestones

- [ ] Release 0.3.3
- [ ] Resolve Issue #21
- [X] Resolve Issue #21
- [X] Save `tile` variable as type NC_CHAR in ocean grid
- [ ] Add test for regular grid with land
- [ ] Add test for regular grid without land
- [ ] Add test for regular grid specified to be without land
- [ ] Add test for regular grid with land but error out with noLandMosaic=True specified
- [X] Add test for regular grid with land
- [X] Add test for regular grid without land
- [X] Add test for regular grid specified to be without land
- [X] Add test for regular grid with land and noLandMosaic=True specified (expected failure)
- [ ] Resolve Issue #19
- [X] add a keyword option to trigger alternate calcuation of angle_dx
in gridutils.grids.mom6 and gridutils.
10 changes: 9 additions & 1 deletion gridtools/grids/mom6.py
Original file line number Diff line number Diff line change
@@ -561,6 +561,11 @@ def write_MOM6_exchange_grid_files(self, grd, **kwargs):
if name1 == 'land' or name2 == 'land':
continue

# If there are no land grid points, skip atmos/land exchange file
if not(kwargs['haveLandPoints']):
if name1 == 'atmos' and name2 == 'land':
continue

ds = xr.Dataset()

# calculate the exchange grid for name1 X name2
@@ -742,9 +747,12 @@ def add_string_var_2d(ds, var_name, dim_name, standard_name, value, strVarMap):

add_string_var_2d(ds, 'aXo_file', 'nfile_aXo', 'atmXocn_exchange_grid_file', self.mom6_grid['filenames']['atmos_ocean_exchange'], strVarMap)
if not(kwargs['noLandMosaic']):
add_string_var_2d(ds, 'aXl_file', 'nfile_aXl', 'atmXlnd_exchange_grid_file', self.mom6_grid['filenames']['atmos_land_exchange'] , strVarMap)
# If no land points exist in the ocean grid, do not use atmos to land exchange
if kwargs['haveLandPoints']:
add_string_var_2d(ds, 'aXl_file', 'nfile_aXl', 'atmXlnd_exchange_grid_file', self.mom6_grid['filenames']['atmos_land_exchange'] , strVarMap)
add_string_var_2d(ds, 'lXo_file', 'nfile_lXo', 'lndXocn_exchange_grid_file', self.mom6_grid['filenames']['land_ocean_exchange'] , strVarMap)
else:
# If no land mosaic is specified, the atmos mosaic is used instead
add_string_var_2d(ds, 'lXo_file', 'nfile_lXo', 'lndXocn_exchange_grid_file', self.mom6_grid['filenames']['atmos_ocean_exchange'], strVarMap)

# Global attributes
124 changes: 124 additions & 0 deletions tests/mosaics/test_smallgrid_regular_land_point.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
# TEST: Generate a regular grid with one land point.
# Expected result: Three exchange grids should be generated.

import pytest
xfail = pytest.mark.xfail

from gridtools.gridutils import GridUtils
import os, glob
import numpy as np
import xarray as xr

def clean_output_directory(tmp_path, glob_patterns):

# Make sure no matching glob patterns are found
# before running this test

for p in glob_patterns:
mfiles = glob.glob(p)
for mfile in mfiles:
if os.path.isfile(mfile):
os.unlink(mfile)

def test_regular_grid_with_land_point(tmp_path):

# Create a gridtools object
grd = GridUtils()

# Define a small 5x5 grid with 1000.0 meter resolution
grd.setGridParameters({
'projection': {
'name': "Stereographic",
'ellps': 'WGS84',
'lon_0': 0.0,
'lat_0': 90.0,
'lat_ts': 75.0,
},
'centerX': 5.0,
'centerY': 75.0,
'cneterUnits': 'degrees',
'dx': 5000.0,
'dy': 5000.0,
'dxUnits': 'meters',
'dyUnits': 'meters',
'gridResolution': 1000.0,
'gridResolutionUnits': 'meters',
'tilt': 0.0,
'gridMode': 2,
'gridType': 'MOM6',
'ensureEvenI': True,
'ensureEvenJ': True
})

grd.makeGrid()
wrkDir = tmp_path
clean_output_directory(tmp_path, ["*.nc"])
grd.saveGrid(filename=os.path.join(wrkDir, "ocean_hgrid.nc"))

h = 30. * np.ones((len(grd.grid.ny)//2, len(grd.grid.nx)//2))
# Create a land point
h[1,1] = 0.0
m = np.ones((len(grd.grid.ny)//2, len(grd.grid.nx)//2))
bathyGrids = xr.Dataset({'depth': (['ny', 'nx'], h),
'mask': (['ny', 'nx'], m) })
grd.writeOceanmask(bathyGrids, 'depth', 'mask',
os.path.join(wrkDir, 'ocean_mask.nc'),
MASKING_DEPTH=0.0)
grd.writeLandmask(bathyGrids, 'depth', 'mask',
os.path.join(wrkDir, 'land_mask.nc'),
MASKING_DEPTH=0.0)
bathyGrids.to_netcdf(os.path.join(wrkDir, 'ocean_topog.nc'),
encoding=grd.removeFillValueAttributes(data=bathyGrids))
grd.makeSoloMosaic(
topographyGrid=bathyGrids['depth'],
writeLandmask=True,
writeOceanmask=True,
inputDirectory=wrkDir,
overwrite=True,
noLandMosaic=False
)

# Examine output

# Test fails if one of these files is missing
files_present = ["atmos_mosaic_tile1Xland_mosaic_tile1.nc",\
"atmos_mosaic_tile1Xocean_mosaic_tile1.nc",\
"land_mosaic_tile1Xocean_mosaic_tile1.nc"]

for file_to_check in files_present:
full_path_file = os.path.join(tmp_path, file_to_check)
if not(os.path.isfile(full_path_file)):
assert False

# Do some variable size and content checks
file_to_check = "atmos_mosaic_tile1Xland_mosaic_tile1.nc"
ds = xr.open_dataset(os.path.join(tmp_path,file_to_check))
mosaic_length = (ds['tile1_cell'].shape)[0]
if mosaic_length != 1:
print("%s should have a length of 1 (%d)." % (file_to_check,\
mosaic_length))
assert False
ds.close()

file_to_check = "land_mosaic_tile1Xocean_mosaic_tile1.nc"
ds = xr.open_dataset(os.path.join(tmp_path,file_to_check))
mosaic_length = (ds['tile1_cell'].shape)[0]
if mosaic_length != 24:
print("%s should have a length of 24 (%d)." % (file_to_check,\
mosaic_length))
assert False
ds.close()

file_to_check = "ocean_hgrid.nc"
ds = xr.open_dataset(os.path.join(tmp_path,file_to_check))
print(ds)
if ds['tile'].dtype.kind != 'S':
print("Variable 'tile' in %s should be of kind (S) string." %\
(file_to_check))
assert False
if ds['tile'].dtype.itemsize != 5:
print("Variable 'tile' in %s should have an itemsize of 5." %\
(file_to_check))
assert False

assert True
127 changes: 127 additions & 0 deletions tests/mosaics/test_smallgrid_regular_land_point_noLandMosaic_True.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,127 @@
# TEST: Generate a regular grid with one land point but with the attempt of
# using the noLandMosaic=True option.
# Expected result: An ERROR should be generated as there is a land point
# in the grid but asked to not generate a land mosaic file. A land
# mosaic file should exist.

import pytest
xfail = pytest.mark.xfail

from gridtools.gridutils import GridUtils
import os, glob
import numpy as np
import xarray as xr

def clean_output_directory(tmp_path, glob_patterns):

# Make sure no matching glob patterns are found
# before running this test

for p in glob_patterns:
mfiles = glob.glob(p)
for mfile in mfiles:
if os.path.isfile(mfile):
os.unlink(mfile)

@xfail
def test_regular_grid_with_land_point(tmp_path):

# Create a gridtools object
grd = GridUtils()

# Define a small 5x5 grid with 1000.0 meter resolution
grd.setGridParameters({
'projection': {
'name': "Stereographic",
'ellps': 'WGS84',
'lon_0': 0.0,
'lat_0': 90.0,
'lat_ts': 75.0,
},
'centerX': 5.0,
'centerY': 75.0,
'cneterUnits': 'degrees',
'dx': 5000.0,
'dy': 5000.0,
'dxUnits': 'meters',
'dyUnits': 'meters',
'gridResolution': 1000.0,
'gridResolutionUnits': 'meters',
'tilt': 0.0,
'gridMode': 2,
'gridType': 'MOM6',
'ensureEvenI': True,
'ensureEvenJ': True
})

grd.makeGrid()
wrkDir = tmp_path
clean_output_directory(tmp_path, ["*.nc"])
grd.saveGrid(filename=os.path.join(wrkDir, "ocean_hgrid.nc"))

h = 30. * np.ones((len(grd.grid.ny)//2, len(grd.grid.nx)//2))
# Create a land point
h[1,1] = 0.0
m = np.ones((len(grd.grid.ny)//2, len(grd.grid.nx)//2))
bathyGrids = xr.Dataset({'depth': (['ny', 'nx'], h),
'mask': (['ny', 'nx'], m) })
grd.writeOceanmask(bathyGrids, 'depth', 'mask',
os.path.join(wrkDir, 'ocean_mask.nc'),
MASKING_DEPTH=0.0)
grd.writeLandmask(bathyGrids, 'depth', 'mask',
os.path.join(wrkDir, 'land_mask.nc'),
MASKING_DEPTH=0.0)
bathyGrids.to_netcdf(os.path.join(wrkDir, 'ocean_topog.nc'),
encoding=grd.removeFillValueAttributes(data=bathyGrids))
grd.makeSoloMosaic(
topographyGrid=bathyGrids['depth'],
writeLandmask=True,
writeOceanmask=True,
inputDirectory=wrkDir,
overwrite=True,
noLandMosaic=True
)

# Examine output

# Test fails if one of these files is missing
files_present = ["atmos_mosaic_tile1Xland_mosaic_tile1.nc",\
"atmos_mosaic_tile1Xocean_mosaic_tile1.nc",\
"land_mosaic_tile1Xocean_mosaic_tile1.nc"]

for file_to_check in files_present:
full_path_file = os.path.join(tmp_path, file_to_check)
if not(os.path.isfile(full_path_file)):
assert False

# Do some variable size and content checks
file_to_check = "atmos_mosaic_tile1Xland_mosaic_tile1.nc"
ds = xr.open_dataset(os.path.join(tmp_path,file_to_check))
mosaic_length = (ds['tile1_cell'].shape)[0]
if mosaic_length != 1:
print("%s should have a length of 1 (%d)." % (file_to_check,\
mosaic_length))
assert False
ds.close()

file_to_check = "land_mosaic_tile1Xocean_mosaic_tile1.nc"
ds = xr.open_dataset(os.path.join(tmp_path,file_to_check))
mosaic_length = (ds['tile1_cell'].shape)[0]
if mosaic_length != 24:
print("%s should have a length of 24 (%d)." % (file_to_check,\
mosaic_length))
assert False
ds.close()

file_to_check = "ocean_hgrid.nc"
ds = xr.open_dataset(os.path.join(tmp_path,file_to_check))
if ds['tile'].dtype.kind != 'S':
print("Variable 'tile' in %s should be of kind (S) string." %\
(file_to_check))
assert False
if ds['tile'].dtype.itemsize != 5:
print("Variable 'tile' in %s should have an itemsize of 5." %\
(file_to_check))
assert False

assert True
Loading

0 comments on commit 81d771a

Please sign in to comment.