Skip to content

Commit 9943058

Browse files
mgdudamgduda
committed
Merge branch 'hexagon_projections' (PR #606)
This merge adds a tool for creating a limited-area hexagonal mesh on a Lambert conformal projection. The tool works by projecting a rectangular mesh composed of perfect hexagons to the sphere using a Lambert conformal projection. The hex_projection directory in the mesh_tools/ directory is new and the utility is placed in hex_projection directory. The README file provides directions for building and running the utility, and a default case is specified in the namelist.projections file used as input to the utility. The utility creates an MPAS grid file (netcdf) and a graph partition file (ascii) for the mesh. The grid file includes the masks needed for regional simulations in MPAS-Atmosphere, along with a mesh_density function. The current mask configuration has 5 levels in the relaxation zone and 2 levels in the specified zone. This is hardwired in the utility, consistent with the hardwired configuration in MPAS-Atmosphere. * hexagon_projections: This commit adds the capability of creating a regional MPAS horizontal mesh by...
2 parents 00920f0 + db60d93 commit 9943058

File tree

12 files changed

+3968
-0
lines changed

12 files changed

+3968
-0
lines changed

‎mesh_tools/hex_projection/Makefile‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
FC = $(shell nf-config --fc)
2+
FFLAGS = -O3
3+
FCINCLUDES = $(shell nf-config --fflags)
4+
RPATH_FLAGS = $(shell nf-config --flibs | grep -o -e '-L\S\+\( \|$$\)' | sed 's/^-L/-Wl,-rpath,/' | tr -d '\n')
5+
RPATH_FLAGS += $(shell nc-config --libs | grep -o -e '-L\S\+\( \|$$\)' | sed 's/^-L/-Wl,-rpath,/' | tr -d '\n')
6+
FCLIBS = -L$(shell nc-config --libdir) $(shell nf-config --flibs) $(RPATH_FLAGS)
7+
8+
OBJS = project_hexes.o \
9+
mpas_geometry_utilities.o \
10+
ph_utils.o \
11+
ph_write_mesh.o \
12+
projection_setup.o \
13+
projections.o \
14+
precision.o
15+
16+
EXE = project_hexes
17+
18+
$(EXE) : $(OBJS)
19+
$(FC) -o $@ $(OBJS) $(FCLIBS)
20+
21+
project_hexes.o : mpas_geometry_utilities.o ph_utils.o ph_write_mesh.o projection_setup.o projections.o precision.o
22+
23+
mpas_geometry_utilities.o : precision.o
24+
25+
ph_utils.o: mpas_geometry_utilities.o precision.o
26+
27+
ph_write_mesh.o : precision.o
28+
29+
projection_setup.o : precision.o
30+
31+
projections.o : projection_setup.o precision.o
32+
33+
clean:
34+
$(RM) $(OBJS) *.mod $(EXE)
35+
36+
37+
%.o : %.mod
38+
39+
%.o : %.f90
40+
$(RM) $@ $*.mod
41+
$(FC) $(FFLAGS) -c $< $(FCINCLUDES)

‎mesh_tools/hex_projection/README‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
MPAS-Atmosphere utility for creating regional meshes by projecting a
2+
Cartesian mesh of perfect hexagons to the sphere using standard map
3+
projections
4+
5+
The current implementation uses a rectangular Cartesian mesh, with the
6+
corners slightly rounded, and projects this mesh to the sphere using a
7+
Lambert Conformal projection. Executing "make" with the included
8+
Makefile produces an executable named "project_hexes". The inputs
9+
to this executable are given in the ascii file "namelist.projections"
10+
(a standard Fortran namelist file). The inputs variables should be
11+
obvious from their names, and some further information is given in
12+
"Notes" below. The utility produces a netcdf file "mpas_hex_mesh.nc"
13+
describing the regional mesh in the form of a standard MPAS grid file,
14+
with the projection existing on the unit sphere (as is the case with
15+
all MPAS grid files for spheres). The utility also produces a
16+
graph.info file that serves as input to METIS to produce a partion
17+
file for use by MPAS as described in the MPAS-A Users Guide.
18+
19+
Notes:
20+
21+
The (x,y) extent of the domain is in meters on the earth-radius sphere,
22+
but as noted the MPAS grid file that is produced is on the unit sphere.
23+
24+
Following WRF, the longitudinal center of the project is given by the
25+
namelist parameter "standard_longitude_degrees". The center of the
26+
projected mesh is given by the namelist parameter
27+
"reference_longitude_degrees".
28+
29+
The hexagons are oriented pointing west-east.
30+
31+
A python script is included that can plot the resulting mesh.
32+
An NCL script (a depricated language) is also provided to plot
33+
the mesh.
34+
35+
The mesh can be moved and rotated to any place on the sphere using
36+
the MPAS grid_rotate tool (see the MPAS Users Guide).
37+
38+
WCS, 4 March 2025

‎mesh_tools/hex_projection/mesh_boundary_plot.py‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
import cartopy.crs as ccrs
2+
import matplotlib.pyplot as plt
3+
import numpy as np
4+
import math
5+
try:
6+
from numba import njit
7+
except ImportError:
8+
njit = lambda f: f
9+
10+
@njit
11+
def boundary_vertices(lonVertex, latVertex, bdyMaskVertex, edgesOnVertex, verticesOnEdge, bdyMask):
12+
bdyLats = []
13+
bdyLons = []
14+
15+
for startVertex in range(bdyMaskVertex.size):
16+
if bdyMaskVertex[startVertex] == bdyMask:
17+
break
18+
19+
for edge in edgesOnVertex[startVertex][:]:
20+
if edge == -1:
21+
continue
22+
23+
if verticesOnEdge[edge][0] != startVertex:
24+
neighbor = verticesOnEdge[edge][0]
25+
else:
26+
neighbor = verticesOnEdge[edge][1]
27+
28+
if bdyMaskVertex[neighbor] == bdyMask:
29+
nextVertex = neighbor
30+
bdyLons.append(lonVertex[nextVertex])
31+
bdyLats.append(latVertex[nextVertex])
32+
break
33+
34+
prevVertex = startVertex
35+
36+
while nextVertex != startVertex:
37+
for edge in edgesOnVertex[nextVertex][:]:
38+
if edge == -1:
39+
continue
40+
41+
if verticesOnEdge[edge][0] != nextVertex:
42+
neighbor = verticesOnEdge[edge][0]
43+
else:
44+
neighbor = verticesOnEdge[edge][1]
45+
46+
if bdyMaskVertex[neighbor] == bdyMask and neighbor != prevVertex:
47+
prevVertex = nextVertex
48+
nextVertex = neighbor
49+
bdyLons.append(lonVertex[nextVertex])
50+
bdyLats.append(latVertex[nextVertex])
51+
break
52+
53+
return np.asarray(bdyLons), np.asarray(bdyLats)
54+
55+
56+
def mesh_extent(latField, lonField):
57+
earthRadius = 6371229.0
58+
extentLat = (max(latField) - min(latField)) * earthRadius * math.pi / 180.0
59+
extentLon = (max(lonField) - min(lonField)) * math.cos(cenLat * math.pi / 180.0) * earthRadius * math.pi / 180.0
60+
return max(extentLat, extentLon)
61+
62+
63+
def map_background(cenLat, cenLon, extent):
64+
import cartopy.feature as cfeature
65+
import matplotlib.ticker as mticker
66+
67+
proj = ccrs.Stereographic(cenLat, cenLon)
68+
ax = plt.axes(projection=proj)
69+
70+
#### MGD
71+
# print(proj.transform_point(-110.0, 40.0, ccrs.PlateCarree()))
72+
73+
ax.set_title('Limited-area mesh')
74+
75+
scaling = 1.5
76+
ax.set_extent([-scaling * extent, scaling * extent, -scaling * extent, scaling * extent], crs=proj)
77+
78+
ax.add_feature(cfeature.LAND)
79+
ax.add_feature(cfeature.OCEAN)
80+
ax.add_feature(cfeature.COASTLINE, linewidth=0.1)
81+
ax.add_feature(cfeature.BORDERS, linewidth=0.1)
82+
ax.add_feature(cfeature.LAKES, linewidth=0.1)
83+
ax.add_feature(cfeature.RIVERS, linewidth=0.1)
84+
ax.add_feature(cfeature.STATES, linewidth=0.1)
85+
86+
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
87+
linewidth=0.1, color='black', alpha=1.0, linestyle='--')
88+
89+
xticks = np.arange(-180, 180, 15)
90+
yticks = np.arange(-90, 90, 15)
91+
92+
gl.ylocator = mticker.FixedLocator(yticks)
93+
gl.xlocator = mticker.FixedLocator(xticks)
94+
95+
return proj, ax
96+
97+
98+
def draw_domain(ax, bdyLats, bdyLons, bdyLatsSpec, bdyLonsSpec):
99+
import matplotlib.patches as mpatches
100+
101+
poly_corners = np.zeros((len(bdyLons), 2), np.float64)
102+
poly_corners[:,0] = np.asarray(bdyLons)
103+
poly_corners[:,1] = np.asarray(bdyLats)
104+
105+
poly = mpatches.Polygon(poly_corners, closed=True, ec='black', fill=True, lw=0.1, fc='black', alpha=0.3, transform=ccrs.Geodetic())
106+
ax.add_patch(poly)
107+
108+
poly_corners = np.zeros((len(bdyLonsSpec), 2), np.float64)
109+
poly_corners[:,0] = np.asarray(bdyLonsSpec)
110+
poly_corners[:,1] = np.asarray(bdyLatsSpec)
111+
112+
poly = mpatches.Polygon(poly_corners, closed=True, ec='black', fill=True, lw=0.1, fc='black', alpha=0.3, transform=ccrs.Geodetic())
113+
ax.add_patch(poly)
114+
115+
116+
if __name__ == '__main__':
117+
from time import time
118+
from netCDF4 import Dataset
119+
120+
t1 = time()
121+
f = Dataset('mpas_hex_mesh.nc')
122+
123+
latVertex = np.ma.getdata(f.variables['latVertex'][:]) * 180.0 / math.pi
124+
lonVertex = np.ma.getdata(f.variables['lonVertex'][:]) * 180.0 / math.pi
125+
bdyMaskVertex = np.ma.getdata(f.variables['bdyMaskVertex'][:])
126+
edgesOnVertex = np.ma.getdata(f.variables['edgesOnVertex'][:]) - 1
127+
verticesOnEdge = np.ma.getdata(f.variables['verticesOnEdge'][:]) - 1
128+
129+
f.close()
130+
t2 = time()
131+
print('Time to reading input fields: ', (t2 - t1))
132+
133+
t1 = time()
134+
cenLat = np.sum(latVertex) / latVertex.size
135+
cenLon = np.sum(lonVertex) / lonVertex.size
136+
extent = mesh_extent(latVertex, lonVertex)
137+
t2 = time()
138+
print('Time to calculate domain extents ', (t2 - t1))
139+
140+
t1 = time()
141+
proj, ax = map_background(cenLat, cenLon, extent)
142+
t2 = time()
143+
print('Time to set up map background: ', (t2 - t1))
144+
145+
t1 = time()
146+
bdyLons, bdyLats = boundary_vertices(lonVertex, latVertex, bdyMaskVertex, edgesOnVertex, verticesOnEdge, 1)
147+
bdyLonsSpec, bdyLatsSpec = boundary_vertices(lonVertex, latVertex, bdyMaskVertex, edgesOnVertex, verticesOnEdge, 7)
148+
t2 = time()
149+
print('Time to find domain boundaries: ', (t2 - t1))
150+
151+
t1 = time()
152+
draw_domain(ax, bdyLats, bdyLons, bdyLatsSpec, bdyLonsSpec)
153+
t2 = time()
154+
print('Time to draw boundaries: ', (t2 - t1))
155+
156+
t1 = time()
157+
plt.savefig('mesh.pdf')
158+
t2 = time()
159+
print('Time to save figure: ', (t2 - t1))

‎mesh_tools/hex_projection/mesh_plot.ncl‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
2+
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
3+
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
4+
5+
begin
6+
7+
;******************************************************************************
8+
9+
wks = gsn_open_wks("pdf","mesh_ncl_plot")
10+
11+
f_in = addfile("mpas_hex_mesh.nc","r")
12+
13+
;******************************************************************************
14+
15+
res = True
16+
; res@gsnMaximize = True
17+
res@gsnDraw = False
18+
res@gsnFrame = False
19+
res@gsnSpreadColors = True
20+
21+
; res@mpProjection = "CylindricalEquidistant"
22+
res@mpProjection = "Orthographic"
23+
res@mpDataBaseVersion = "MediumRes"
24+
res@mpCenterLatF = 38.5
25+
; res@mpCenterLatF = 0.0
26+
res@mpCenterLonF = -97.5
27+
28+
res@mpMinLonF = -157.5
29+
res@mpMaxLonF = -37.5
30+
res@mpMinLatF = 10.
31+
res@mpMaxLatF = 70.
32+
33+
res@cnFillOn = True
34+
res@cnFillMode = "AreaFill"
35+
res@cnLinesOn = False
36+
;res@cnLineLabelsOn = True
37+
;res@cnLineThicknessF = 2.
38+
39+
res@cnLevelSelectionMode = "ManualLevels"
40+
41+
res@tiMainFontHeightF = 0.016
42+
res@tiMainString = "MPAS regional mesh"
43+
res@gsnLeftStringFontHeightF = 0.015
44+
res@gsnRightStringFontHeightF = 0.015
45+
46+
;res@lbLabelStride = 2
47+
48+
;-----1. contour for global z500 -------------------------------------------
49+
50+
; gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
51+
; gsn_define_colormap(wks,"BlAqGrYeOrReVi200")
52+
gsn_define_colormap(wks,"perc2_9lev")
53+
54+
;res@gsnSpreadColors = False
55+
;res@gsnSpreadColorStart = 2
56+
;res@gsnSpreadColorEnd = 201
57+
58+
r2d = 180.D/(4.*atan(1.))
59+
lat1d = f_in->latCell(:)*r2d
60+
lon1d = f_in->lonCell(:)*r2d
61+
62+
res@sfXArray = lon1d
63+
res@sfYArray = lat1d
64+
65+
66+
ibdymask = f_in->bdyMaskCell
67+
bdymask = int2dble(ibdymask)
68+
nCells = dimsizes(lat1d)
69+
70+
res@cnMinLevelValF = 0.5
71+
res@cnMaxLevelValF = 7.5
72+
res@cnLevelSpacingF = 1.0
73+
74+
75+
; plot = gsn_csm_contour_map_ce(wks,bdymask,res)
76+
plot = gsn_csm_contour_map(wks,bdymask,res)
77+
draw(plot)
78+
frame(wks)
79+
80+
81+
end
82+

0 commit comments

Comments
 (0)