From c4c56d2895417ea16d19161c5267c90b0b3def10 Mon Sep 17 00:00:00 2001 From: Thomas Saillour Date: Mon, 4 Dec 2023 19:45:11 +0100 Subject: [PATCH] Added support for global meshes (stereographic projection) (#68) --- README.md | 87 +++++++++++++++++--- oceanmesh/__init__.py | 13 ++- oceanmesh/boundary.py | 2 +- oceanmesh/clean.py | 2 +- oceanmesh/edgefx.py | 61 ++++++++++++-- oceanmesh/filterfx.py | 6 +- oceanmesh/fix_mesh.py | 3 +- oceanmesh/geodata.py | 36 ++++++-- oceanmesh/grid.py | 46 +++++++++-- oceanmesh/mesh_generator.py | 63 +++++++++++--- oceanmesh/region.py | 57 +++++++++++++ tests/global/global_latlon.cpg | 1 + tests/global/global_latlon.dbf | Bin 0 -> 1225 bytes tests/global/global_latlon.prj | 1 + tests/global/global_latlon.shp | Bin 0 -> 40236 bytes tests/global/global_latlon.shx | Bin 0 -> 588 bytes tests/global/global_stereo.cpg | 1 + tests/global/global_stereo.dbf | Bin 0 -> 8622 bytes tests/global/global_stereo.prj | 1 + tests/global/global_stereo.shp | Bin 0 -> 40116 bytes tests/global/global_stereo.shx | Bin 0 -> 580 bytes tests/test_README.py | 8 -- tests/test_bathymetric_gradient_function.py | 4 +- tests/test_edgefx.py | 10 +-- tests/test_geodata.py | 2 +- tests/test_global_stereo.py | 60 ++++++++++++++ tests/test_grade.py | 2 +- tox.ini | 2 +- 28 files changed, 398 insertions(+), 70 deletions(-) create mode 100644 tests/global/global_latlon.cpg create mode 100644 tests/global/global_latlon.dbf create mode 100644 tests/global/global_latlon.prj create mode 100644 tests/global/global_latlon.shp create mode 100644 tests/global/global_latlon.shx create mode 100644 tests/global/global_stereo.cpg create mode 100644 tests/global/global_stereo.dbf create mode 100644 tests/global/global_stereo.prj create mode 100644 tests/global/global_stereo.shp create mode 100644 tests/global/global_stereo.shx delete mode 100644 tests/test_README.py create mode 100644 tests/test_global_stereo.py diff --git a/README.md b/README.md index 1882184..8451919 100644 --- a/README.md +++ b/README.md @@ -275,12 +275,12 @@ extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG) min_edge_length = 0.01 # minimum mesh size in domain in projection shoreline = om.Shoreline(fname, extent.bbox, min_edge_length) edge_length = om.distance_sizing_function(shoreline, rate=0.15) -ax = edge_length.plot( +fig, ax, pc = edge_length.plot( xlabel="longitude (WGS84 degrees)", ylabel="latitude (WGS84 degrees)", title="Distance sizing function", cbarlabel="mesh size (degrees)", - hold=True, + holding=True, ) shoreline.plot(ax=ax) ``` @@ -303,12 +303,12 @@ sdf = om.signed_distance_function(shoreline) edge_length = om.feature_sizing_function( shoreline, sdf, max_edge_length=0.05, plot=True ) -ax = edge_length.plot( +fig, ax, pc = edge_length.plot( xlabel="longitude (WGS84 degrees)", ylabel="latitude (WGS84 degrees)", title="Feature sizing function", cbarlabel="mesh size (degrees)", - hold=True, + holding=True, xlim=[-74.3, -73.8], ylim=[40.3, 40.8], ) @@ -332,12 +332,12 @@ shoreline = om.Shoreline(fname, extent.bbox, min_edge_length) sdf = om.signed_distance_function(shoreline) edge_length = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05) edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.15) -ax = edge_length.plot( +fig, ax, pc = edge_length.plot( xlabel="longitude (WGS84 degrees)", ylabel="latitude (WGS84 degrees)", title="Feature sizing function with gradation bound", cbarlabel="mesh size (degrees)", - hold=True, + holding=True, xlim=[-74.3, -73.8], ylim=[40.3, 40.8], ) @@ -358,7 +358,8 @@ fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp" min_edge_length = 0.01 -dem = om.DEM(fdem, crs=4326) +extent = om.Region(extent=(-74.3, -73.8, 40.3,40.8), crs=4326) +dem = om.DEM(fdem, bbox=extent,crs=4326) shoreline = om.Shoreline(fname, dem.bbox, min_edge_length) sdf = om.signed_distance_function(shoreline) edge_length1 = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05) @@ -368,12 +369,12 @@ edge_length2 = om.wavelength_sizing_function( # Compute the minimum of the sizing functions edge_length = om.compute_minimum([edge_length1, edge_length2]) edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.15) -ax = edge_length.plot( +fig, ax, pc = edge_length.plot( xlabel="longitude (WGS84 degrees)", ylabel="latitude (WGS84 degrees)", title="Feature sizing function + wavelength + gradation bound", cbarlabel="mesh size (degrees)", - hold=True, + holding=True, xlim=[-74.3, -73.8], ylim=[40.3, 40.8], ) @@ -385,6 +386,7 @@ shoreline.plot(ax=ax) The distance, feature size, and/or wavelength mesh size functions can lead to coarse mesh resolution in deeper waters that under-resolve and smooth over the sharp topographic gradients that characterize the continental shelf break. These slope features can be important for coastal ocean models in order to capture dissipative effects driven by the internal tides, transmissional reflection at the shelf break that controls the astronomical tides, and trapped shelf waves. The bathymetry field contains often excessive details that are not relevant for most flows, thus the bathymetry can be smoothed by a variety of filters (e.g., lowpass, bandpass, and highpass filters) before calculating the mesh sizes. + ```python import oceanmesh as om @@ -394,7 +396,7 @@ fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp" EPSG = 4326 # EPSG:4326 or WGS84 bbox = (-74.4, -73.4, 40.2, 41.2) extent = om.Region(extent=bbox, crs=EPSG) -dem = om.DEM(fdem, crs=4326) +dem = om.DEM(fdem, crs=EPSG) min_edge_length = 0.0025 # minimum mesh size in domain in projection max_edge_length = 0.10 # maximum mesh size in domain in projection @@ -414,7 +416,7 @@ edge_length2 = om.bathymetric_gradient_sizing_function( min_edge_length=min_edge_length, max_edge_length=max_edge_length, crs=EPSG, -) +) # will be reactivated edge_length3 = om.compute_minimum([edge_length1, edge_length2]) edge_length3 = om.enforce_mesh_gradation(edge_length3, gradation=0.15) ``` @@ -591,6 +593,69 @@ plt.show() ``` ![Multiscale](https://user-images.githubusercontent.com/18619644/136119785-8746552d-4ff6-44c3-9aa1-3e4981ba3518.png) +Global mesh generation +---------------------- +Using oceanmesh is now possible for global meshes. +The process is done in two steps: + * first the definition of the sizing functions in WGS84 coordinates, + * then the mesh generation is done in the stereographic projection + +```python +import os +import numpy as np +import oceanmesh as om +from oceanmesh.region import to_lat_lon +import matplotlib.pyplot as plt +# utilities functions for plotting +def crosses_dateline(lon1, lon2): + return abs(lon1 - lon2) > 180 + +def filter_triangles(points, cells): + filtered_cells = [] + for cell in cells: + p1, p2, p3 = points[cell[0]], points[cell[1]], points[cell[2]] + if not (crosses_dateline(p1[0], p2[0]) or crosses_dateline(p2[0], p3[0]) or crosses_dateline(p3[0], p1[0])): + filtered_cells.append(cell) + return filtered_cells + +# Note: global_stereo.shp has been generated using global_tag() function in pyposeidon +# https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452 +fname = "tests/global/global_latlon.shp" +fname2 = "tests/global/global_stereo.shp" + +EPSG = 4326 # EPSG:4326 or WGS84 +bbox = (-180.00, 180.00, -89.00, 90.00) +extent = om.Region(extent=bbox, crs=4326) + +min_edge_length = 0.5 # minimum mesh size in domain in meters +max_edge_length = 2 # maximum mesh size in domain in meters +shoreline = om.Shoreline(fname, extent.bbox, min_edge_length) +sdf = om.signed_distance_function(shoreline) +edge_length = om.distance_sizing_function(shoreline, rate=0.11) + +# once the size functions have been defined, wed need to mesh inside domain in +# stereographic projections. This is way we use another coastline which is +# already translated in a sterographic projection +shoreline_stereo = om.Shoreline(fname2, extent.bbox, min_edge_length, stereo=True) +domain = om.signed_distance_function(shoreline_stereo) + +points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=100) + +# remove degenerate mesh faces and other common problems in the mesh +points, cells = om.make_mesh_boundaries_traversable(points, cells) +points, cells = om.delete_faces_connected_to_one_face(points, cells) + +# apply a Laplacian smoother +points, cells = om.laplacian2(points, cells, max_iter=100) +lon, lat = to_lat_lon(points[:, 0], points[:, 1]) +trin = filter_triangles(np.array([lon,lat]).T, cells) + +fig, ax, pc = edge_length.plot(holding = True, plot_colorbar=True, cbarlabel = "Resolution in °", cmap='magma') +ax.triplot(lon, lat, trin, color="w", linewidth=0.25) +plt.tight_layout() +plt.show() +``` +![Global](https://github.com/tomsail/oceanmesh/assets/18373442/a9c45416-78d5-4f3b-a0a3-1afd061d8dbd) See the tests inside the `testing/` folder for more inspiration. Work is ongoing on this package. diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index 5440f5e..40d138e 100644 --- a/oceanmesh/__init__.py +++ b/oceanmesh/__init__.py @@ -38,7 +38,14 @@ get_polygon_coordinates, ) from oceanmesh.grid import Grid, compute_minimum -from oceanmesh.region import Region, warp_coordinates +from oceanmesh.region import ( + Region, + stereo_to_3d, + to_3d, + to_lat_lon, + to_stereo, + warp_coordinates, +) from oceanmesh.signed_distance_function import ( Difference, Domain, @@ -63,6 +70,10 @@ __all__ = [ "create_bbox", "Region", + "stereo_to_3d", + "to_lat_lon", + "to_3d", + "to_stereo", "compute_minimum", "create_circle_coords", "bathymetric_gradient_sizing_function", diff --git a/oceanmesh/boundary.py b/oceanmesh/boundary.py index ea3824a..8d66784 100644 --- a/oceanmesh/boundary.py +++ b/oceanmesh/boundary.py @@ -1,7 +1,7 @@ +import matplotlib.pyplot as plt import numpy as np from .edges import get_winded_boundary_edges -import matplotlib.pyplot as plt __all__ = ["identify_ocean_boundary_sections"] diff --git a/oceanmesh/clean.py b/oceanmesh/clean.py index 3335432..d1d4f98 100644 --- a/oceanmesh/clean.py +++ b/oceanmesh/clean.py @@ -256,7 +256,7 @@ def delete_exterior_faces(vertices, faces, min_disconnected_area): # Delete where nflag == 1 from tmp t1 mesh t1 = np.delete(t1, nflag == 1, axis=0) logger.info( - f"ACCEPTED: Deleting {int(np.sum(nflag==0))} faces outside the main mesh" + f"ACCEPTED: Deleting {int(np.sum(nflag == 0))} faces outside the main mesh" ) # Calculate the remaining area diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py index fd70331..2f536ec 100644 --- a/oceanmesh/edgefx.py +++ b/oceanmesh/edgefx.py @@ -7,15 +7,16 @@ import numpy as np import scipy.spatial import skfmm -from _HamiltonJacobi import gradient_limit from inpoly import inpoly2 from shapely.geometry import LineString from skimage.morphology import medial_axis +from _HamiltonJacobi import gradient_limit from oceanmesh.filterfx import filt2 from . import edges from .grid import Grid +from .region import to_lat_lon, to_stereo logger = logging.getLogger(__name__) @@ -87,7 +88,7 @@ def enforce_mesh_size_bounds_elevation(grid, dem, bounds): return grid -def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326"): +def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326", stereo=False): """Enforce a mesh size gradation bound `gradation` on a :class:`grid` Parameters @@ -122,6 +123,46 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326"): cell_size = cell_size.flatten("F") tmp = gradient_limit([*sz], elen, gradation, 10000, cell_size) tmp = np.reshape(tmp, (sz[0], sz[1]), "F") + if stereo: + logger.info("Global mesh: fixing gradient on the north pole...") + # max distortion at the pole: 2 / 180 * PI / (1 - cos(lat))**2 + dx_stereo = grid.dx * 1 / 180 * np.pi / 2 + # in stereo projection, all north hemisphere is contained in the unit sphere + # we want to fix the gradient close to the north pole, + # so we extract all the coordinates between -1 and 1 in stereographic projection + + us, vs = np.meshgrid( + np.arange(-1, 1, dx_stereo), np.arange(-1, 1, dx_stereo), indexing="ij" + ) + ulon, vlat = to_lat_lon(us.ravel(), vs.ravel()) + utmp = grid.eval((ulon, vlat)) + utmp = np.reshape(utmp, us.shape) + szs = utmp.shape + szs = (szs[0], szs[1], 1) + # we choose an excessively large number for the gradiation = 10 + # this is merely to fix the north pole gradient + vtmp = gradient_limit([*szs], dx_stereo, 10, 10000, utmp.flatten("F")) + vtmp = np.reshape(vtmp, (szs[0], szs[1]), "F") + # construct stereo interpolating function + grid_stereo = Grid( + bbox=(-1, 1, -1, 1), + dx=dx_stereo, + values=vtmp, + hmin=grid.hmin, + extrapolate=grid.extrapolate, + crs=crs, + ) + grid_stereo.build_interpolant() + # reinject back into the original grid and redo the gradient computation + xg, yg = grid.create_grid() + tmp[yg > 0] = grid_stereo.eval(to_stereo(xg[yg > 0], yg[yg > 0])) + logger.info( + "Global mesh: reinject back stereographic gradient and recomputing gradient..." + ) + cell_size = tmp.flatten("F") + tmp = gradient_limit([*sz], elen, gradation, 10000, cell_size) + tmp = np.reshape(tmp, (sz[0], sz[1]), "F") + grid_limited = Grid( bbox=grid.bbox, dx=grid.dx, @@ -453,7 +494,7 @@ def bathymetric_gradient_sizing_function( nx, ny = dem.nx, dem.ny coords = (xg, yg) - if crs == "EPSG:4326": + if crs == "EPSG:4326" or crs == 4326: mean_latitude = np.mean(dem.bbox[2:]) meters_per_degree = ( 111132.92 @@ -463,7 +504,6 @@ def bathymetric_gradient_sizing_function( ) dy *= meters_per_degree dx *= meters_per_degree - grid_details = (nx, ny, dx, dy) if type_of_filter == "barotropic" and filter_quotient > 0: @@ -500,7 +540,7 @@ def bathymetric_gradient_sizing_function( grid.values = (2 * np.pi / slope_parameter) * np.abs(dp) / (bs + eps) # Convert back to degrees from meters (if geographic) - if crs == "EPSG:4326": + if crs == "EPSG:4326" or crs == 4326: grid.values /= meters_per_degree if max_edge_length is not None: @@ -818,7 +858,9 @@ def wavelength_sizing_function( lon, lat = dem.create_grid() tmpz = dem.eval((lon, lat)) - if crs == "EPSG:4326": + dx, dy = dem.dx, dem.dy # for gradient function + + if crs == "EPSG:4326" or crs == 4326: mean_latitude = np.mean(dem.bbox[2:]) meters_per_degree = ( 111132.92 @@ -826,7 +868,8 @@ def wavelength_sizing_function( + 1.175 * np.cos(4 * mean_latitude) - 0.0023 * np.cos(6 * mean_latitude) ) - + dy *= meters_per_degree + dx *= meters_per_degree grid = Grid( bbox=dem.bbox, dx=dem.dx, dy=dem.dy, extrapolate=True, values=0.0, crs=crs ) @@ -834,7 +877,7 @@ def wavelength_sizing_function( grid.values = period * np.sqrt(gravity * np.abs(tmpz)) / wl # Convert back to degrees from meters (if geographic) - if crs == "EPSG:4326": + if crs == "EPSG:4326" or crs == 4326: grid.values /= meters_per_degree if min_edgelength is None: @@ -892,7 +935,7 @@ def multiscale_sizing_function( # interpolate all finer nests onto coarse func and enforce gradation rate for k, finer in enumerate(list_of_grids[idx1 + 1 :]): logger.info( - f" Interpolating sizing function #{idx1+1 + k} onto sizing function #{idx1}" + f" Interpolating sizing function #{idx1 + 1 + k} onto sizing function #{idx1}" ) _wkt = finer.crs.to_dict() if "units" in _wkt: diff --git a/oceanmesh/filterfx.py b/oceanmesh/filterfx.py index ce84826..efef4a3 100644 --- a/oceanmesh/filterfx.py +++ b/oceanmesh/filterfx.py @@ -26,7 +26,7 @@ def filt2(Z, res, wl, filtertype, truncate=2.6): highpass, bandpass or bandstop" ) - if hasattr(wl, "__len__") and type(wl) != np.ndarray: + if hasattr(wl, "__len__") and ~isinstance(type(wl), np.ndarray): wl = np.array(wl) if np.any(wl <= 2 * res): @@ -37,12 +37,12 @@ def filt2(Z, res, wl, filtertype, truncate=2.6): if filtertype in ["bandpass", "bandstop"]: if hasattr(wl, "__len__"): - if len(wl) != 2 or type(wl) == str: + if len(wl) != 2 or isinstance(wl, str): raise TypeError( "Wavelength lambda must be a two-element array for a bandpass filter." ) - if type(wl) != np.array: + if ~isinstance(wl, np.ndarray): wl = np.array(list(wl)) else: diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py index 525b677..8e3be27 100644 --- a/oceanmesh/fix_mesh.py +++ b/oceanmesh/fix_mesh.py @@ -1,6 +1,7 @@ -import numpy as np import warnings +import numpy as np + def simp_qual(p, t): """Simplex quality radius-to-edge ratio diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py index 27d74ce..d7268d6 100644 --- a/oceanmesh/geodata.py +++ b/oceanmesh/geodata.py @@ -18,7 +18,7 @@ from rasterio.windows import from_bounds from .grid import Grid -from .region import Region +from .region import Region, to_lat_lon nan = np.nan fiona_version = fiona.__version__ @@ -174,7 +174,7 @@ def _poly_length(coords): return np.sum(np.sqrt(np.sum(np.diff(c, axis=0) ** 2, axis=1))) -def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult): +def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult, stereo=False): """Classify segments in numpy.array `polys` as either `inner` or `mainland`. (1) The `mainland` category contains segments that are not totally enclosed inside the `bbox`. (2) The `inner` (i.e., islands) category contains segments totally enclosed inside the `bbox`. @@ -210,13 +210,21 @@ def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult): for poly in polyL: pSGP = shapely.geometry.Polygon(poly[:-2, :]) if bSGP.contains(pSGP): - if pSGP.area >= _AREAMIN: + if stereo: + # convert back to Lat/Lon coordinates for the area testing + area = _poly_area(*to_lat_lon(*np.asarray(pSGP.exterior.xy))) + else: + area = pSGP.area + if area >= _AREAMIN: inner = np.append(inner, poly, axis=0) elif pSGP.overlaps(bSGP): - # Append polygon segment to mainland - mainland = np.vstack((mainland, poly)) - # Clip polygon segment from boubox and regenerate path - bSGP = bSGP.difference(pSGP) + if stereo: + bSGP = pSGP + else: + bSGP = bSGP.difference(pSGP) + # Append polygon segment to mainland + mainland = np.vstack((mainland, poly)) + # Clip polygon segment from boubox and regenerate path out = np.empty(shape=(0, 2)) @@ -515,6 +523,7 @@ def __init__( refinements=1, minimum_area_mult=4.0, smooth_shoreline=True, + stereo=False, ): if isinstance(shp, str): shp = Path(shp) @@ -545,6 +554,15 @@ def __init__( polys = self._read() + if stereo: + self.bbox = ( + np.nanmin(polys[:, 0] * 0.99), + np.nanmax(polys[:, 0] * 0.99), + np.nanmin(polys[:, 1] * 0.99), + np.nanmax(polys[:, 1] * 0.99), + ) # so that bbox overlaps with antarctica > and becomes the outer boundary + self.boubox = np.asarray(_create_boubox(self.bbox)) + if smooth_shoreline: polys = _smooth_shoreline(polys, self.refinements) @@ -553,7 +571,7 @@ def __init__( polys = _clip_polys(polys, self.bbox) self.inner, self.mainland, self.boubox = _classify_shoreline( - self.bbox, self.boubox, polys, self.h0 / 2, self.minimum_area_mult + self.bbox, self.boubox, polys, self.h0 / 2, self.minimum_area_mult, stereo ) @property @@ -789,7 +807,7 @@ def __init__(self, dem, crs="EPSG:4326", bbox=None, extrapolate=False): dy=abs( self.meta["transform"][4] ), # Note: grid spacing in y-direction is negative. - values=topobathy, + values=np.fliplr(topobathy), # we need to flip the array extrapolate=extrapolate, # user-specified potentially "dangerous" option ) super().build_interpolant() diff --git a/oceanmesh/grid.py b/oceanmesh/grid.py index 670103b..5c63af0 100644 --- a/oceanmesh/grid.py +++ b/oceanmesh/grid.py @@ -6,7 +6,7 @@ from scipy.interpolate import RegularGridInterpolator from .idw import Invdisttree -from .region import Region +from .region import Region, to_stereo logger = logging.getLogger(__name__) @@ -160,7 +160,7 @@ def create_vectors(self): """ x = self.x0y0[0] + np.arange(0, self.nx) * self.dx # ascending monotonically y = self.x0y0[1] + np.arange(0, self.ny) * abs(self.dy) - y = y[::-1] # descending monotonically + # y = y[::-1] # descending monotonically return x, y def create_grid(self): @@ -343,9 +343,18 @@ def blend_into(self, coarse, blend_width=10, p=1, nnear=6, eps=0.0): def plot( self, + ax=None, + xlabel=None, + ylabel=None, + title=None, holding=False, coarsen=1, plot_colorbar=False, + cbarlabel=None, + stereo=False, + xlim=None, + ylim=None, + filename=None, **kwargs, ): """Visualize the values in :obj:`Grid` @@ -363,16 +372,37 @@ def plot( """ _xg, _yg = self.create_grid() - fig, ax = plt.subplots() - ax.axis("equal") + if stereo: + _xg, _yg = to_stereo(_xg, _yg) + if ax is None: + fig, ax = plt.subplots() + ax.axis("equal") pc = ax.pcolor( _xg[::coarsen, ::coarsen], _yg[::coarsen, ::coarsen], self.values[::coarsen, ::coarsen], **kwargs, ) - if plot_colorbar: - fig.colorbar(pc) + + if xlabel is not None: + ax.set_xlabel(xlabel) + if ylabel is not None: + ax.set_ylabel(ylabel) + if title is not None: + ax.set_title(title) + if xlim is not None: + ax.set_xlim(xlim) + if ylim is not None: + ax.set_ylim(ylim) + + if cbarlabel is not None: + plot_colorbar = True + if plot_colorbar or cbarlabel: + cbar = fig.colorbar(pc) + cbar.set_label(cbarlabel) + + if filename is not None: + plt.savefig(filename) if holding is False: plt.show() return fig, ax, pc @@ -393,6 +423,10 @@ def build_interpolant(self): else: _FILL = 999999 + # for global mesh make it cyclical (from MatLab) + if (abs(self.bbox[0]) == 180) & (abs(self.bbox[1]) == 180): + self.values[[0, -1], :] = self.values[[-1, 0], :] + fp = RegularGridInterpolator( (lon1, lat1), self.values, diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py index ceba4d6..06ae670 100644 --- a/oceanmesh/mesh_generator.py +++ b/oceanmesh/mesh_generator.py @@ -7,6 +7,7 @@ import matplotlib.tri as tri import numpy as np import scipy.sparse as spsparse + from _delaunay_class import DelaunayTriangulation as DT from _fast_geometry import unique_edges @@ -14,6 +15,7 @@ from .edgefx import multiscale_sizing_function from .fix_mesh import fix_mesh from .grid import Grid +from .region import to_lat_lon, to_stereo from .signed_distance_function import Domain, multiscale_signed_distance_function logger = logging.getLogger(__name__) @@ -287,6 +289,7 @@ def _parse_kwargs(kwargs): "blend_nnear", "lock_boundary", "pseudo_dt", + "stereo", }: pass else: @@ -400,7 +403,7 @@ def generate_mesh(domain, edge_length, **kwargs): * *max_iter* (``float``) -- Maximum number of meshing iterations. (default==50) * *seed* (``float`` or ``int``) -- - Psuedo-random seed to initialize meshing points. (default==0) + Pseudo-random seed to initialize meshing points. (default==0) * *pfix* (`array-like`) -- An array of points to constrain in the mesh. (default==None) * *min_edge_length* (``float``) -- @@ -409,6 +412,8 @@ def generate_mesh(domain, edge_length, **kwargs): The mesh is visualized every `plot` meshing iterations. * *pseudo_dt* (``float``) -- The pseudo time step for the meshing algorithm. (default==0.2) + * *stereo* (``bool``) -- + To mesh the whole world (default==False) Returns ------- @@ -428,6 +433,7 @@ def generate_mesh(domain, edge_length, **kwargs): "plot": 999999, "lock_boundary": False, "pseudo_dt": 0.2, + "stereo": False, } opts.update(kwargs) _parse_kwargs(kwargs) @@ -461,6 +467,7 @@ def generate_mesh(domain, edge_length, **kwargs): fh, fd, pfix, + opts["stereo"], ) else: p = opts["points"] @@ -472,7 +479,6 @@ def generate_mesh(domain, edge_length, **kwargs): logger.info( f"Commencing mesh generation with {N} vertices will perform {max_iter} iterations." ) - for count in range(max_iter): start = time.time() @@ -507,7 +513,7 @@ def generate_mesh(domain, edge_length, **kwargs): return p, t # Compute the forces on the bars - Ftot = _compute_forces(p, t, fh, min_edge_length, L0mult) + Ftot = _compute_forces(p, t, fh, min_edge_length, L0mult, opts) # Force = 0 at fixed points Ftot[:nfix] = 0 @@ -522,11 +528,11 @@ def generate_mesh(domain, edge_length, **kwargs): maxdp = delta_t * np.sqrt((Ftot**2).sum(1)).max() logger.info( - f"Iteration #{count+1}, max movement is {maxdp}, there are {len(p)} vertices and {len(t)}" + f"Iteration #{count + 1}, max movement is {maxdp}, there are {len(p)} vertices and {len(t)}" ) end = time.time() - logger.info(f"Elapsed wall-clock time {end-start} seconds") + logger.info(f"Elapsed wall-clock time {end - start} seconds") def _unpack_sizing(edge_length, opts): @@ -564,15 +570,21 @@ def _get_bars(t): # Persson-Strang -def _compute_forces(p, t, fh, min_edge_length, L0mult): +def _compute_forces(p, t, fh, min_edge_length, L0mult, opts): """Compute the forces on each edge based on the sizing function""" N = p.shape[0] bars = _get_bars(t) barvec = p[bars[:, 0]] - p[bars[:, 1]] # List of bar vectors L = np.sqrt((barvec**2).sum(1)) # L = Bar lengths L[L == 0] = np.finfo(float).eps - hbars = fh(p[bars].sum(1) / 2) - L0 = hbars * L0mult * ((L**2).sum() / (hbars**2).sum()) ** (1.0 / 2) + if opts["stereo"]: + p1 = p[bars].sum(1) / 2 + x, y = to_lat_lon(p1[:, 0], p1[:, 1]) + p2 = np.asarray([x, y]).T + hbars = fh(p2) * _stereo_distortion_dist(y) + else: + hbars = fh(p[bars].sum(1) / 2) + L0 = hbars * L0mult * (np.nanmedian(L) / np.nanmedian(hbars)) F = L0 - L F[F < 0] = 0 # Bar forces (scalars) Fvec = ( @@ -666,13 +678,42 @@ def _deps_vec(i): return p -def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix): +def _stereo_distortion(lat): + # we use here Stereographic projection of the sphere + # from the north pole onto the plane + # https://en.wikipedia.org/wiki/Stereographic_projection + lat0 = 90 + ll = lat + lat0 + lrad = ll / 180 * np.pi + res = 2 / (1 + np.sin(lrad)) + return res + + +def _stereo_distortion_dist(lat): + lrad = np.radians(lat) + # Calculate the scale factor for the stereographic projection + res = 2 / (1 + np.sin(lrad)) / 180 * np.pi + return res + + +def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=False): """Create initial distribution in bounding box (equilateral triangles)""" + if stereo: + bbox = np.array([[-180, 180], [-89, 89]]) p = np.mgrid[ tuple(slice(min, max + min_edge_length, min_edge_length) for min, max in bbox) ].astype(float) - p = p.reshape(2, -1).T - r0 = fh(p) + if stereo: + # for global meshes in stereographic projections, + # we need to reproject the points from lon/lat to stereo projection + # then, we need to rectify their coordinates to lat/lon for the sizing function + p0 = p.reshape(2, -1).T + x, y = to_stereo(p0[:, 0], p0[:, 1]) + p = np.asarray([x, y]).T + r0 = fh(to_lat_lon(p[:, 0], p[:, 1])) * _stereo_distortion(p0[:, 1]) + else: + p = p.reshape(2, -1).T + r0 = fh(p) r0m = np.min(r0[r0 >= min_edge_length]) p = p[np.random.rand(p.shape[0]) < r0m**2 / r0**2] p = p[fd(p) < geps] # Keep only d<0 points diff --git a/oceanmesh/region.py b/oceanmesh/region.py index 668dd8b..8714610 100644 --- a/oceanmesh/region.py +++ b/oceanmesh/region.py @@ -12,6 +12,63 @@ def warp_coordinates(points, src_crs, dst_crs): return np.asarray(points).T +def stereo_to_3d(u, v, R=1): + # to 3D + # c=4*R**2/(u**2+v**2+4*R**2) + # x=c*u + # y=c*v + # z=2*c*R-R + + rp2 = u**2 + v**2 + x = -2 * R * u / (1 + rp2) + y = -2 * R * v / (1 + rp2) + z = R * (1 - rp2) / (1 + rp2) + + return x, y, z + + +def to_lat_lon(x, y, z=None, R=1): + if z is None: + x, y, z = stereo_to_3d(x, y, R=R) + + # to lat/lon + rad = x**2 + y**2 + z**2 + rad = np.sqrt(rad) + + rad[rad == 0] = rad.max() + + rlat = np.arcsin(z / rad) + rlon = np.arctan2(y, x) + + rlat = rlat * 180 / np.pi + rlon = rlon * 180 / np.pi + + return rlon, rlat + + +def to_3d(x, y, R=1): + lon = np.array(x) + lat = np.array(y) + # to 3D + kx = np.cos(lat / 180 * np.pi) * np.cos(lon / 180 * np.pi) * R + ky = np.cos(lat / 180 * np.pi) * np.sin(lon / 180 * np.pi) * R + kz = np.sin(lat / 180 * np.pi) * R + + return kx, ky, kz + + +def to_stereo(x, y, R=1): + kx, ky, kz = to_3d(x, y, R) + + # to 2D in stereo + # u = 2*R*kx/(R+kz) + # v = 2*R*ky/(R+kz) + u = -kx / (R + kz) + v = -ky / (R + kz) + + return u, v + + class Region: def __init__(self, extent, crs): self.bbox = extent diff --git a/tests/global/global_latlon.cpg b/tests/global/global_latlon.cpg new file mode 100644 index 0000000..cd89cb9 --- /dev/null +++ b/tests/global/global_latlon.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/tests/global/global_latlon.dbf b/tests/global/global_latlon.dbf new file mode 100644 index 0000000000000000000000000000000000000000..21e576fb7eeaf42cc0a39d3660dd1c5bb8a7537e GIT binary patch literal 1225 zcmZY3J!(Qx5P;DKmSST%L2wDq-&2befe^59LvG##Vs)mNx%aEi=D9t7F3a+@eEe?b z=l=bUkH5$M{rkVW6LY4V!jz^mwP{SN>0-K?Zlq8j>O*~~5A~rw)Q9>|AL>(ms!#Q) zKGmoCRG;coeX1|@rM}dc`chx&OMR&?^`*YmxB6D!>RWxQZ}qLd)wlXlKk7&Qs2}yC Oe$wEp**YEnB%ZuxFpZoPXuXE0Q*6ZAKbc(C!{vW~RFMZre*T&3W#Q~M$CXgz;o#Qyf8&5|=5Hi= zHk#w`E9TggS7aa%>6ZTdBH(DkS%q3zH(Zu{PJCURALfV zyprgK58n1}a+kNqj05#C9jDw-Zmg$Z?4%toEG(xcdbwd}N32u0q#gQ~ln=P3xM4eg z)u}hRws_*q<{rBzZdiPWo%4f&Ez0_D-E<(>4X-{w^qFzO2EX*)p0Uw$!$6f{@AOC; zJhdr(b+@S-&Rp252XCygLjEQLnV?UW!o`2q-Wq8e67|9(-7x#~{S>wVE6f&%b4uCn zh86{EvAaoDcyH?aY`B{nT3?jv41Q*bv>VSNp3~g0itcOF4P8s59=%mcPwdZCV%OCHk*BVo@8W9$pL@|v za`a(r=L>G6btM1i3Ovc99 zeIsE{bg<^WXXTwtGA1|$v}j(?#zO&7g&&TRQBhRt;G0A(Eal>Ce3wsLH|KY5@zg}c zEzCpKZ^`J#JX$PQx&cGAZkrYJxM56>^Weup4b%`C^R1+~VQ+A+`4xW+jNg-y>qzK> z&usd-E7#RAT2|?S$7?ryoV%*j)=3RNNWZ(*$wI;Dqf9fgCslAKM~etago4uD>o0bB zuSY$zajsoJL5hC;a!d3&9M4KnNZUxk2h3kT&hKA~(wf_*lhY|!UqWAD6(Nrw(k>jH zETUipf8|)XlPnt7I6shWq+oTXa(}g%G-it@udDb*LFUxMr+V@vFx#&jyXB}DBVi`| z@T?e)^?Iz7nNv~gy7uFvn}o2c@%My@6BXBPc+Vi$&x<2FV^g;aQ8CA^c!{-y10U)B z>fAX&K}+?D;njVNIJ5GTn@)#{#qo{D!JQtEMWY%!>BZFNU1-o6|h#hid zg2LL?4;S1i`1;yF(Wn&%q(n*dB@@@{cU7$!(c_2a*{;*ZZWPQp;8gC{Dh!uSikwjo zrQl2Qy6N+#tKoFT-r<343i9sQG^*1o31aO>4o%&pphKUnluVioHIcOG+*?4cGI~E?2 zPrJE60cyGR4J}W&V@Q1hAfNcFt zeEo z9$%KPlA4gqKA9jr?~NPnU(~iYYl8iYOlvt~A0*v)X8GGy3rf_)R-F#_L6yJ{eL^5?noJL381!MkqG`ViZ?{|7g-AO0Tiom&6{0M50U57!d(-AHlg zax#SOZ;kaCTm5ia_wg#{PD99+aCtL#%MVA*HFi@D8v!+RVJSb&59OI<1uD~x!7b>R z(L|jeUOs&6N>`u>6eYU+(3mHsrTB1=CPj1Wtw!e~j8N_FMWUK)j~gob+{njCN9uX=5^n0qyehq~HFi#UQ2kLeCsh zjJwW-=?9=%)C;bG5_9O5-JZ8`YXCAxSC@jg1xS?74Z)cJWG$h1gY^FC?*Yip>T|07oF$w~HM18M3B<$l-fENbRv>?*(N>8m5OsTBes8I? zf;*3SQkTR6@$nTdN{^*A)X4@Ov*HRwj_?PL{$1AayxVVw5P@G^XDPlj!UpzS$S>tz z3BZ5@`_dQ|ZGdC>vw$Z@AeymXSW|t(7JSd$(A?G=fLHwQ+;Zo$18RfVK<2dojNT`$ zBa>wd8$+fY-sT4&$CL09VOcwv8=$|@-57xKFRN~tU$p~?{dqBNih*dcxO}#A+75Q) z`7=rE3dE{jp2j#YduaaMeYmPK5HEf_%C@P>9!8q{s(e-jVf{g?+r}Fl;BinvS-fx% z$~vz~Sw;si)m;eRWgdj)(V`SyWk9%hzLC|=aJ%2<$>Ky1b}c%&iJf$Y{??m^ z4++!oZT)BVvl7nmO@8THyDbf8pG94XPIHFRqoM|?30MnV7QX*C`$!2jvQuE6wQd||);R^?$A zSSdJ}5nLRMTZRw1Ze((SCvQaGg^UHGZ0?hg+j1^&xA}%jM@KMP`-`3~DszD`PW^dN zwh*L^hFra7?E>wEH;zPVgqIOiC9kEDgb`a|4C1Bb|Z%s;79ubO`cK&ho)#XNWJA$T2qvMInY-LCtt) zm1kLR7?DeN@PKM;y+Pw(an#JIpX&xp9`lTfsd7YN>sE|6g|zf0#* zDAx2G2|8%y3cQ-PMD6@Tu_UmLqvNx6y3|%N~Kgt@Hck9`#?RUKavD}7B~2v?}cLW z+tW9H=aN8!oW?;82t|pT;?%h`5}egoDiUQ4MW-CEW~s9zSbff2<(hjadX=VSX+9!B z*3b14qGv*Jytmdf><0-F?FW}*ABUovqjQRwJ{i>oMca<#j$-pxCy2NLFFb+8H9#!5;hL@$U_AMa|OAmhf$r?|FpxyK0 z;(;_g9{FBkvaAZkB+p5jehL{z8_i!0(+8r}#?}@-0)I>TRAaAB09KUis?|r5VW(;y zpVR?=T=hV#uEK&0$2tm>SZe(+c~f63mo6Dz+*!`$KkA1}dB5W;CrEJa&HM3fExtIH z6!7x&6bWidLXtKJ_~PP+;MnMoBpBXGeaZIR2f5BMcubxm0ss4(l-bWd=vg{4xxJ4B zbo0e0S{MlbV}G}pPKgW?n%uu5Equ_1&Ed*VK{D6}aEMi%^2Ry2JsmHJ`Kyl0)30=S zVXDL)>d8qGcsOTWms0Y=2R_dSMfu5aaLo(z#Fw6^BYL`n=M@PolHNr&9Pz}eq=?xT zAu_z?O}qIy#uHyL8q@X!kb!@qp!}JX7xHxeG*>)C@Y~w3l10+%UoN=+k&%L=A%P#= z8pG8_hFgV?B<00D@rm&lo?U5VFb(cge_Q5(h7txR^KHoR{*%-31B>pce!n13MxG38 zC&t3|op;CZ7T3p~gCtl=ow|5(vpX7hPp#8fPX-PC6{86v9+U_>dn{!e87{ULafTkH zqCP&FEF_cR2xY+f%w8%ck`FO$W+MZ0_ptQ+?NoftX>-S!KaS4}RP^+1-VmKh zf*s=N4Iby*(em`s;5$)-KC9$Uq+5AlV>NwPmOBZk54T@%ruV=GxZ}#sPl9*LU*gRV zdmwwez)it+S0Ej`Wb{_m6TJ?Pe-;~dg*Mr}uFG#d(DW+TWv5N9z^-)-Z#sJ7!H}=< z3@omYu{~Q!m*6k^qE+>#R2MjOUs6Qzt|w~qNDbVaafSo~_O^hdo=9yyaI8Aq893LK zXjyA|;%VdV{LpzP;Pj9yxs4u#Juy2SUE%~l<}BqT=#`7SRw@+b1d<9z1!O0i2dStg(Y$;6H3!(UDlFbEl!|lu>vZ(jIKaXckD2YVRP5~Aw$Gr*9z@n?E9c&( zVB+JQIcHk!;lzfk9>w=4$P<^HvO3TKBA(kEN~xnD$KhzJ!tV~S$TqQW-~t8Z^Qiv1 zsg58pr8EE9oPyqD_OEqAj&SoCu#@D7YZH-2iE&QQPkN`&Z%aXj)S9VQI%jwJ>^YGa zQy4@h;%^=9xkD~+A+6O(JdKRY`&VDBBmCvrx$M9U7cz>hKeRl*jj;bqJm=27CSj#* z;BvN-GYnrD4e<*jVY1!kbIZP z_xBf%!6`?m5;W5oRdvS5EoI>Yl8)fw#tmcNoiSm*b@uKj4lq5EyV`o4Gjglm-ySaR z2y<_IMxC!N<&`~5tF4`g z)OW_inzt+4uh@f`blI1ioz9qQ#yoey&>q$t&+!w!KkDbJ2<3EYEsxoMysRC&q+15a9JnHf3+bQ54{gtC{wfpNBcLY7(bBE zb@P<UP0L?o~Jaj@S_S-p|Z6g{~OwCp_)rX9LFGopB}dB;4)$g2RK?rCD+9^a2xSt++bhTSFA?$Vbl|$$1QOa(K78&yZ3Aor z97RK2B&1J%rhdWO1~%txe6W+)zgCYjhnk@cq>k;(-AvSHR0O~H6bx8H;>?+^7b?jp z&t9Y^<7W-OxpwrcJRoBj`C(LFvlYBoXH-5()TjLF*4ynVvVz3$b1zwl`d8vtYn}d$ zR#47l#kl^k8?volERp+Y3I1`9+$Oq+dK>)e>e*)rS7H7x-B&lv%vg>+HDdt{a>{+aBF7pM1?>IBg}yCWV4rpJx(b5#X~94xDPHV5Rf)IT0d%aTTzS zjbo?c>6;2ceS~@y9=y^~;Xp;@mkh@)c9=n+1KU3OKq_LzU{m2;Q`m4k-)pxc6$MRB zpH;~>0q?T)GHn~E$m5vn_9oRB45U96nFtc~mblBWZ$=w|R!iE6K_)5&6gl%8bTNb) zleSmuEve{5_dI#!gaOQ-VmqO~lZu=MW84gJhR|_*^_7}pDy~VJO!bd2g0q58&5m57 z;*dUjyhf%md@ji6?SD(f?K;KcFH%k5?yZsM66@S?_*wF<`Hv<*m-S*fH^Ci0jVdWL z>zTq2*}+fiTM7Fuw!_DHiz!s5-u^x;?STvL>oz^7nt~C@k|in50|nKOAJE}5h0Zb<^uDzt*pkx(4!JctkariCauk__+ga?e_{X%HQ30W^xXWeKv7i0{^ zim|_fJ-jgKz)b(+BqKQR;$5PjyBB5+)jam|Fa&u&uFJwvUPOK2CC_6M1GxV%efUwF z7wYcxRB{*7hl2?lY7fPD;g5cXIjvb;SYi28hg-bx_9sC)wU;_DDmK&K;pBynHNX0C zKh}l~w76#_5?;7;KJU=XM=fAK7>x<@o@mB=cD?xrO-O7~4&HUk6O-((@@|RNggDph zlPbZUXqie)o7}7ko37}aJr?%Fo8yZ5fn}Po%{QSsq}2nr`rWE;T+xIFxijac&U@gE z{o5bCAJ+o$Zyj9$Sss|Zw%MnhTL)f;-Zsd!_CV>`le@PP^%9YvmMY>u-H}t>IMOgj z2W~vtX5)6t9T_}|HT30lVCKZ;kH0^Wj0PV(W zeKrqN(p^*yTBQkd<5_#!`8=@oV7ipcum-&8*y#Is#T`HQi48^Ps>7x0!_Ju%?)YRO z=G2c9>Tu`WN&dD7cT}D$@R=m)NmXv2-ZZUtM^=O8m)9?=1ItFeoq6x6Sg|x?xqLtk zBxf97Fu)e5POi^`-_S2KRfnGoUl~@-##wG z;q9&{BtE*?I7|Vq%v{xPVQ@u9k9*TO-{b(}`)WIrTu}ble2~YiEZFvNr{s&e;EQH{ z?+0cwa8-RHx9BNntdN~y5|ot&KlYbiEGEte1ZTg+w>sgepfR0aSH<92 zY@T_9wG*l-L=G0Wi@<>GZ?mb}j`+raU6E%{1j?5>A3c(E#0!$24M_PSpent8hs7xe zf`2!zWIq>z9IhvPRpJiF+P`kt$V&hgZjPJPT(rmXHAn6-HSxm1t|;@0RC~;F^Zk-} zjt$oCc(Iu1WRI6VM{p-|(7{5O?pLj~_UI@0>(JU~V`g1jN;rDb?J?hNT%50-2kqrE}6-OS1}>y6Pj#rzhu2#-=-f8cfYPnAGNVTiA~YxY}l3? zxT@F>eKobgDY7`fv@#pMf7oLj{?Ho3<|Y~AA{p_0jEGc(wl$Jobmr8~G2;QdU2iik zSs~Y^6xRJSk!7~pfR)xH-4jgBqM5x5saHm zUa9jUMbr188gcGaV{lFK9lZb=r?lkMh*+Rh3-!iDnjk({oxYmyj5!K9DOr7bDvWI= zXNq*<&2hGD@Su#|>VIc>cem5U{A`lKP1(}*8wthp!!7V?y`GjygaWp@j`fQ@u)z4YW7^KPiYWRm-7Am967#c;?7eTLgp&t0 zh1Dim;`;YfG}*iBakxdg%Xz;g(%*Y_Q+!1QeMD2Hl6P66(oJW*k9O)PnD~lSCfgE| zGq@kyJk>yM2BxGdWtR9(VTbX%FB^zD#3}8K6A)x!yqE{4)78}zd|Pv`Jj5BZOWk2|~BVwsNk zX;}$yP1j>R z{msx#%e(l!sslRiJl?`&1Sm4Y{AKlC2V5!Kx^<-%a8T-v*4r-*SP`|RYeCB#e@x9L zOJ@-AzrNb7nizAuq%Ou%Zs>$L-v#phgUnGpRVwq!uoKSt`BYCEn`7p?4>Wo!XUtn5 z=4k-@DYXBcrLr@=(@oncNH)iIg`jH{r=4(Db5vIJWplh}Te^1A#|b|+aR{(zSfHep z)X3^%PN>JUvAmzqqwSSzuCHQr#^!XsRZD7?=>2UhJgveRqu-s`QW<24Q&x!*f?Hg0 zW38;VYKkRVuy!=ZaJu5q=er+it(GYC;hQc^+7+#ZJ4NjRtT45`^2FuWuDDnwd95ec z3V%|gte)DE@b{dlUg{PrWZ~hh*ltC_Ue?eDqAON7_%d%pgD?rh*|J80_gdo}j~o&& zg@m<}$|C~v*7)$}sR_0%B)n(abCZ^4gF?SH-uT@@!UJ~%y3f|zpvtRe-|IyFXZBzt z<8zWN(*26R#THM(YRB!=kQ`eqy)Gp%m`%cVQgZZmX*+zJ{;-a7h=dJOUl-C3*r8x? z?*>(QG9JlhK5F&b4n5d+2yU|{qifjlyTXU;aTR6lp-Cd&0kEYpUS>67f)sG+RBc1y~#s})m^nCUXIEDo)te$SjyJU4xDAfV~dtOf* z;WH3Fq(A5WzDJbFABnSDzS!*FBP-7E(At*y{0~+9+e7f@$lsP2f8e3UrqWkiY|ypc z=-Cw;7jOunuv8qf#U1?y8x405dAVS-^u^tFc$cAkC7_L{`;I<2pDk;TUE8?*I~-i$ z;)JN^r`z^e;A?&&yv-FpM?K8vShC0T^OWp5HKK0Tx~g`_#~w8$qR#YeB|&9CD6~~N z-~%=;@5=oo2=_eFnRkk6>45v#qO93E zNf4?d?zJt(0rx7zueghnA?ae)ndmk4_*-hD-FhIyqgbPXdjs|;T+4i-HNALQEg~m}bY_Ay=d-co_gZIl; zzUw1Hy|$sr$xV*fu+Filah)3!Hd!{!4>)4S#}@<4F>W9&%kiaQw-ahaDmeR8xdDAV zC1&*rXLOubvHCXW23uPkl7rv6VhrP9Ln;RawmrJIa#)g#>^*xwkXBK^>zA1p3le$Y zyw_N&Jq7MaT9!M{#ObAnw8-9TZ$pJ~!R{wh zU$t4L*sD2F+@>r*deP469 z=Llqtm+(!*Vl5diKc)CL37O&mzgfZ! zqRwlUS{bM84|wOED>ad*8}|%8=ceZ}$EdT<+mmjQAlkYoxo))u(tmTaR!Sj(oS~F? zMx-UKP0#-x8$yEE+(E{qWGhTNeZjU^mju@9I-g#$w89@x5Bomsc7?g7X*o4DYg8y? z%`4d93g`FD4R?fDW6*+Ii1rs3IJ=SQ_WE#Z^f(w>olbFqb2ZJy<=!^`XM6r*UkH;! zbUHDf_+zK;-`Fn4PV<8!SA*;}rCMT;>=n)?tq}Nl^t(%6;#&M;r?wK~4?KT-&!NOP z7%UW=+fHw^#PYS1J?jO6p>HT}R8GMXQ)#OP&WzAN-BZB5FxV11d}xDm^)!&p8*5IJ zx5NMj$f3`ofzOWn*5xY}7`W|XMEgz}Y(B)yUU9$z&EGTZeG*NBPf5c{6SL-cXzWAZ zYD*f(HU4s1YBR^46sx|s#xz)y_Y-PO%(44L{-V%o8suGPUwV7a92uLOO2*V_FwOIh zQBuMjMXBvJZ0-3T$PANn`k$=i= z3+GosaG5~)@1QNNQ zU*;v7;pL}?=ZlPkV9I(^7X4{6bU5F4zLh@+E#D4aN$UnO%V80>X=K{dY zS|Gmn2jHY_hU+V<062f~mUH?z;B}eadoRcR!FZ)trHIuWh4l&)g^T^+4*iL;OmA~M zz1?G(?BWmmbPLaR5%whM{pk^f1wW`uo#xx*VUBhWM66bI`@xCDCkH>Z06OV-3v8M6 zgN{wJ^{IZ7`T8p&`#<8GTexY4FCg3*14l)r&m*id78*e6sdmQVD`fzmEP$ zFI2%=GuLCQtb<@M=KGd{{`L5Ztj)5X8U&0y>Uz|sb?ExWlopp51Ooo-D+&kJ;n|y{ zAB>qnkk9%mFgblK-a5pxA)+h@=;@Mf=N?x?p*2i4_lfsI{t>^{N$GXS-m0a4eJ%*J ztY^PWq_4wQjY^?zf;1@kwr+$nsf05+X`%uvL7=B#)Xe)@8AF#@ll&)wVC738cj|yD z1`D|Pb#l_6WA(GtSaWsU;&WN-yc!L%61EPWU8jL5agssLTxf74_4QMSXBzn7g@_mL zRvM%NOVInU4Y)8m*R7F91HNERAFCr9u=aj;=-#U|h#58I!J``Jk#X|1OF0cR?F0?Y zdel(Nc0sN74h_VWLzt>_RPfE~7uje#idFA_x7;T9JA zWzYAtf41;Y<;SZ+gEaVB^~5Wpd>wv^rR86qrh#?8^rf7q>yZ{k*L>qQ4X)KmlUwIi z@MN}Bw4`h>1nFOH=RB>7@V#R$Ryi1Q6zMuAUa2E@oW4-LX)sJbmO1_5iw4mzc{0$^ zJQzlg79LOJ&_pp8QwIxDFqrjrPcVvW;f*)USN&pxp|v*V{>KU}{FOMvdp|cA<~-*p zH`TQxFi^y;;0c(dO9fHZ5CN`Cm3ElVUnD1(#4)wS-Y-xg#C_uls$G{7d5h$ zhrbRI^E>gGw^r+6pPK|<0DTA?>C(9EnW~RFkKYpIS`z{~K@}k;9tOBsacz~8TnKEv z4|in}3~|lnguGIv5ZKLn-+tw^(Ldkuec8SLWR)?lzP2`N3vs^;Q@n>dOt3cQu==7* z2q^8BYQ8sPiW~;~{PTPvux?XA_6rs$Yf#3| zbt)Ki*Oq^Z9=5=IHTKKjHU~q`ugjlGBrWl!LA&UfUGP8O_h0`;xI890I(n}^esVK02-EM~pvkHG+FE@;M2QGT&^oWp@Sl zq?^gnM4c+Pod@9eZhF5?{=`^m&O)W5_~QE;(q&$ z8}#bA@-u`+&caH=KczL^;d3Y$Ob#~}%1BC?Jx)^8Nh`!o_7(V$` zK_rk^!r}eyZpaY(o5AoJ2_gc`6;a*|Th`6&5GUdjA&-t5@fX}Mtk$sl6mbq7Y|FI6 zfs2B3!a?c95oCCk!+MV4BLy2TW=91jlA)g`Y4x6RDuxu*F23DGTpPeqmNa+V5vUiI zb$|?id6N^dNSSy*2Cg0*FD`9&+?;%>+oXlKzw8z8RL>nFCCgrC6Y-r+Q=Pw;fcrn* zh2<+{v&$3bs(Ei7&|A5|S*GZ^ru#%6>i43uR(dyhAD5%-{+Z|t{$bv?B+^}LkK|-d;RZX`U$?$!fqtSvQ2@6AacQ8IB{gZWL zEklX%eIlP4TKi;+qAPajK7mc!NYMCo7f&V;|J&~M7#MaTffL;=tB_#goa&hAuxK6$ zT6Wg(exf7eaBY6a`&bg_{mxn8Ce}N$dQ$qV0SSyo73qwLbG3(8ZEhzmxkAz7Qt}He z5=J}dkmoiM=inc5t_cu%@rs{E#&sh83oB1aoZUnEr;Gp9H$nvom8H-Bqqn5&<4=l- z80Deg^u(gC0$z8XeNAc(hB0~0o&Mwx4S)1jgBX9{@4dgcl*2drQZ5IB(3LesXENk5 zHOo1F*IwfN616&MU@MRJ*yHV4LV{sc(sG-^oGj6I_Ws_yV=&-b*P%7*WzYyL_A06p z_DM2Fy^N@j?aBRIHNg=KjvSj4>mRSdgcy6jZQp5dX7;(m?<8qVKKqg{td#~F!_K90 z$Mq;xB8m2!EVn@3pdVo*4^oSHHRaQl*XvJzm-H#lXpw29qrDVEg4~!vPj$ zt>7RU1m4R+0BL{KRHOJzR}=*N$W8|XAXQ9afxck9Sp*hTQUs6+;}xQUW29*441Yk?%vhT zhZl@^=30${!TsVXlH{Tg))YMem4IN79PMS9yC{N*Iy*1#P9pB7UN)p{6h}?bJtv$9 zf7r6GCUO5ONj$AL_}2GOFhtMc3qEZbbo6&;&bkr|f^27&ugA*b!wUbo8Nxr%;`&`m zROJ7?3xD;Ucm-I9ap=!`^O`=(PHw9+zH|)qOANLJnG30(`g|RbMI+yTlje^^z-pCC)7E9&cv|O7XHrVz{a>&)7}m^#oo7GnsP$v zt&yFYkL=(@rL`n?nGO?N{10#o>h6M-Cd-(%XZH>T|9` zmX5f0USIF+13L(`x|#T0&=Eg3rkIDS+Wn(dto6$AJc6E9hKxy}mk$5-@_24=5;z^v zbKX-$bNB}++cQ=SdE-~MlI9C8t0Q-f zg|9SRq#KRk5BCL+`k)Zy=i?101Fq(7vG;|X2cinLPa4er(rPi~T{SN?Yj&!+ewez? z7gnz?T|Qv`%j}Z%)V$ntU*K30&}QUZG1E9JD^oh_3qj}3dPT0VKz;SrWqV^kDB!l- zKE=uoWvYH9AyIy?MQ-cr#P>X~_OP6$YOx=<8WrB(zr+uB^O{RFdi=nS`_k2p0U>Z| zmKAB`^@os(=C%kWF?caGQqk=Bf4us@sP{hos6W&uIn1B15rMzFob*v;^@$V$+2weT zZJGW=AMVfY)0F~nYn#IEBuL(ZW9k=XP+hAit-2UGo$YA)-0f5saG|Az#q)i z;=HWa{xZA8aA<5zg+BTkHU^>U?~hCj5Bk7P7j{b;BQ-IkLh>klU`4R2v7X2eWZF@_F0f9QE~G-m!82TIv` z4{)I$teQ(NZsz92|H_+C0XAav{n49Vd4r#~HwQsL?ORROrvRZ9C3;8dX>cfV*S&!Q zR{zr*F#_)ooE0tJ?Kb&;xGh3=oi&KOK_sZ4=c<=ERLrc|QTH$i_>?zM4mz7d?ws}O zyi-AtGoD{GDQOOW*P4Bh+G(s~4(E(MmPrzMyGeF$aoS09*zv+=vo4X>fA7~yC80U! zK6cB|{~QEQ-iqwl`^+2ytgfwk!%BnU2}x6lwHDC3ExLD9l?Gxrsj?b7E#RbGN(!DEf?6#359Fc_iV|!_^H!4x2?Xe|XAa%aol0t)(K}v$A zhZP)4i+mPx}H0>-%aK)x(AA<}lMRa+aHp@N?_GR3bVHv(#P7@y zdzjDmRMM(s#T`jYHcY%A;#pLXxdgmycED$ zv1gFzmr`8AZmu{+xIx!MR{c&kS(p2~3y1muwJaK+5k#f0ivy@YwsS)4^&kJpYu>F`Cf{;ywyUAI@58 z@GI2KEADiJf-Lgh)$67j7JO@pLYy3di;am(T)n0t!!+_?(Q}7?yi0$#*(f{K0oL}G zYs;1q=e8w^_N_zq(0|T@lWii<%wxbq{z<4k{Fk2ui-~U;>Hh?KbUV!Kthmox!)ygf z`NM7lypg0S$F;^D6w)N~}-5an&EWT-zh7AZ!G{)CG z)H0j!f7pAC>rG1kUFc5Y!7o)=5lO7_vCPzrm7n5yOdxh zM)aFycpSO(C`$!RP1pH&(ClDU?%>_l_Vs9-ekM=1(eB^g+ZJ-Opg&VNGo{HN)|Bbri;ovZvb{mmf^r~8iZ?AB zdo*MAM-L>3@dvUvijQ3G_Xo$1BV0eK_))E>oO17zKcq>JtqrU=@$X!Z7h4i@vKVn> ztt&de_J{pjhn9D-{b<6$OvjhacDzxA0R z2mD+6U2OGA{b`~;CKfIFV0qHa`I-;?t$_fT|6NS2lOHqtORJmuz-WBRSF^|kdQXRM z0q|VrjA(achgm_2?s2K20O)+$Fn&DYlbJ?l^owWS0kDl}?e3a;y=LmMzAyWE1AsSn zeC)!FH)f6}Z8J1i5%s*-xYx64qYaH!OWto>0$?9c&0;V2OhdJVEQfne06gOzkC9#) zX-K9!w&QYl090~qX5aH;xZ#&o-dpO002sPmu4~CL*Ra`4D>-Ca0MPAi?RU3kz+iX% z^=Gn)`XJ?6;P!p1P`8-vS7CktT#VI>e9*y%qLVW-3BnR5SoAW5aJ0^z-3IX$IjB5&{2Y$C|B_e!bTp_Q)L3NTCzN|H_LH zAty0X|Hyq^`>;-GfIc=@x7yh}B5Dj~KL$Cil(FcI?_|X}(*Nu-7lHQ&-ZRdgveVdr z`L|wX%YAl*)vtbS-Ev3k-?NNleW|!b9SnF`TrSK+oCg|R&}3`WL7D2!BHGR*_`8;x z(ZP6f7zr}BT`H|Vq5E&GAnSuJxj(g0Xf5w2V;TwgX1MC~LNrk`{F;8xZW3HN(vhfV zuZ~>`Vp(a&NdIW%-?U-#ZeJ24yyuN}*r$qDi=T!*5+%Vyha)sKI>K(qD>^n3^!}HA zVoQIv_P>0W9Xyw`yF%eZ!uonGQ(JUPP%~9(3udgErMZ-#P`uRMBkXk34}qA z68D{pzpU|bM;T|6X&4c2oFpaPvPR*3joUVbgn>uM!L`1Q);QsDCn|7v7>u4|-i)nQ zxJ9d=fbY`((e{_tPZv(S%IUJggY~_Fqh(>x4d-8Sd0OLUx+aE;=fc2fb7kwherq&; z`REr$hru-`@x;W{HYn0c6lQ$GVCl-5J4Hn{s3H1ej0M7=sETA)AZ?5P-7mocB88!6 z`=dXr0$h&QfA9oqS1Nlamk2ifeC)Mtw-2a#UEJbyUig3dBl(AE1h^K?4Rud=!?xrd z%#Py%c+7V9RQ7FekhPqAJ{ZD>(pH}q9rt)c!k5AfT@g+kRdC%9Tq^wVcdvif5;mG+fxAL@npfu+_Q{>l`2GQub8MbI281fzv zQ{2n^Pgd4nUDhh3aO0lX0D;g1A2?*FGu_uAfaOa5U2Cj-{>fJDQqmj?oe%6=F0fT| z7Q+9^lMpcvF{=OQ#Y%C)>Mzb7$fQv({+pWuzwLsa@vSG}-Y7f$fEU#N>4hMH_Xi%x zN>j$Ry8XNN7q^?c`0BZ8Dkxr?qI@?aVQHq&AayGh!gvMm@Q9Jo#pA5Q2URNklA8zd zQ8G4^Xsc(9Q{bBMw>?jEDCoy#o-W=<0crIkxhdO-@4uN+jUVr(K-}|gPD|o@=)@`^ z<~BD9sO+5aP~oHgqcz0;jKd>qcXYCv_fzMmfM8kQ6?S4Th5PfoW(_GYX&j=sfv8Jr zrZ?ovgizpzP515dc~s0j&9xF$Oo6DJC9@eTD&DU(yr(usf%e=~SqtL(7@o@-d)u|C z@Rv7m1HxRES>O72-tnM+a7Q9M(yr!s(gu`xDj_JlGd= zmA#1y|K%6K!at^|{14yT^K);?~=VzdjI3Q0D<=h zW}A;+;isoy+~tK7*B@G-EGp~B;!64VtWZ?qn@Zxl_Hy63_nG&#VXcWQ|E?Me`md5x zpq)Gt&BrKXkyM$c*TfrD9o})0=42g~3J}Gu6uk{+@38CJsGN zm!bQ7$%EkY|8Fm+L{{BxXQd+jPu9IbL?5%=O7h)|M-;5tK);+BstaOM&h>Y8QE)^m z(@b_$2dc{kzwzo*Q1M}th@6uS#L~`R-tgTGGohhP~W zF=GDkS@&~i$eCPuoQ(Cp%Kg?JDqkh2Can?1l9c>SuP;0PkI#RCt3U84cCYQhY3YBs zf1fpOn`@hyv4?#p@9bH4Cxah~P=@N|0579^dB$JMp-$@FZTGtzK%MiutmuL?YO8Hg zy<*@9u>;Mv7cWTsTkG%IMe3zX?PN!I8n4l@#zPdHc39t0Z6v;b6}rqVwnYq+FTIOe z)^dOWUW^=FlEUSr3}TCqHFSfe_Gdpa%nx6uM63)E{XNB(?MiR&^hM4LSMxLzH;}y|nbJD!gGbLk ztaM;-1J{;aN+}6G|7iadVj&l#=Ywr`pLUNgk>Tv?x;6hFvfcxp%J+R7mr4rE_P_W`G={Px1-Qy!LqBQ zYO~RAh3fLsfg;yGZPoccBfr_@UI!1hNx}A^) zicbrX9s@+67k2u=IVbg}?xeQ(x4x*&ek~m~zs6zNS&nqyDGy9JbYX@{!vum(-+BC! z#si}ka-VGU`%mABl#st40iiSWGuJ&Z7Z;uR{8@DF`afI}?E2S2{bNxc*tW%pbrVHX z$R{4WcTm{_dp=)InT!3UP2we$a)eJf$OCI&7-^b0f%K>U&__&JQTe^kt`px<_a>eQ zrgj?K*&*!U(46P>sMrIZ4iD{+dn@?+`;zwA55`3=GSYg0SILxh4J{|IY8xwA(Rlz6 zbC|c~^Fq4e0iF6MM37u4F=Zm-0}}1xZmS>LWq3iLL3VP-T_Rk594Fjr$PXUPRo9uw zJm3r0ZPkY!{E)P-?Iypn2TZ=s<;~wB2rnk|cXgT@!X zpYG*`n5th{?Df;$aA#gXZNl%refBS<(XS^x$Ax>pvNw=d+I?d>X^1V$9kUb?^@e3! zCr8Gr0Y;)l!ab&{^{vs)8!ku2gzP?NfSJUT1?I$iL$o9xYqP(>A6^W5K4{w_{LJUi zeQt+%!_^}Fm~-KVSX)%!@`~p53m$uX-1zOqU~`O0M8p463?7qL)Yj8Sze%+Am;Lk) zP&3EOdOgIdlJHnm(!M}gG(&MG+E`bRe;J9ktevsK1PdN}x8N;g7GZ`B-F8nA8n(g4 z^!r@(R?Ys%LXz!Y{zxxNqrZ;t%|SjmGkia0_ip)b zpY}@``|BA$d9KZIn*!8(-LG*8vxNo23-LQYio=GL=%PxVEs$s}&-3=1jM;+KgX317 zXg(p)BGLYM{^2cLi!J1s#Q6vvR0Kxz>(~6yzHF}=bt}h+>>pYF<%dN0rRx1Hcd&Pn zmyMbyx<@TLL`J3yo$8sDC2^i`#lf%G@eTYQPb9h7yl*nb7)K&}na(R~iy@xKe#l_*<*Y6| zBsdOv*#DK~Uw%ktE>wQc^KP^hPEq%~;I(m7)r%iSAXV9Lz>d-zc50ooJUnFfKl8k6 zeCm^-N9dl6=d!}SWMkMG+?lzn%?nHq)Gxo`F@g6cUk+uWd(BrLIkdg2O<>#S2KMs^ zK37IOl1ASYNW6}l&}fzPn8IP}9Ueya(LH6BIIp#z89ZausXH*_^@sP~pv=(#WQVu@ zxR$NS(hTk!ZHd(X>;=S$y-IEOO`+_jTb|;o7gQXjw5gRc{UbZgiWSFTjtOY*;#91E z?FBb->~4L)nZRM@NjkDYFHjeuU@ERP23xxq#dOH-rLIh{e@@C6eiDlLcVK_p;@|qB zHseMm|5ra9932WB9do7PKisScq`?Wdy^>j=>?P;)g z1HAuA3gJ3GX1aQZCzaL>rjpukeY?hw%{l9y$b94qX_{k?qvQE8E5YdM@5fxBdp5t= z<`5rd*)e|Xm7*)`mb7NO*Srf8=uvQ_n0A5W+^@>}Lv~?#`ghOc^IeeNSr#!bjTa-y zwnI5&rw}t2_VuIX)Z87HKeFY;eUqi#bb&Pe%YmXCoLGWwZQpr4R~Vgoow~uojXkQ~ zQLLWg3QTH~ck1?WVUujDcJpNDS9M-D5P%|6?<8*Qsv0pU(8>+^ zhN6vIl+8h6NZ#J=7y)WN+^}5Pfnp9^b$;H( zfNYczAZgzE7d9x=EK7UI9Hx%8U#Gwb@c7z~&B()+P$k1y(1h%zPZoLEKCiX@LyJWF z%W~o2#GSSvd|2zj9GMLoV4g-5vjmtKtByacXa!wIj>cv+5s)7t$9Dz+OE?oG*S~QB z*^1j;+`*-{{Lu^XfMj|bV_nEgq^M`RN^<# zI6E-De;d-32PdXnS-OMl$BxWryU7^A@rhYFt|qTPycZ*0^)e14`&H|5HU%>y2-gn! z6t>+P#vH4BwSF4@ku{y?v>5K1Aq>n;R3}h)Lo3sxD_utop?Q#D%>s?X_W>)qG-QVG zul%U`zf>x}_dn~Vwi!P)%^>`*YI%hfS~VNX4`n z#J#LCebtT!67IJ&eahY)1n@6Dw^%)8{s*_ZrR_|_dvjoH!*w9fKBr~<;00DcfdQWfsz5)yt;TznQ@z~nV{=uKva&v)_3r;*C88a+7mFeP<;p8t)MXWcL`%?u5!!kx(lbSoe>5}U zh|@boxEUB(w2CAAq5UP=>G(IaR`VpXM9`Eg!TTAJ$*6KX*1Sm&DeO>!cHIsWm%S+d zU*|U|UOA;dvMnYK$d=nHf`7}9cT_AMoBq7Ne`c>D@U9xYOzXg7cT^ciCN~w}UwKg# z_)*#WtKJ0Rs7p*YMLnSn(Dr)APqDg z|CI*`DS!(6s@JYt%w!ymR+vjrYfAsB5twZKMkhLFiBWMhk{Kl!|Nfr8;C?-%cO+O6 z?t`@VwzD)Sp1=zV-B1&&Ke&Qw=h=M_tzOC{yS2@~v`P4gq#!E4=d)|3+zdKPwvZjf z@6z^@9nxJm`)VWYK=G1s$%(7~?W;-inV2j4A#HLyWWu0)yhDT&NVxR`|6KlSyueh< zN%3LL7Dl`d>bP+WfMSb9-Fd~owA^;kFSI(@!MUT2E%9aie|XVe*5Z0D$_v-nrdOiS z`xm}PnK(8%pm_4rcoKu%AKCupkLn$O5T9*`v96O&u(m7{Y(2_d?bS~D!=zt zPLw5rL&)}}h3V6P;aYY{(oqDd$;ROAv;p?s-QNzc2?|| zGQr~7h&NQdP2Bf-lO6k)9}-1`BzHr&>X)= z|3~g;tbK*s-v5QIRS%hwyXOmfttVWo`cyDiRodYz*}l*>kY+|FisEnnA8it?ts%Ae zoNK;7>8}zO*Q$)2x5lluelu+@!$OY%ed&*^R06tYGKUj4=f4yE=7#O4O!&R zq=|jw82yX;v03whjkOcD#q2h?Y7qX=j-zLo5}(Jzzw{6bF;srb#lGB_?|r}>KK)z| zB(@5`iB^gF_tWm6Tt_zl^^@@L{z#IG1iSvVV<6K)5VRiU6>*fggL>~f7vTY+{|C3J zGOWD264`?kv#TUt5`xk_(bkkPe|h~&56L7>l53j0{3c=21MLlJk4dWaH;o28f044y z6VjMC$!VN6e&17)JMBn@VI7^ zz2m=mp0o~=U~#^(C#x3-K;-(wA|UAY2UgJtV{%<|!5ngL#Whmc{-H&pEfjf{amO=P zY`0MeE#FH5=sC(!?izE$#IX7KZD?QiFMTAVBuP#dWs$L|8;)3+wr3LiX2ZR4{Tk>VDI5=D6Xg7`_}wzc7JFU z07c(P%obc;shlm(aQMTESLAH*1d6kC%_eg9ttO;b|36ts^bu<*RDREgW-nAf>_s;H zhrUu!t)R7QflAsDLRvKj_}Z`F}fs12zxUUREa|=+j|)+u{G8 z{Sm3ZHGOTmwLHNOdH9y4ke4_h;l8~RndI(e@c-bRP@y;YVJ!q=2d>ca_TeC)bgqIh z`Ipzf^pH$4B)Lq5%eoFoxx)Cd53EPs6k)j|D~{o|D}3eLS^I2@^8bua5^S1FSpnWn z2_g=b=(h&_g(cx0wC=ydou&l$#7|ljp!0~+a?A6ZS&A@nUeuHV{_^^l9+F8GmEYsw z(z!!UA?WV>P=nl@!JG;RXaD5B6F>yM9jC6bN~!-({~F%FDwan`61f}$3<0G z@mkwajPAYhLNCM2nAHCtTDNra^D59ixva-CLk@a1;NN_(Vy{Mo#h6p#RAj2~FMTAV z9I5@fkCKs_$`E0D!a{Z8LjjE6QEQSCBZ4c(w3bd0KlZDCk-n>Buh$b+q)&5v*zv|d zTmVZ;>p6iB`HM@UMZ(zbE&5D|=AIcto5^8caDp1}39vw|P{mmXpzkIL`$7t7Ob zy?xmp%VjU4U8%DMrpkw@<+^s*om=v}edhn&$CK7yhHa~APb}>)E3NLq@CPUknvbbA z`x`qfdu6k-XR8fxEzE~y``KeB+xMo+pu2q%?b-Y@d>UfNUr_ng+`BFIe`v3Uo>%Kf zaoe`2I3*_vpg4H{$D1U}zxfENyc5g*`{(x2^JxLFY zqejl6-GTrH1N)^~wEx9N zGAN<)yI*|69Y6nbi3rOEUlt9|31cU>FxITLptX%p@^IAm|N1MbUyxuwKTsr$uL)vl zcx_^*nghhGfHK~ zp75k@N&Sod*8dr2Da|sUsw`uULt?sfRHHYS}|&q#r&oHFFuk%mDE1x4ry)C z9QOvNr*Xr5t&G^s^IQDC^7+8FoVOi=y8qigQYG{=n=C$mU`G!YpUXk}<$W5LJ6D@% zF%m5j?mp>)E?P<-C=a6??C;tBhqnLC&XZU7Fkt`ELo%t6_^W9<#K%j_?+vPszdRtNU$W_2YM8t6bI1wm9Dm= zrgp@Bl(vpiA{~WttV;Vv<6mC?(nB(-qw;%wd1n}|aO9K&hTCs+mpJbvP_OTKE@S1rMBuoAZ-CzDgAIYeJN(^aVLKfqPy(WDeu)??di)tSH<<&-8r*=Ee2HX~mW$DiE z!NPNCAG0FA@!O{uWAu~v{E^M%=D1Vuz#c4hLOYG9i2Tl9>lqQb{4mF{tgCl?ZNN(F zx^srVAErXQXSyTV27b%=N8W#baf{@?|N7s5(7WIN12$l?d^O+E&JWZEjV)<byDa=|&zWKs;m7EKAHLc)n^6RbHRlabOqie~M-3p#&ny<+R`@)d? z@ObiDOE~TG=##IdF9ckd2v}v(iT}ha#DRr;mQ1@fkj*l=j&S41E3N-@fpH z>BGm(T@uY;UR8KF;;|1X9zU~lXvq{@3v?=C&ilZC%XnN&qbX<(r?Ds``M{_3ej*Es z8_5)LytF>r2YimMija?(f&KgIkw3$HpvGe+H06LfG%{Ne$PW2HU!Kc|IHd*5ydS5c zx#$Di6KkCmk}V*2*Hh`l`##_;x@Y^&bqk>T?mPT)-UsH7dMkZ$wS?`1_vJ40_(G_p zul;uvH`=bc?6|0@FR0~TiYqm+f^zocbK-|leU+`n&eoy$BU|*YQiS=!E%kn$8B=R8 zC{G=~V&)6y%c7m2+{yj^Oa>x($ zrtyaH^-~qhHI`7rW_e~`EAmS=dh%E_$rAW?OY3zbzu{XICSwd-C~oXss;aagFDPiU z`EGyB0*bwBRX^x?!RNE)4%?B0ZBAP>O6^+}dppzIIhjhdz3O{oULo4KH&z zRlCjJjLi#7%8Z)Cksg&YX~6Ic@>}i_b&?~tAU$npVZ)O|vpA8oSTur>zWRdQz* zbOz%xI4UXKX9T%?`cuzok)PIwmG}4kj6fY*{m_T}Ej%tK%D$X6glP|4A6W;z7!tGM?~b$1{0p1%2>%xcug)m?u~~8cFlN zss|%_@8#K6J%HHet8aBo5AL6Q_;`Pz2fUOGkUQn52l&F(@a}34;L7{)Jbg+Jxb419 z`rq{c%IC>?BS-b2Wyi_1wO{)vqte|1v-)7vuH4;Y;sIe-md2+2a3K1%=Bw~752&OT zuNj-hf%D|Up7|ysgw(ffe0ha~pLT@3EuKVZe(|A9Z~zBgImWxc?k2*mOWSz|ybU1A zI65?M#T{x6?qe(8z(Mo-rXHrd$Pb4yYSlIk2kh?F{dvB~&t{qGdZ8H(xI9Ah+UnfF zgg1mfZyOE<3SP_?B78a#r?bR7eURO;j@`2+LUB9CtYL>f6s;J9k!27e#+UzbWS~AU zm)$+xi}+vR%7{CG(}%H?iZiFy+<|B0EYh9gMu1bWY2Q*`2UI6P4CWGR=iEy3hPE-C^7q(sg)YFIJBs|M*{WdnE z3pK$LEH_;d{*PPgr^0oii|(HBueb>-bV0#XcpY#&S~BYwh1zFULW9e8H}qL~?3pzu z!oy;f_sc%oFz{)P%0h?;=~DT@48B^>%NO#^aorslFGnU1d1%6g>DjfDPu<~4Iy_2? zL^iJHY-}H1c8AW3QGwza>TvGn%B~uBcSs4nmoN2H4T4MTd@Y3Bp`naFwwqiHXfGXG zZ0jB^0*3W(Je|7=Mz_63t zN95sEuFHN-L01&~S?M#HjiH=tACcP76*iuBmuacVKybwiQ6p%a9V)xGMvyz}NaxvR{lnBVt;7LrX_EcS>w zL!I|!+iNsTKySkHFq_L6xSL9E>YZf+!A%kRPMg6EYK5!A>{)miPxQASk(M=L9zq{joNBACODycbIhe`)v#$YMTo~TV!gv0 zAymjnb8}G|4FKP(gzpD&oyIbCxuQ>o+@@TPO zs|rMx@(&w4cYu?9afNL4s?hHz(6&3-0rvTvBIglM1zUhNNww zR~Q8K^O}w|+QL1j=0bNDA#iZ`k=K}P3%3-%YD?V_0B)MTUUDm2c$ZUn=$J1*NLq6Z zGV|!MW*tvtU`_m;sSSIAd;Vc3K+l03kbxOm|!wa`l4-;U5rS@X}J6W_J zyO3e)fbLOcjyVp6$%B(e(%~920!V(A7@bO20JCEUr=kJ9A6S|%tXG5#^1h@jNdhz} zT)Hg#N(qiKIu=o@5+JNeYoF_qG7PB9vH8#=ef`|ib?PM)XE;NJ&Fmu{_C7Q@6}qAf zm$C>=bD4PHdHcz6d0ZKiKATrzo_NqU+>x1nQw2&Ns-pB6NZ-}i#T@!d6?&&_li$kV z;l;rK&&4-tAoBKAPxpHpxV}4d=JuQ#jDuw1#R?lxs>nX}wpSgbg;MzLqw`p3LXV=A zn+6ocTk^Fr+rUmA+QI-64WPs^k7p!WL%~I=p+o8#AeE=1`E|?+>Yuoh^9pG|_{|U7 zT#T(i#O>znYN!SXp6kq9Dn{||&vaAMqV-X@gaIq0u!Mby&UcrxH6XN7$Il*}3nuM( z#!b+ATIPO2Sb)g_y7zahsEBC+Bd63;W|U7L!EWZ3VZJ7uE$A+D95jP+MP+5DJDPwU znAPzon88mU-9zW>wZOwEA^rho3ip-VIlh0;0&bCObC-inA;9xlt@@BQtd;J%Fve{P zqNz31)x*1C?Q6`t)i4v#(rVSFpwfXh#gIE@)5hR%nz%j9eK+uv#j@=;GyzZF3&Cl7 zwV{~*rS6L~6R6v6+xFQ{3%qiER`QXXf>CvzTvm!E)TpeUHBdpmjdcq=TOu^zg~&m2 z=WbIdXk}eC=G1_&mep&Ko)Bv}tUBpvo30>N^zKE2lg7c5$p0;XB zu&{W-Kpmq3j7rZ+r_C%uXtR%!`YM{A<}$ysq5M3G0p^MayHMQh&O}ekL<@Lsl62(? zvV$Nz82$cO+5-5~s0BsPe8JnWJL_tRIV_(Spll^7K!9Mf040k#5X|4Fmu*pizTz{F zP27C(>T`9~!$ae)B(mIsVi<5_*D`Z zMnaX1dW_(bqHyr7DoNOaKoo<2 zMimY$ZY~;VPAWs+Lc#ot!#LQdrF-k!stORO1uszh<3LEQhZxVR2KG^pH%Do4fHS49 zZJ<+!x}wW9+4=e)CCR;9Hmwfz_+3xhkq-IAnKiDe1r1=c;kW-rj?N9om>NT7G{N(s zgev!SJt8NDgQN?~`cO4RSqr+*XEndaVVoe3S`gj?qj89m$cwZ zl493+YCSOfLKEkdpb6AC*F!szpSSI^Bf^TY8sKT^x4fLG2d`Lz89p6S2ZfJumhs5X zU9)i5Q-Y`lyIh{tZb<6G6|1$OwZp35;2)IP8->nYJp9_r$I*KE+Hv+wraqJgw(3c^ zD}#tY%pRH>`ml5`H{h*?66`%o&d{_S2Zdpee@NIWqWDloO2mUWI5!aZ_}dW$xV0dD znk)+kY1NPT=!o*b>a^>#>N3i2upsd6@{}xmu@5?s9B6>@;hC?^gvr2L{?QhZHwIw0 zXJ(RTSPF(1IUie_BHw1$YKC4HNjPO1XdqLK`cI(mF8Lec5V_+9M{1iPd~$R9%;ql& z;a^s#-y?tIeNRMJkn! z8bkfm3Qdu=C=5l%l;9&xz&|qI+$l;F4lDtIIP6*y|*EI!dn1Z>XP#5zv0qE6ndrCiJ0*VK83L5Y6LF3jD z7Q0dt7~p*(7hS{+AwMFB%%LVgNV{x$)Q=sk9Q$gixJ^KF%cA#}LS}Fgt6i(8D~dNS663t7;edb6<-kYNZ_a<2XU#XY)Sp&rX2iiXgWwwgLhKGI~wrQUpG zzyM@-DlsI;t~Ax<3W<8{Hvmo(bEk;8Ef`PX- zN@zXpn1nNCc-__5KT{3 z-#rRrwV=dRu_3t-unc+W_j8~Lq3m+O6GMQhY0UBXeNu<~vgWoN6o+y#fVo!@t&8i) zLyoiSfGszQ3hhCDV+jWa`F_yoVFLCyvVsPc;c<;0eY_@$o1N0tmv>wVNV?_#!rgX} z8L%5R2gZ*|D8jvf5_Vfpz{o~5A9eF9!Zo2?7l}E#Sme2ly}J7qph>#b zbLy%N*7^XyqDH2G@~sq{*^;swBlgW5mUWbebDrv9XU=M4SM)cFPSDH2HCNddy)rF~ zEmNgH(o7bn7q;rj0n!!$5_Te>g>w+wiE6}EhAua0>xM*!!H zH1yxQD~3~0!?bDM52uDn!^KSdpVW-17=QnXX)A_F$H^>8*+@b4Gm(2QXhkrN59ad*Zjx~HOB?U0 zB|(fXUC~PtlY~bthVm>o`LTYcYd-IvNWiJ3WR*cpUQALm{6p45aag}+9Kt2UiTS@t zs%cRXhXj0Rh+HKLcDn4dOTnNh>^C}a4!$v9){3vQ$dRwY;1h4>jBT`-T;k{3@5_V$ zH_h}c?~VBu}eZoNZX;CcO(Xx9K0*s%ujES=y4 zhv@yA`XQ^j&21;A`bRiGGl#`Xu>XVZj(VSKcwTnUn0zyLIdP;(@L+xT`3P3Hp2Hl> zD7MzL?T1Ok-X&%b3hCm~HQI^|*W)3clLZQxhxhm6>9P6Ksmle6EU@vV+>u6y1!LXQ zd;41g3o!LQc{|w1ft6hH-0hUg3~wfRzYCM|V#Qq{Ne>4YA>Z1voc1w4wtP8nsVjpC z%2=-FmOT>0y4XLsZl7U+OggJqSAvBx%2GT#J2MA_Kl%JP*FzLDe?g((c9IK3a= zw)t*cRTwf(xMuJ7k-_#EM^n`@i-J?qfL467ES6YgHN^Zy46aT@U+zpm@i2SYQ+6H| z2hl{+2fi=lF+TU!o}aanKtr^Yx@sne^+aSom)($rAXm}Nad^1SSJOfyQ*_K z(yM~gGVq{u$cGBibLgcy-aR1)7FHefu6puV*O=oQGe!k;r)_dbni}QPE{>5&R8WMz ztgCDNhvhMD$t~JDTc1dVRA~fBSt}L5N}}deeSF*W{cY;TjH(`BBLj~V$x+W3Z*CY zNx|rRA@MDf)ma*|qUcg&*{2B;JC@6_>yjAD2RXV+W?E1){)wXDrZ^UbzqKozMjIL* zmK_nRHg&5~oau@l(4BIW z*keG6IX*oYb2LQ{c#_sEE+>C&qB}gbr=m#@SdLAr)84+VJ3Oy2V)I=OCM`}7qUV=% zIoo`dcj)T_Tz_-EvXcrllh#l1+Uvtz4i=Tb7zPlU;ixVn6 z?v>M8wqFP4^>)A8C7=Xt z=fdhJ40XU}IQTQ2ureHH+n}4K(1lA0FYgSnDMJ+%=fGJOfc>nK-zayhz>2SA-aT)C z&UdS+8(FIGQRKCkm=A#31Dr7*{mmXtJ$6llT5J#|+eMO|>(RX+S^S{rt|D!yqL)`9Es7lL+5>A)&W$Bc`i zE_jNj&Fr66 z{q*4Z`C@xp-N%tE?$UV}mk1F`R9`p#m{gf@G2S++3rTAy{z|hX0 z|NTRBzq5MEExKM0b_HK%tvAz$r_AS)FCzL_cy<}Tm_DTWx-!^006d**_G;a#59?CI zk%@R+u(j>ro^;g*T7#=2Ag2TNEqnYF5WKZ~F_UYLHYh#Q%=59;hj+6b@^Oc>Aot;& zBOWpOVBn_&tb>~HApMd^85ItkQV;G}z%(HwVCKovWgJ|I#_l;>(*)(WHh2?B12|r@ zSN4jiHpIDf96X|A01TD%()Uz$!=}o(ly$BFNIYus({Ip$pl2!LAI}&-wD}ceW{#8%pHzIz`XG$qq52il1RNbI3gsT-J+r>roK{AWhHm)4) z({lq|ddKv^G$pt)Bgp_9!!sJpuHk^rcCVLaq5;wmMQ$5HHjd(A!j#w23_$ln#>IQ+ ze0A*FV)3;dhLB&WCX*Io04g!e^}R`kAno#|guWdIM=U#X8s8a$(Cl^<0TCQjWpgg^ zDI+`nSDC6_<)~dpTYEIhji5)RR*eJQTT!YOojp!t40lfN=Qo(t2TRy25==3M`+aZD zGzg+`m-+K;uaz-8bIPx%3D$#vD39B^6vptRn?0TteSh!c+ql)xxrMzX)heqCKzdKZ za8Rrf;FukEwHE*+a$Tw3!e|U{JVTjNTXa#J(D(*A8)JCQaOS6_kPg)9|8Pi_HwG^C zk%Y2}-4MIVP}@~x1ULA1nX}&3hPxAc`Y*~F0n1gdFMWGHtJDTN8?u9B9@h(= zRKUSPhzE!?C_`c9{d%pVIB1W}cTh&YH;LRnX*%d!Z8ss+F?&!6=BoWg@5LBE+NIf_ zyO@>WNwz)xcREAxQ9Mw4y+;9Z7b`m~k)9={i{{Lg^YZX+d{2__Lqo6(7JB49g3jg7 z4S3~6jKF1CKty1RG$>HtY3M`u#mf61-Ei&|2j1<>0XFeQFlibTwALm8>WL1GTq}lP z?pUrg(I5qjuj^V~#2bRA+z*R~Q8HlMI8!yc-4J+0GA$(0y?`05*2SxJ2GD!W?S8nc z0?JSBGo2Y`0J}d-DD7rbgr%3)cUXPJLFc#Q(S^C_o?tS<@;0pjFxoek^76_9V?4{) z@mL&mG}t?i)XRb~G3)!%CJvVR)~%Kcq~XptN3Wy+18|~`C@K9U0iv%vHGdm?Rs*Xt3hl3MnRd;Lta6y}=Nu z$KGdTWm3b_dU2nn1Vbnc_v?LKLIYuZdIepJ2EajBzGLFsw>FVe3PE7@$3RxY%D1 zTbH>T_0<#yIX#csRXGJQug~YszD9X2WA0SZO?U}lm#H17EsD5t@ZBQ zLu4Otma}!*$y*RRb4Md9KgAFlm2?I55Ab0kHUaCp1BOsvk0YB!dIH6$92tQQLs%s1 zD`psH!WvK4eInjLet(%xD>vk5u+lxLLBf@WU?D(nKd-ymv~B99vGaXHm^2t3vAq1X zX`uF=p$GEYbDtu=;5zfQN$o+n`Of!7P3sN+gtg{fl z^rMv@qB8#DYV*;uLe6YfqLw@8KSNkkl&)Gw=!i19kLGmXq6f4jC#Guc4HSTMKxQw zzcmAym+Kbi<%QtXc>V5?XJ%k}WTjIyUjQmTQ(ttjGY5va`#X%C_@Tq*hCgxG9F``7 zZ~J7T_yH<+D4!m%fD{Uqo)JtKu1f^>KR9IpU)bWJ2IGYx(!SE&V%Z$jxifj1Kce$f z?bQ_vUvr?isaW$4f-3?zAC^$@}*lx0?omOzk_iL@=bqqA;q}N*&J4EKLnY| zsDLD|zmQdzIgr`&FFio{A{UO{;BG{Dd-I)=X?cuPL8|j)){A}%s6ExV*^TtF#(9JY zP9%TsPpTkRq}ROsq|Qo|-4fb8qiWsC)L=TOQM&c2B|J-LDGhj~0ls6pKaJWg!T;s@ zQ5{xom~_tBE7WEQIvtE3&*W-@Fx%9~B+8%C#XBFC?yLz)GN1eAhe3&TQAQf@q0xbL|opxXQS;vkT>C8rq$-nk}sYbW!ZA%|vS; z1}f=E=Ba{+(p!hh*VZT=a>xCxZt8ISOzZY5f;RBnpfK~%s5-0_y9uTR+CXLh(@*5P zRN-?d)29V=Pjps)FsU(51q43wK7Y?+18;Q_13lj>f!vdH!^SRaSX1WKTmPT{=e_T3 z{Mc^|?0t-x%bp5wBJ26G!-y5kZbU{f7AnAnp#7>b)>e?9VwB86r3BY?ir#q+SVHM@ zTk*z71^9G=!))JvOZdTTt~?ef4^(cuc6EHQfPhtMlPjJ`56}CEQNY~-p7?zpu}qMG z$D{dGEn6&L;?83iwm}IfR?SuyePj-YC2aZ5Ka_w&+}m_oZ7rbf$V!d+4Kd`Pzm@T7 zj|E)t+jHQ|0a3U>HM8Ht(-MN%Ed%Jbihz{P@rcv2mLO4PXhob8gyP}P*WU$L!eAu} zhqE;we5bOvP^GqngyOyO0-umy-qw<^Tfvq<)7;*?g^3L!#NE$5Td@RbU$d`ACRiYV z%RB&o+6pQ@C0v}E=7mZM=Y_d?D;T7UtIQ*?09~v|cVB}QXs#}XQd06j`y1yw%FNcV zqOo_-;21x|b~EW@T|)MVcVzD3*M&e}v-QHab!+gqY?2I|7lq}sA5W8O+aSHCyLW@H zBv>6c+QJoK1N{w6hsxumK(jOd4)sGD$nq=Qozx-=OEP83DQh;czeMya1G=}ln{&tg zn<5@&blAQhMf>}scR$+>I^toG-KME_rz#Xmy*D3A#Y4J1%O$6S8gORcc73}Gco>XT zzi9YS3o;|W8`(XByHTb3cK*$+B zw7)yCy7tM@28M-sCN4hJ1K->4ZsSorKam}x9HU!tuuC>|iXFA*Hk(8Huh!|q%dR+~ z5kng&h55H(%mzR+gZCE2*#PzY{R!z>WJhM{VLOn3+E=uE=OZSh&%RMwL7ia(y}MX1 z$u1c}sp7Sr0>d^Ccxm}f@d+ak7LQo162iloqlZJkMi>FJ`LQ$Z1E~G&c65(-qJ3)C znPb)W5I*a{R*$PDP%?h$XU>WZR2|4!FPt(3;~f<-Ec|$&Fn;w;(#IUq+eU(VY4G4h zY`l_FZvmYLcr3=`P?qDmN(UDo^$B;Q{4B4J-e)7;na=>jZi#n>IwEHo7n)} zjjn|#D?C(@UG`LFLq38hM<^yC$!+N07wP?~Dg?Mx>{q0~fPAaYs~z|>i-+>8Q${ak z2=Hmq>ySzj%GW>^YmMh6fZ}=14{R=Y7{eQ%l_U^gNz%U$$BYNFr_|4?=m>D_{(N~X zlH<^mlZK`!&lYvUXlw=gzQ?q-t~3}EpgX$mU;@f##dE-tfu@!K(h){WQFAs>$Hsj| z>oEbiVpY!WRl$S&$2k_3XyiAhkoxp4J08>*H*&tL65uie&#;FzYM+d4M!P~$o(F9n zsv9TqpsD%MYO5#mW4;!v__`YpM|=<3IE@mZGw+lmhWy$4A6m4g1|oX)NyoWSdmOFJ zC>JV0rVEW#F@hWf_}UQP9OaMlrX_XccYQ*(Kg^QT(<82L9gxwcloJ^N&dcfA7N7Gxb0M1m%R3(j2jeEwlnEzNN?? zyD+Jm1{{Eb?;`n~cPJkPb#|MXwj(SaJtA|P+!jRMp@s`~0;YrqUJod3fk?Z^O|9e+@Tz=qm_aM=r2SmcyW ztBoVTX*)qt^-edi5WK}cf$CFOILFdkh5R&)vm-)G2~gFxsTW&JgvhVai(52N?1%$O zRBnkL5G;8v{kbav=I-;(nTn!#+5P7t29Q5$EVj@?!|RDsh@TR>KY<4!H|p?>eaIfv zI_QPUO+3sR-ex}D?Fo+HYwPtI4+6JKZ!Q0_pB;?1AP(c`LF^f!#ao8iTqI98*vAcvXjo`3tQ>kFNPZP9UeoLM!TTF9q*hA@hZX`J}vu zS$MFWEL}gO;0u24W5cJM@GxN!+>B-Vg4`$RDP~_hq|qt*^=#P#Bd)Uhp7`NmAAVo8 zmct$hG<-U%u8oJT1-sXSj(dQ${&D8nMC23ocE_nlkM=-w4cj2&Q9RrZB-5A|^8=%Z z$;b?ZPxkAHA=$tG`}?<@y%vJ_?Z09d^Sap)#rSNPH@ENGtB*a%3iVX%L$>cWyL?|n z>HYSr`DK6q>mj`(!5Uo&)ug?@2XrHY4rNg5V_uzx^Ee5@y)Pusv`o!$ICB#sJ%XIX&rwwI6JlRXH6ALb-JiAj!``%bA{O2FXZz zH_nq7ZMcErlq!ZwXIJQ>_=jshZwvXscyq^{?G(5_v{uup6U|@gV>yljRA+YkA-j!l zV-?Og%u^w%kqW^VUOdk%+=Ig^k1i_7?(l=A_XE`E-VRG}qGwrs|5rAxtCV;2hxfq$ z*VTDPMRBbG9Hj)q0~Mr5F-S*=N+;9{DBvzEGqV&~m(J1@rHK@25d@Tuh(wg85D_IH zMN#C7fPx{GAT2aO5d;+x8*P0z@4R>3IqrXR&diHQ8J@kUHYV@#UbobrMeGrh;S)&*jv$}R7# z@d8$}tVY~x>>-UP&sZ1X30eXsaWdEk)O5IK>IIJ{++WOT;htEI~Z2F`%|_Qd^)>QBI@7+SPD+^si$h&vEk!lOuly#J|Az0-aFQQxQops*_Uv8inC*ewk+(0vc6v3 zNWoM?Z~4(Y%n(d^Af1^&!A`*qRoBK{pgHknN3kIVReWmF-F_~xFFIsWcA12XL;Z&S zG5dfwa+s^Mgt;40?<~L6J40j6{;|`F6gcNzZk@a91ZlaVHWkto#3&!=-!tt9!fSc5 zvWzL{WC&?wzjJ__bh`|;4h5~M`BPepOfaRkd~PtIpww*MmCnnAF`oF!YE24ke~&)L zTu*^ma#ZHH3I$__zSNy)qhPwVPu%7c3HH`Sh=yKU1s~>$$evs< zyNJC!^Xo!pZWG|xW2rr8LtuX6_d98^*ayYKJ+ge-9&V{}H$UBrIjG|jp+@-(V5{po zH%So?vs~-eXpg<6|IM|h+|b?2t#%NlA16~zBjBu&l z)5^)E`xNMBW>nJ~Y{BW2BH#Tl6b!!m+7gO6XHVSe*?ik^eowr-q|#>tT@xxd3>}zY zj3m$4bXkGd(zYj6$C(hr(UI#D+6{KLcOg?neMT4 z3K;ph-il+|Q0P6llz;IUSB9)arH`^PJ&F`{jCf&mJeCsr5&e z(RL|R_g(3|f(#OThqm;wo+==n=bI;+F*il~&C zxogfONPWqXpyyg5QTFD|hYphv{dTKWoDrU1K0+oQ(IBB~g-bt+=Z{@YYrXd1=V4Y+ ze62LDOBtVXmU#V&04u`&D%TKm@#oeF>9rG(o-bFXNbW`Cs?x;uJOYG97iL^4XvpuU zorZ!j1YESs5*p;OMe)tNrGeO&Wtj3etNp7ja&^?&cCCYelk1vy%KI}=SW@8A{$T=~ zRbSRUeQA#r=6kd0;{-hDQ)A`h@oJFJrpZI|1oU)n3x5(wB854jMPC6D*c~G|WtbQ0 z_+I?gH3JgDPqs7oG4tVh&ET>+zF%^m;;)6qG7(2N|Mh?|eqWdb*e1;}(VIJRzn#H8 zuY?Uv4IdL6(Eel1ABo2VxbGD;w7uYfoLmlTj;0XsVma+-Aui4-ebeOe^#Xq0v#!Ue z%sZmyh#5iS1OjLy0W60|Cv@Y{jp4Rh0z^37k-oc}k^97}J=VVxaBM)w_2(gHgp1%s zTQdka@!PWc#@K!6GA$*F#WlJ#$_P|5%UJGG6KZ;#SG$gSGl8ZTN~mCsQ$~! z-`)sw==8JGI?oYMcxCC{_S2q7O8lpemKxmukxcr;JIwcTt5~z#NI*oDd#@G)`ww^V zY4YOj4w_ezd!iRQZXAB<`2qpK_V0^52E34otg_XjI0<7R_Z%M@c_XKxGPQgu5>!e8 zW0f<#(FTPbvPzaD*ob*bi@fzlygeUwzjY(wR`SNBI5Qt~h3jZ-n@Pe$wJKriSG-Ox zSYCc@Id1QlsQMoSIVjay{ccYu2`NHrJIwGn^Y@X4X4$WpBY8`8biEP>&GbL9ST@J+ zC->DQ&ki57vgw7?g)j<6`B#N&;(bu*p=rSrb@=~v;uvj!<%24F`Kahc3c8j9zdOkB zLFHW*=2tA3F#LVh2MIn7vU0Uq-|ox=9iIMfnQ9KwSGLG7WB-3_)6a~qGfHH_x8s?A E12GUg8UO$Q literal 0 HcmV?d00001 diff --git a/tests/global/global_latlon.shx b/tests/global/global_latlon.shx new file mode 100644 index 0000000000000000000000000000000000000000..2ecf8a6d6219192e310fa2da2bdb8104ac37ff0d GIT binary patch literal 588 zcmaLTy)Q#y6vy#j6{!wLhnAvVT5a{^wj>tQO$Wgs5{6FFsS$&SL5w#IO=A!-AO1{XB;ePHp&w0){w^G8dKYyP1r&sD|ojK)t%zM+hcd}e> zq@1$7``R__dt0$vCI2`j_54V^mejT+1uByFPznYmPf-eO;i=T$O^+I$5%)>qTjF77ig_Kyu5}I>IO3otvl4jo>X@c+58dlIoGp{72%zWb12GT$E z;+l@*oRsGOX!zio@gsdRW}fjHZTv>+vxTMy$a(g|jV%l5Z|x&{Sf<{Vax-o|cZ}T6 knR@d1=KOqfe%{Qd-(&zCq&{>h%@>Gmth)IH_G_>G0%NsEC;$Ke literal 0 HcmV?d00001 diff --git a/tests/global/global_stereo.cpg b/tests/global/global_stereo.cpg new file mode 100644 index 0000000..cd89cb9 --- /dev/null +++ b/tests/global/global_stereo.cpg @@ -0,0 +1 @@ +ISO-8859-1 \ No newline at end of file diff --git a/tests/global/global_stereo.dbf b/tests/global/global_stereo.dbf new file mode 100644 index 0000000000000000000000000000000000000000..5a8e79aec173e25b977b85926119fc49a99bf98e GIT binary patch literal 8622 zcmcgxL5f{R4E&H~2!X6ZAkY`+PfIPeB=FArkSiE$9K$dK>_s4F$RTo5*BfRQmS!2{ zFaFH$w|Ir5DwWjV{r&x~Z*I5SzqfzzT|ZymeSA55d;RnK%dh9Z4<8;se);3+$8Xp7 z&p$ptc>eS7>GA!Czn*?PJ^$(3$G68nKfj#3`uvCS>GQ+8$M@sE{;d1j%5jb#K8MZo z3nt~0A)R1cfrQ;XXhnYYkvS@AFl3zV{KIsJ`YN2Y#S8z>Jf!7PhnFA>b zGoZhbWzU1Ur&vMcWMiwtw?qbd*4%0Hc)-Qbpr* zrb$H(kSrsKmE$lA*?ZIETaT?HB4GR7tJQI0xj#YDHHeW8!Y2lOgc5*%$BRV4^>+&-PO#${Z>kL-~?i4)- z6pRPL%cNL{$vJVok%<)!O9L;nX%%X7Rp3tBi7o)`T~10`MwMB;Un{jE&>;EHfhUhl z*QS019q!zyd^pWET9uG6~i#=Bn5}l2Zm8J=#nU4JFnV$rUo|lBCm+lc8BFxo=D?lH)LU@xXA!eKliO z&0tOey`jOmYwx1K@P2(@T%mgxqs&Dzcjm4N+~{Y>&ovsLj17sZk;!BqohuOvLaq41 z5(4XpTp@F3c&*gFe7ImZ@W2YI4snwCH98d-OV`ENc7{4sI9sV%s%Uhw31HYl9O`T6 z-kD8BlLst%8=LqnK`S`x;!`c)qCBIM+8l_7K=R=9_+%na$%z|q#ch$0nd}3M>jM*5qN(UIhPTSTZfa3r7|~8JuFH3@(bD2) z%JkjZ4zQ#LiLU#Dr$tR>ljc|rY`5AX;|q6R>;k#L3N}~BEU{{36XF5eX$O7#wYjPA zYhcD~i%e{caaaYLKptkRpbC@^HODKh&Wx-qMug7=bL zyS8R{)V777bw3i=gj^Dgi#X5RnVJ;_a&38h2v$Hij@tWn~ksr_V|NChkYvkwQA^-6h$ytuPSDh8* zwypa9R0Z*<)3O^o&rO)i=yXZgsWfgTh`IA9u>CW(xu;)v`4SWv=es|=< zqJ;+bwN6QBFF3(m)uRAtMS=gd<5+a4Yf9h&A!~*S0A?&)WoUl~`UVbW;?d4Gh zm%jh#x-+KXKkGKQQS192i88nmckIFVM-9kn%c;w0+n!N!RdQMbzM9cA!DM|Im-guP zgom-D?4j%Sxo}6L!s`&+p;$k@Lq7)Jjj0~LATIz9SB3OnPewK$lxg}cC>8&6S~>(} zmUKSyT^NnxT0g>Cbc4BdIc*W{ABGpqfp%5arl{%u!hZ#JKt z{k5Rumrz=4IrePzSZ?IygL2F0c{uPV(}@Z=$>}_xvd5(&qHFctz%p zW1T)YnjYba%3k*TS&Er_G>~zEs(|7nukVFJ^1TCpD<^xv!!Cmk*A-Fheb1S{HW4wA zeJpow-IrV5oJd6%*m>LsTkXPuLyOh%YJ7fr1n9rgdNe1@iAyV3?AuCzmJ2*Ous1xV zB$C^gL8Iu}DHkJQG=0LOFo#dZ+`!YCw@tB2h(9ypWemLfqIgnm{Tev5uU$29Tl|0c zbm?E&mvZI-a(?5~Ez?!u3qBeL@#pnrON_K|`n4$L%gM2@t+{O2Ir;^cHgAb>u%R^gn!Ze&cgqV41Ki5H?;DTw-wT9B|q`ij-Y-M(qNgn`)lcJTnkDwA*$Wzt#L4 z3}cP+u9uEg$DRR%8+m=IXQO&v)E%9k!`h2C`K0QlrUuKhMYVq85!1~p~{_DP{0f(0@d;Ut@ zEKG;d{Sg)>lEU7vSFX?b6VRIias@T}c_!q*Vb#>`s`gaa;M|^~AeHx@`Q`j36n`ko zKQ0IMFvQNp=B9Gn9-5Ky_wvM4*w^`R+X2@c;Lv_0zQ*dsnOyj)ctamI|4X~)LAri% zV=AQ1J6O8wdM>wb4t>JoBvRP(YrXdC@BA6zsIE44wb67(Y&hp#bKptNKwgo_yG;tLiE8 zsRqUJkMkd&R{$nAhUZWHR?X&%Bd2uynr{Y(7jWwqmGj>^rI`&p(@NLqs?}c6qDJqZYk+s$8<}(0ns9cG-YKDsT6nX5 zd#H|niW>mZmiYjs$=M>Np$ z&M|n#nFsL0eIOZeKHJ#SC03^MMb%7Zx7Xj06rZh4X5ql4nK1^>g=QrcCk${`Vzu$wN_~ zaNmWm{hmLF%yKrm`7jzr^FvsbAcZZ*1Pv~<1qv6y$vG7&BX>&CLBlCI-l`C?FF)Qe zKc|>I-yAuXRF{{rd`%(zaZ5CEyIqQ5@5*d$l@`Ejt0~)id`pl+t8#Np|L}_vlz;O+ z)Yhl~-i6;>usOdNU+r&bn7;R4+GU^LegFQl5SWV7qK-S1aQo))GP*y)vm`0(^;mLW zM0jMC74{^a>=?cq0p}-1of-P`08FNSI3uq2-}MOjRDJP;{!J)<6)-TBCyHB_)3#@C zztN=Rtzk9`6_9P(+ zDkiMpS6OR{#ojSXWwu9yn6b>8;-Txv;pHHIRZ{k$Tac_C;i6v>1*7>PW0E3;z5jE< zYp2dkt^?uWAoH1Dn(_U*A5VKn>S5KgC{>l-rvL2!Ya=N;EOP1~VBM7n#SR)>!_WO0kZ})7i@`|o!5av8@?R0Ut9;~T^5(^CwBaY zFV1i4VlG9$Rqen@OMPuSTx#Lt{2vabsU0|)R$fqAldn}R%vY2teDQ|RK6ts0@&C{s zO^@(2os{cTUhth>|MbOhKDJ%C_PI>o9=5UXi(E8!m3&}h;RRG)Q013?oVqdWx%>sY z40mL_j>zZM-8)f7?DUa*Jnv|myCmxu9I*NLnvKG;F9#u{sWqMTwdIh_aDmb&BMhSD)S609pS9rS#?e4d@L^(`S5y& z6NvUm-j8d{$H&y&L9o9-@GR+)K@Wx62>=?2Z}TKA4Kj;QA+d z(LmT-o)Ot4n}Vm>I&|Mv_`{&U6CP=mWG-z^U;R7Z{*g2AhZDCLQ~KJHFr=R;4;Py?h?H%^C+*G zy$?C}ui?b*dV!)?WbE7W;^6vLE##A zrM3OquF$0)qL8g39ckhAB{lAqr0cHoAYh zG^DWCQ~BMSk#*nlU{dSAGSO!_cnJuT`2~>i^+M}T#Z2}%IO|#TadCL6S^;D#-<`jD zdp0@=Ts)p0nhz@{E3S{o$;Qg?wVnAB3%Inx*as{hq^=&EQY@5%ORZK2T=|+0YuyC; ztkkpdvq3(8!|Z&>G%;TN)i(<{ynJ^k-8pS;7Va=P|7YIWe7Lb_-x=N6S!gwD>c?+H z51hQETu3$`6Q5j-vnWsdmshVX4u6C`6o6e<Ez_;ai0#m`KA1?^B-}0JYjJw5!-J;2cy?9MUz z?)SHY3|mjD#?fu(4pB>{9UPxo2R=MsjVVf{Wmb0Wu*OaG@Nu6i?2u63^Vh1KOWVli zSYUx<2gnJ{OY=9X;`aUW;w>ScRn<7UKf?2DQts0Aw6trQmq7NE@Pwn~&C2!Qva0_* zDn?r$PyNBiWwQ_BR1DqHTM$)H56-KWclpjwLG9aa%4N zOPcz!%l&SdsGGTD+@AgQz~NPh>Sk~2VXJEu>&dH3tdF>SKE|sa4CZzx$xO(=8Qa|B zU3eN`G{0mBb4X#YSMY0#!I|Uhp^RvT6XQGZ=2cfIrd|VF_|WzM0@~Q;nX_J;y3gm$ z+M{^B1KG=VxgIt?sY$d7YsYT3-dqnO70xg?x&8mpQWF2UJW#X&n)aPI*r(Nw?dBN6vxg_`I_$?;O}& zq?XMula75$#?~9n$bksU!MX*EH1sn{I__zl4b2jdeIg&H;;Ze!@3*nCfLN1xnO{&J3rTgAV*^S9qWRfd-a3CAb&iM4j^+yUW>pS@%+b#zRjxP01+>(l` zuM0JQ+?EUFAn|Q;RVvyJw{Cl*o(t2|COw#NC>_sfTE@pU*hPi1ItG;zV5PGJ zkF|*lckEwI37*^X%^eCYes8P%o&ZD27K#V^z2UOuP-w1g0)&;k89LwZ3qNv0MJ{7J zjGiAdW%Ec0rt_>4?fPAlVUNBBA+I{(++k+P%~032y_ooP>77CMNcKEaqn3WFoNGkG zcs}6vGb;}+U*rKeKQL_-;|sSfhgJb~fkb=5m9h5yqIPy@u0a(fKe#MCiAQIwEyfHp zLatT@0N5{Fl&$_5FZ#TATb3IH-;NqTzftxXMQ*siExi~Di~XJ#+#CN6g<}_0B&`pJ z4HfC01$7^B!5Yiy6K+R>g=TU^|J3JPUiF9LPcWsA;d$AXpZSL)z&?4?B17J1D44p# zVqaDmM3gw+J#pYIx`j^6|EL`TJmIQ03tUWaj3TqjIx7$YYfR2)X1%}-TS`~lKJ3q( z3p0nyJ(nN(gOHrvMTKwsQ6X@<+PO9#=+S$-d0Y4k%!!W?RW|a07wNfr;`iTT$_{Vh zRr7$+^F^j`J}C@3@8x8t?t3;22T~vTeO!`;o7XCzyOuiwyH6@txF#pER6>P z<5eF_pepmVA{_hdY4T7W-{iRWX}zkQGH zkMO*J6t;f1;%@fg_YLK6>*P60*9s=ywR}6EnN|gX?79EPf5p-7_rTr%t^F6j=OHqBx_H%S-$@!9}&Fs^E&0 z>Kx}S4P0JVy4RSdIW^$U$l`{CTUGG2TJV9f6cgiRsiUb1*vITA6G!t;=3pTy?EB?t zGjo^PuR_pSVw0iyq8?x6NN&n7EQ41&&4|Ta%U(au{c=d)@`}4sWuQgH;ktTU5I)95 zRkRe2Mvb{vzN;QrdUjg3%9L-fBu-o4;2t zhP;$#r<$3yXudSH@ci9kxR~&*sa3lUB{lx;&9naB-*t7ED7NwION$bi-|o}nC|r-? z*E73aOiSTy{pUlQYwEaT|3|;6$I;^@bGC?-BzpeZHjYiQ%ql?Zt9%!pO?iyxt&Lyb zZOO)6BjyL3t?b!yk2;s=cfdBE&oHe3Ew;)#>s&ShPFweHu)4d>9EdSCeoJNnE=!o+ z!*jz9cHH`W^x4@0%$nOG9yaF_81&$Yw-NMjyLUA#V$uu|%wA25akBfd%y;nC0N{n;%RXrhEjG^`b|=kITh% zd*`Y~PkIGKxAoM{KFLM>cY*5ywtR%}ni|)#xLh2$;dA6(@Kbn5)e7A_>}(HOB0k;# z%B;t|xG9p4uB&Btq(9ySW&!B8fyfHZTviD&_a`nYz`6M|;wzTlL7NX7+vkw^n|VX! zMCiq)$kFpjrg$+a?0cVXkrk?mSVdj8HFY4+_`&b>b%IA7_F zdpzn=RWPO=_Ds)jVid%3>*|~h6?KCcGzqyf?zcldRGksC4vC6FksZ{DQ4cezi60k% zy*rM+SbL}*f_-)#&=U#8+0)85wrbSF(xc+l^PK`w=6TEB;92#scH1FiX%jyrLUn@D z>3S}2t51e*pYWs}hPDUnPCVy}i@*QQNJenS`eu)=z^qR&4bdx?2%f4&|#b9c`%M02i3J&c|oGRI3usiOP8Th*)~X-5G>Y&7ic z+*6Em`10YEb}{%~7Hzn*EQc*OsMCyo0X(aBSQ;0R`+16kvqTXl@`l;0uP=n(-!`07 ze_4p{YG&NNb)pbld;xdt+e=JU8=e-z>*dDN`gL;ffWBs{<4_??qG1ogT%4b= zcU$7YB8Z%Er*cJjF3L|xS#48R1iCR1EAeb3#%fj-G$A zBuhzgpzG)tntwCLFc&oL=#Htm-HA8F0#*lYssI6xsQ39EEcQCm;1v7&J~YYuP&ugF z|6Aca!a}`Z;`r_0x+zbBxpCtLSvi`KmD<6YBKWZ6G+$6V z3)QFzu%HB_Zx&AJie_QhJmu4eElRm#o&6?x`ss!;AP&oM3I#0gTpWEo+wJ%e3wcJ3 zpUjB@4-YR}50W)}W;6R>C&`=6(%N8bjjF2bG<^&T>VW3pUkCIiBpD3 z{xm@%RCGMGVzGW;9J_u3$Q^PloBY!OAB2XU^VyXE#|oDfiIjdur!y;>j~OPw*QJkl z?ELZyEgvrRI-!@qor|j$d(TT5#KTaR#=S^yGkiZmInZu)JRH^gW*d0wPuFHh1iGyT`rHs4Fg+;}56kT zhtd0uY{6xuoTlf$l){I6{zBl<-jI@%9-Mm1aQnHQVQ?~>WI4m2i+%n%{@Q=^j$v-y zQ%h9^qV5iXoJ^1IA)_AT(2DcT`98I$8%bcqi~qt97*qeyv2ILS8Rw(mImD%%Wulf_ z{CXG!H~BsqSJuVtyY}xc62<7on`H}Uy5lg&Qae1N2S@i$*k4Y{03CFtHgi0{q zvdrnowRW^`e*dNVMGZ_}dtu9^j(_4s)LBoz=80!7jGtKz6aD!%A5v<^if5a~%!sK1 zWeR0*Xv2NmZ6D3wR0V#8VNkKS4Fh@YWwws30uqgpwJ>SLzweC)Ralj<4UJ~5Wws!v z@0~Xe78|Iw;Fb5=w`P#`f1W;OR>U^rgd951Dqxc7V45`J!r&u&pI)kl0s8yqW_&}$ z+T8AHF0c8AT9<{nRD;T_RePi#He>wQ>gRo)U*@|~|2oZ#2 z4Fq1Gnq@0@u6w`sM0T!e!^o+JKN%HP0}YV{X|&!}f^m>scTtVdSSi0`j9We#E^#Sno z+OkU~5}jzguPwJzDi=<;eE6JF+REO))TvLu=5NQgOtr{@{<<@c1wAbER*#&b{yGC9 zcm^Bu@>wW3Yw+k2wKOQ(%70h*I}0bLdm1jDk^(jKqQAky@|pW4yuX+P{KYNhGRs&v zMd#xB$#DsgwtuGnsb=leLo1~+54ppIb-FF zZHn5Gvqfh@!G7AZayh{^e7pE5ffyo~(6BoX`CFMr5~87tcuChAwBde+M$Gf|u@K)- zTe8Qe4GZ}UHNLNo|DPYV9hHv!625Dl2uW;AOgm<$p6Q<%n+zL?_d3zA9d+6Mxm2L= zlDy1zG|1KY`oc9Gdb^$Hh;Qw{BX>_OQf$ZsHEkLM>%fEa|L%j1EGVPCiZLv_r4zPE zW=RfcKJ1&TCdWc4IKU9x9I2e!?wHhjom3}Y!o zZ)XRZljweZG>PZ`uP^zZ8Y_gN3m{0e&f%GK2SycCxcth>hbvDg2)G?rcC1y~)|Cef zrp~DnUDS>oUWFA7<(;-^!?T7xF9y}~;YbP93K+!PCyfWqYih)67F27sjPq2xxW;?PR|wl ze76bxAtQ-lrxjTWd~;|ZqX|!czH?b2s1$7v!mx!~(RADj%L8rY zu>Z?i;%#b04;n9eP!1N8?kl>fwxPv5UZR4Q!y$Qx*WT4_xWp`FS6)sz9HpUpvvv$i z@{&2rUk={nmI<9s&gH1}CR<@g@$z)zek#;=MRFaVNwpr4jhr!{!dH#yA?k9K{7ua#GL_D)YE<3ex1fKQ= z>jh{g;%NTK@K=*4JR660U$P^y%`yp>&k-tBwfP7`SA4f_*y4zm;2dyCH;Jv!ZX*ro zceBVSkyVt4Q;B*n8IlM!x#x}xeoerg1PSc1D+#7vG*Dl`i08JQaNhOJ&z@M^Z2xX$ z4PO%MG&)-HGb9Fok=gUsPlS&=YXWpVqEU@r57r4#xWKny&*f;GW`B5nuX;Qz+vFp@ zdTA6kq)!X?N{NSdU#p@AC!(-ZT%v5!-8kskeck7qU=#|4_6hXF#DjvY|Ex%jNSq=& zUPmK34u+pxX>WWOf!Cy!)?_`1g}l!cP!x$pu_qUtUKaz>%OqYT5&YN9j%htVf}&yh z{I`SeG(zy(+(%#H-$p~M!rq8A-eLHfAO;JDqTq4m3n{OAA(&z9Pa*w6G9NyzYxnLo}FaY;*8iiB%FM3d&3`Qp#uT7q~D zhe>R?;*Fliy_z;|41=K|R(Iu7H>`^MTI@;u&(-8}gdx%-qhQLEMNVPUslN)bt+l=dSe`nS$v7F~0D2MC(tNmLqQe81dulGH)n*_;vZ2)jyE$ z>lK2mcLxS<$fcr>u0%iGyQhBAPxvZyj95%R@b7{%ye>K4VYO@AZIi{naBkoe=Xr;| z!8U!W*twx(V%DLQLK_evQQom;Uby;&_w6yM4ls3}O?vHeU)=1mro&>`7ACUo5I^o- zv^)E-c*8w^R5&JQZt89e(-Xucx77!t@|i%-8aoH#M;X85(ZUeC{OH$o*3BsJ#H#HHWh>FHnA0aL@#!0||M-)`WA&= zL$+<%(_LXqN@avoLKJ@K6Fg<5;Q=1OBT2U!qp<#A2(cu5VD=rV_{ZS3@2kHgp792K zlI`%|R4lqubGgkEzT6ogrzZ~e_x403p7sIPW7e#lLosM(`FHnlyASNuf4NJ$AsUaL ze)syLjvp+Qzi1IFAB+FsI(#9AL~)jB#-drCA@O7cKp;WqXYM|9H+jRR-%li_ERV+us@vc2gDb=Z9?PtFfdFw2kZk&1$8FaKy2 z*~H^e@jHStiw5d3IAR))%H*2J4UPe;zY?}f0uu1jgZac05eu4h*zpNym_3Ui`{LjZ zTec)(<8FEp#REsy2=g2)3yx2KF8;&!_|7EZC*9l1^JJ4?w0t2VXALQA`^Al(rHMT# zIQjb6gNp?NE;+Npm z<+Oc4LBM`O43Av(o%O^A@+ZxlXN|ZJgLoZw@=ayoa;oqg6lC1uql$AUx9_R!ctjW4 z6PSimfFOgz+q7XD_r1G?7^C|qtSJ(Gg@@id`|lAn7H=1RQoXg&XlEDhRC>MfbbS-v zz5ip$<{9nmdxwH-=*Oum`d4q>xY@4&!po(jI#?WptWoK{&vJJQqSJF^`%DNeSgjmzt<*K+rz)!Hyk zyrOhBzont2U|2=e{V(kNx}inAuOC;v?#jfdK(%`&0b$JK-RaG017s>A3?lObNpeI`POSqr+oHI?(^*g!(OEooGQK%Bfdbs7})LJcYV&>3fHP z6%8!z+6@wP{}laBEEXx1IIq}+f_wXU)yet@{O{k<>qnMUi4=A`hBsl=#l};-3={Ql z#7k6zs!!}-m4I_vkDXJsf%&BV9-jE;ebT@UY6p|M}O zY12Y}F0EZOK9#`7=v0s-$eb#CZkPTm^EMwt^m%4MbaXk|({-INh;oY+ZZym#Vxx7kJ0Q_7& z#`smvr%nB`aFIay$s8*_hHkR5&&xYGXh8RNI3Ht*`whjJj|*`1RRZK(#m`vJme)o2 zk3G)M7{g<>!O*jmI~Sbsk~v#T&JTOulI{7cst#QV`lHUF4dw}x#UpXDlgVZzxvdSp zuyH{3Xh?R`@$GHk)tR&F?E40!Fw`8Oc34gA0{sSjPGLYl+Fnp586sN6+WDM$BDKHvG62(6~<=rG88#)_(rY zea&#Sh8he^Oy2#=NX4!R>hqNawY`{lU3=E=r@T#2vi2 z7aCDwkSZ2T=(_VP_ladArfCXpOJ2<6&b2IM%(wwOD6Xfi5s&>|lKC{O0ougwncv#g zh*xb4zJ8nA05*?P2KucVk+W7v@##-VJ_U7^ey@k&=%||=)l9rbqAeFzH9*7Rn~zSL zG5?9T6gGe)dmv0?n+bcEU_@c>%}i{$o5nc*mkB=9mp#ZtjZ*Z@k!u8^+;p4mXW|kH zZ{6AmS(IF#-J50j<+o=`caEp*na3XkYvwmEwa@CvOjXEU^`ZeS^MZ^&<=8JL%J3xNgwLm z+RzSmkLJg`aH+)#wtu1n*x_LFT0BMJ6bo6PMg7zpYSDO{n7iK906sLBKw}+ zJ6d~)@ZMPq&$M~Pf^QEZg(#VJ48l$tiW2-M4@Qq38uD4KCLo$$gW&y5Fn4))7jAM| z_sX{E55a>HkdkaSuBHAe`~&tENMcGa@%If)RC_%98_vqm#F-wX$0}&wZ?LBOptuL6 zH!Ej*J{bUZ%1d+)iqm0c_k-$@5$5t!J-D9ji0=o+2jcij>p|I2hsg70eXt*>;n{-| zX`GAT1>!L?__}N_8p9eAN$Z3A6pEtTgT-Trz)k9fQ#n!|PR-polTE)DOqPyOUf|e+ zx_fn=YUuaENG9<|F?!K_|1CL7?188B+7R!>{Zhxr>WlV(I8Uwg>Yg5Ki&CDu!>}6; z5$i5scQ1D}xW`18UR$i5~Y3b>Ny&uM(6Mo!>g_0DQ z+6k9m;`}R9`|(qSQe@wRPS6`m0M(EBG4#W|x=9HvcuixKa|h5=+STNtJ?o!5nEU}e zLh}!2v0x$1J39IsPg3)6pdGGWS|iu)^BYrXxTm5GPHpum%USg2AAjoVHkdwtzkpND zAC&r(H}2w`RtTrLtVag1{W(2&&D^!GY-1m95S-@{}pbmpiOwjKk7Y_VJD+0$}*3^)r&s2iKkP&4%AN9XU&=3i)z*)E0dGz;4T}Q z)PwIyMufIqJzYseylBPJ+GvNbqCtug>#9DS3TvS;n-2|H1A5IK!W}y;=<99Z}6yu-vN(WgO7u1w#EY}QAGp@8w$?rh+vdwRX z`sJriJ{oGrhXR)fGNK87Z?TipTGfuv9%T6)z1<88+4nN>(~hSI>tvlZL) zHqGXn*aEriecOr&8$PSvif)1UVPC_p<*n#%*hM5;E3oLc32DJI3sXfCRodX=Y_r=p zzO-P$pE-_8FSJ1m4T$TuqJA~m3EpiW!p5Gq;W##~vmMy$EYXHZMD?*h-2vHCRlIys zTT%J%8iKNA!EE9M@_yHfvNQ+Ej0N7bJK~3C`_t-Jki1p)WT96Jp0XtuY$ywKD#fEm z3Yt-ZCNAi+;4xJ%7c^tsyA=jYXR}~nmV$Eo$0pP#x7$6t4hUjnL!0pNB_dJhcfbdB zoVO9_cD^*B16DpK6Kc*x=scIEmD&Na^ofNkCi+%4ee%^{LF$p!wzd-*aj&Mo+k0OY zIGfHTS(}ZxDB<^)J!3ne$!R}H(PX0DhVg|?*ED--FcM_J}TE+v(|zxQU>FaY;9vC)pRvwHr}$fFPbeb^|%> z#DPw5NyKYEU91OwH|BlH>|$b7?7T!<^B!16-~$ij8}Sd-OS^ht9*v3UF>zo4am?TD z1r?eu>c&KE$AZ-z1HJH@EE{tl6BCFfrS`awdl!79wux{*e5CU83=@M6RG3}2>4#q# z)Yu~THkCsw2H?<|jd?~qjX0j>iMS2GTKZ>BBbwHA$>c2i4Ut5e%(rgB_iTIMH;e>O z&{Z?0@DT;9{Wla6FC(M41sBrm^a#06NH*Jh{x&3F1%mhb1Fw!Qd!aw49sL#_zH7O2 z5VovcO3+*#+`H=%J04B^$tymJwKWeym`dWK)?OC&l@>m}nJ@^8DQrHe6BS7Kejyuw4sS7a9(fer;6U!%S z3ZL{J!)PMW!r4Piq^r1IY7FDIWjTqAHele&YRf;?{0yok{(e`7PbG~6E^Xy!{Nq1v zz!TK3(#prMlxI zyX5LnV<`Nm`aE8S2-P=qYO%+RT!=wD44X|9Fkg$`(*&)xiU0Cp5QAvOHMpO~$zn#p zn0RB3ET}}$EZeXZe6Vz0E)`d6|Y8BY7*=k02(#ztf>0O--`WE zdNizU)6FXM3Z#yuesEy>YpQU^nVX{UPWn@djYPHbis%6i0r7*luLK1o z7hhg(+5Tx6T6U%J8hV*Sa}1=(0`7+e>;zZ*0*vIhD*v(WcBiKI^G z235P3!RzlcaUt6t?t<%o_gsogOUFetPIs{j;%xZG-0VrkbPB_r(FLc%#zkHtc^6?> zB=S_y3D2L;*8g)q5yOagc9Tmd2o5hi_^mY_H{JCoPhWJx=>^tbOQZ>|Xf9D;1v=q^ zM)b4`xe0iH>Os*gh_;M*s#BPVYl+GB_c{xX=9`nu?|7tU$m0tvP@woHmpF7(3+x`6 z$^!o@n~t3${;0l!s6*pEcfhYb-bcy!_Ci(` zY;1+YGlWhqHAq4Wy}*za5iQWZQ*oBAPXeyCdQm<1YzssfurWhsz+CxnCXo)gQLla z+}GY_wJ&$Lp%)u})(mbf6>s9k{2(t@41Sp7OZ0=jgkz?P4!Dvo(3U3Pv1aCs^L>ZP z6uY~p36#tx&hll=#(sK$=QV+FQ6706;~h-8V7x|XY!g@{ytvjUWCw#^NJ4XLBP&XSh0E-yUt|$!)1xU-v3|gz$2&pwqji( zP&ZM`vN|aB^UT(FNCE=RA(_UtAnivTL&-ql$1z80!9M!x)#`gGAfZQgG=D8@e#1Lu zmqQA?Z1#BUVO|4`pX(lJ?MelInT2y=7u7)2fgejfJ5xb|hy#(>YN#=wk7}d=i)3&6 z9ErD^o~FSAw!K>ge>YPFISr_qYVfcM=IZp_GwDi)7wr6>DsW-ra5CT&y%=n& zU>Q69kpVhWQO;D9YZ1|>2T1Kr675u8UkqAQ)D`@Vc1pGAp6CwL3vI(1TRg| zeZRd1WU_L^b?2tT$;rgw7*+%K_J$FJN*dgGO$M~I7K+NUNA#Db!?01pbAHEK5FlQN zH!>N(R(Cekg5E;ewug^0K#mPqu7SMyFL$ea$N+z;KOC!p)imA7I}=_TE-`B4tAQ`~ zC7yMuW&ttg6AJCBVTQ7G@Ilc`xViNs@xWKZ&_XkQug{sVgJRZttKji^YM5uj#`7ei z`L+s-esotK<;#M9@&>D*L`*I)OehN`hY$tpLM5o{dS5puaq5d=^8M>RR=~|?YHD9K zv%sZg&#vtS<&Z=3EDEx~%#$vBIb6~tKr<3|4m)WxC+boeXnv)@qHL%lD505OO9A@p z)9wglL&PkNH}5h_U>A*V*JZ(M*;AVYgi4@y{);;8nk-<${5KTCWhQwLBrFRqJQDUj zf3FC7guUeMDP%#G?dg%*?S;@Ky+?0rdM5a#2C2F$6+#8gYfZ@jJ2oDp0A^A6Qb9WG zHR)Dz6fFS3-Hqghgfv+5BmB@NgM1L9YTm{)SoPUkF(D%lL?$b_nDeHBOhmP+{*^qK z!}cqrLQV^Hz~_UtZ|X6W~V1%4hnww*>Pu!$bC=b4c8T;}WUnJMt`YSqNNg&A<`ovx=1GZ|tSqOzM4 z)1Z>-SiQ-Rvm}eW7?KLt<#%kFb2S;XhqIL~O-lha0@sjRkWBP`%{eic1jkpcOwxLm z0)g}Nqow~Qfy4&ZiX%5tAa(Y1b;&P@pgH?+^HI?h*f!&ytgm_kytZXc`T8dXTH__# zIyc2b1%c^j`H}guVN8B26$d3Wohl|3eCAFpW`xGTPMS_Rft<$+{s*^fN5ccQpDZ05 zDokv=h@3iqW8_0qRyv$0yr`VDISLNHbaK08kOpKE9NIV`3gjivot~MP4zu0nlxZJ8 zSnG1sp?r2Gtl^zHcV|N+thnJOE2xqM3y3>P#S3A#V`#`gX%;2mgHKlhr=)%xSRt&Kb`zyO5Wn(=&JN}mDsKzs3LCSD>rh$ zpQZti4~47bWeMqQLQg^J#LaDC@ZtXMdVl$B*!PqgMByMl-$duUVHVuRY_sXcVQ~2A z@*qQ!llIR$F~Q)=j)!FdPe1=GDMlbjtK8Z1`brkCueII&aASpF!HvI}(5oK4-gdDs zoU&JL*>9E!N{Nd_H!So5Wx0SoE0r_g%)P@KPS5g&lsUV?`PZidBO}J{&VUcRXZym_ zAb8G;nF-(h;NU!(gr5Tbm+j_uY577Iv-HW2&{Q}hJw3Eo#v63}6NQ^tsj$Sd=9h2B zFYu5RB$41$Fdo|epnZ%B#1q8Oypd!e$ys*hH-Etj^2&rtV=}BEr|xjF7o3VyB`c8( zlG~qL7nte`hlN{Le8@?HX7gWd^9MX(D#>1ayD|yj3r*zs1&edumJBUQ1deTRwZ7Fh z`%nTrGRYg{aL^U}_i(6Xwx}nfzhwz#8*0)iE$lZ*hd;^dJ!XBywfZ+Zb4Je+}c3YY3$9 zNNT5cMZ^6hI&tCf(A@$K6}B`6@HoLJ|=?}QOly^J>jU3 zkfGS-BskZ8Y?ErE2UHQ1y<>L*9G>lHaO>(%h^5!=&qR1UY^(m(><5rVB7v+VIQC7( zx9z9{2qZ2soKlnowT*|Lxd(rSKq{8bCBbiXeXl1O<}ho6v0koyGE5RnQVjiQ0{cYQ zM@C#pf<;4b-isf81qumKDrPc?(E4{?ck#}T@LQt5)L>fzxY+vS$p~7&{n_+|k~qlv zU8kR}@EGRN;m(bR6Q_U3XDb=QTK0ZR0O=(1Y)H*ncroyD=?0rbP>r(gW-2~FrF$R7 z4XjB5ofo?CPQPE_Ti$1r<{e9hpkH72YECl6C}x!Yvk6Jibwa91t<4OltvNMMCM*H? zV+@akoUumvw!9P9mJ__LRB)HNKFO;qA5*#7C<$y7T4nxxw?*yg!Th&Vl3|;IsOG_z z&uBG1nE5+51u_I*tu~6fGIx(42AY zm4*o@8drak%glsLZ3_qj-U32n=dZl^I1AbzT)3R1U}pAX9Q*{0s)-5;9P!M8Xc%=V2VH?$X$q?R8S?&lF4|Hyg~DQ`4gGxup(d@j5e zi2Qb(ynmX@R>C9jMEb#8R_)EUq4Qn3w^w$VyL2naWG5-4#71+3_4Dg;tAL8;t z8S zSX~IhjDg7D1M{+AeQM~#CFC71HcY^o^k-(0-!Kw}R4cz9l?6V*YZoc)jX;^wl~qef zePM#m;1`NSHk_o7$bGiFjl!01jgQ|&XTqyn;35mrXiRZVJF-A>;m7A*UgX`u5VB#q zv%r*nU&kQ(T2RV{2@I!OO{d}z$%}3>JF;OSO^^B+kEZk7O?@n~K~PWKiBCQrjVzou zX0Oiy>4&96BaFwJujfwM@;MtCvJdQTkV(WKvbwQvvcYJh@Hi>0B=jS|Y0xlOmOT@7>O=(^({CL!MI*#NVIn`d)K3tLu6G(Q`s?`bT8R9fj ze~hg677BEV!#E04an1wdXcy&sA7b%qFa;Xs!N?Yhd&@%--tnqPB^)+0QQrtqsgyhF@!v?_|hXErnd!G8|22}2C9Dv6u=ozHQ}tl*njZW z56L?f)IYNz9tTCIFZuet0PfNBu4xHKU$MDK-Z1q^d7grw6Okk(th%zW5OlSi>KT)g zahSqvb{0Y-8|a>dJK1l&klaH4rn9*_67e5Aa}o5t=5OpDibwKN7_q2}V4QHf(E8;G zsAgooPIpcbD6#dS1nju^z<=$z!vDR)K>Q08`n~7GNFI@A`w?A>INZPW1kbb||6f;U z85i~1eSJi*3lJ3S_Ndr_d5wwP9oVfXc8i_Z*qwlaA}z6o?(QCjE(H_&+55u#dGSB* z_?+W$nEA~Wd#&|dDI&zK&LG1%{#Y$@)TzQN>g2C}$NaEkhxL00r&KZ1+SNU1tsm5G zMXPLgris15Joxg%!kZ^{P3@O1ESp^F_2{TCrd~>U-}+&?=$x)>Wiv~Od85^>C$-5C zYveylOkmqHW2hnd1Akm2?F1!9&g;P0bQ}BrnDbU3s1!nkgK7 zj)YpZ_r~e$DQ;z#vIM@ep#=DH=S=aa9@*lB^PlS)zkqa}|9#FD+w+dlPVR%}-oDK` z&r*q5M!XSvfmC90La0Q|dc1FYfFLJyzEO$O_Nif${p`v<-}mKi z_g2nuk^T|74@+ab6Rxl)H&ZR@ShSwHY=|Q&m~VWB{{Bw(RVV#kLdLYWsl}i0PeX0m zJVW(;S}kTMYQLD~br&q-22E0nm7gcAf4t)Y#?`qv@WOnxxH*2-+<|SEVUzVT3zO|? zQQ-Hk@_mRUzLOlkEI};-vdxcqEZ-r9vB2%47Bw9;Y6M-lEaoyM?UGvbOt^UC@Q4Ru z1i4~XS*gW!Z9Q_-J{Q$K8nw`rAcR+nfkIof6q>kb&rbzugS4=(>K-Mgpf5>}Rp3pe6loIQ%f>b~U} zbacSV7hNCMJA87Go9MygAmV0>U+H|ouCuGCL$bJ2NxsNAJ$m4!p3Y(b^JQj|uD*>! zTW^IvUGGl@6bPM9m;P&W)IrenyW7)W1>)V_Icr^%&;Q;33q|b2*iO6q-4_x`Tv;HB z_MUCo-u!|%dw~b^1;T%%WrJ2j_lS$XNRTH+#`yE^HjaL}K_t-lv|&}g=>Iu6X4}}4 z=+1cAF?B6fQ7iMi%E$ChO0AY-Zn^GV?6i>gA zRB<04rO&-U1YVtW{ebCBB=Wr5q)^-@cc@;}bxhJsNOx?0?$NYJTq^3}R44E)HjWN7Ijl!NTi!fgf3GWy zDdJqqv`Fmtn07&Pj2lwea1mT2iX;!o6{Y9cFjOpFY+)d-D}3qknH^Uw_P4mx%GSpX z?|z;7oIa0y0AbH6>%V@78f7Hx_>_o<{!FO#fwb#Zm53!A$FcH7@v`3!{C|~*E%H9# zkAr)8B2gl4l|6V-<5(bQnoRSD60wGJhL?h|vdR3x&LJgY$lCgK>`#TlTIx@wV#>wm z`}+P3h3y71kXx3BKk@Xo_=JH7k+f15i<(5oZL=aA)qd+@u|@(mB0-@n@*foo7uzER z+uujQjz|r9IYr`h-`Rn)_r*Ym#_!);i$uZOSRzHmpp;J2Mnj85N4}R>#o{g77VZ{{ zJtpra?oN)u**7G3^(_|Gt*3lhbto2Wm--PWkKy)q6fB5NESyst=T-;VZ3-x-gq zPH$ecI#eW(b2n@KohM#A;XBH`O7 z{Jz$}WaQH7c2-p=9NL-m%K)YLR_S%d>am>?cSNe%e1k=+|PltRV@s=awUSvUs zjU}Cn#Egk`<|c+`BP?V=N_?Xt@k;BQ|E;AebUF4UFVUk=wBxZI=?A2SZ;Wawti;g<>%! zlbhbl#T?EB)9>4I4%4#oAaTRPibT2OMdzV_$j>_2MPl^!d*c>%%R|&oi%}n)~*$P_!X> z^{7JduOHU8P)IqhNfFqF_Vr1j&}E%VDGL7O5Q@-OE>H!ci=A67*NMf|^Xk3@ zV!bUxql-bIeqs|)E}M1E*(H!U|F!vo-?hbmO3=RAUnSy5uAXLBin(5FpvV`FS4k;gcI`|tnx zr1C}eKJp&53B~z-VS$kL?$Y;|8BPj!yFwA4MS-6GKH$gY_B7!nKY+0PTD<22^x1Fm zxJbAd>3*BP?*qO{y_M$abc;UTumX?z(!^A|M0{p;-`Wb?mzlMa>bbZR0{he32)o+Cmc_d3v-TTJgP+I_d^@JcDxr=>?6pl z#96`%9-8?<{7=hl$CwD?3Q_N!_tQSND{*q$<0U7qRZzZ=2F=4O(S$6u8@^Tu`MC|N zM8S>)lQ$iz6ihXp`Kba=d9LuOQgk0c%3MkXgv?5O6f4fuzBYw?mRx7xQ2UeU#Cft7 zmDp7IpGMI3PhucnjE+{qfb5qh>Q7=DvzT62B1wYMKZ_JzD9G3P*~oLu7lSX=Is?6a zy&jF$9&-62EXc-sK8Ad+^1$>}@RcK z>dpV#$EjrZ;PAx8E#y_kKc67{Cf3QrpHQ>;tiU>yB^RqO!yC6@Bs z{U^PjY@;E)RBrooE2{mt6Ml=C4>?Ho0i&iEFDkA7M+|2lmZ-oG)~nQiM4un^D1rC^ z9CI9H{a4s78s5cwOF7sr*2Vj;sL$uViSOb1^3unv6EzeaJ$R%00FN~vBM08lP#pct z!Z~pv&0PBCG^?R_Mz+z!qz@=`=#+QtehtM8UWjMBN5&s(>(1Jm3Q+toqSgn5t!p)B z$|X&OEemCDE3lNNQvSCz6@?s2%BsK@lJl<)sHt#fKfrBydC=& zwrMGDaO|>W1)lKyyR){U55ECdDi9Oi_4=YTZABRSMJR_=&OVgAwG^iU$f)wU9RJhi zbB?jmwHybhB-xx=S4%O4oNe7!l3#_F<)T`O|Et%Ngg=xJt5Wy(v-jL<=l=@?RLF@P z_a0_xTjv{9C5S4UlfUoJug-P<(v-5l#_s1gnumWTCbDE5?t&M^L)VTWox;u_cu<%n zHd1TG<(tvMf+OPZuOP5#eq5DyjA&bL%;zS3??Lil4#bFazC(&r?mmXj*zrVuj1g0N zd(2p6{{qjBFPq~{|Gv(|*f4M9EBMu_-(~$F`tuXZ$2zxigrW51$B2}!rGxb1o$z!b z-LPB6igte^+KzbRg8gx4XaX54Y8rQ3wqvgg>Z{1PyeCeak9s;KakwjXN(N4xXvn!( zH+T4X1>fqtDPF8Svg=gSn(uIc?LN2Tg@3-Fuf7+05GH5}%}-7ZZ8mP-43Mv$&ov1m z&HelxqjgGjT|IMj3(EvClPH^Gl6}$6zE9Edn0WDtG34gHa2{w??$S11Vo33wVGYEAkYLT{bsYV!Gsi#)=jdryl%m?E@Ot(oHy4{Bdi$ z==y6<4Bk%z`Op|qwP@%JU0n}IJMV@Vv5~Pea7D^L9q4mV6OESpol(1M!vXasM2kb8 z)B*E`zQs5H>vXb-6sB|w(|_s+hXdD%L>(pKCvZT?2|_*ZRZjb8ar>J^f77|Hm`UXRm!|!9NXVX)Gtwzaq`msFcypi^2GrT^dE;laND4}+ z+soG*oBrHPdhQr4rpvjD5*-&UnOgUBlt`5aB_-UY4?Rkh(v)<|SU)Un+&M6~OQiT5 zv)k?M7=O%dc4+H}oe?7Qn#cEXI|5L#cZ_5DdlEU<0lVMt6U zIZuiZPh>td1QC||`fc!u5PH(b7>cCs_xC;ZjSw46Eoh<}4$FXD-!w8KM6m9?)dr6u zFoLE?&&EcH2`PvA2Azt4=ETufqpn1YlIp(~T8zL|M+?S>7+Z(ojh@-jSr1~!2h1R{ zFi_!2>C+Lhf=}zlufs8g>{1_TK0_i?CzCG`s5`8;WktIlkdheOZCWrqsz1&1=v^p$91EM4 zR|9{W#tXeDbK|SaW`Il1*)i&@&6y;8 z=VQ*%IMHieyBbd>CBvAp6Kmpx(~Zq9>~&ILcqXRKz@G79jfRan$e8$6Ivi|@7e`(f zC2ic7f)_SFtfNjP5C@XJ=f6_^|5{V@xLxiHyHt>gW8_kr9}XKf>5|U$G z3?5bMW7)?v7<#;TIrMtGkQ>s_G%#J&v`f5@{QnE7@VQAun&olg$+&i*ugp@h_iFaz ztIueCXu4=*ylyIcE~{3S-nxrkE=j;-jvZsCd;vA$XSH$5F|*PEYsV-zL$)MQ1b?M&IV5-l#>ot)9_bUMBbFL^z0evH`g?tam`iZtvg z@?5l_eXQ78mANClQ925_YT|l~Xv){zr|EF)^SVQ$1u=ptrvd9T5Ys?ib11EkbGZg) zU?v{Ysda+E^)#>KWcLCyIC;o=a+u039X0Hx0sC|%|@>592+lm#yPz`7?O?o zMA6vn7cZ*MT`HWc*QnR&5L)l5c>`E2eKGXs<+(WrF6(Nf&2@?w8D?w<%RyxuGCZD& z7gK4>{JNmSEgcbG%~Lt1#yMYq=PmE!guJ#^sxc_d!v8~OdR{q(+B+9hnn#9x zeH1IQ*Df|`LOyYGzK4>o-J0S9+0n$ynnVI5b=-?P}CXlsD1};o6JygquH`H_l6Xpm~pN1!o?ou z5LgvduODq9MK33Ly) zbwVLLNZlAm&vUg8stAj(3`_o&6C)P$AS=BHZYiDhtp>!1V?>n<9#@PnaZPJX{SYOJ z=~S)zigL_d-|xAa9Vs4A+&rvv2^KWfy13;7@iy28zk#^nbe&Mu4;Q2VYBlR#QjCL) zfN2^k2FW}|G1&cLFg{o?9mjcfF&ubIc`{It*=PF@TrB zBfjD{HAPu`r2CWX;qVn+?W?jU?kYyM%xU|IWqi;6Q4DU{SZ+YRm{a|5`&bEVj=S{P z<>yzOTgWQGdZMd(mQcQxkCVhl-W3t`cz3rzVf2HD0y>|WB81#O3{{O z*Y%;|vdmSNVPJg*Z-xt*zq(b15au$x4i`1#y`&trQpb!GZ5jXJP!1DnX(&Ua_}^K{ z(|>=X9`xtsxO2e!>ibQtP+`D01UjF!-!s}s-y~RQ%MZ95Zzs0@*)}{-RP!#%Q9W-8 z5LLnZ6T^-ZcY=>k5B-JAu{9{e^q=t3fttXW622u>aiEPaC5mn4{8B+{BUtXxZ#DS?i z$SXqLA`Zm96BKGEzgH2qRP(rmwAU6CqD9do*M?7Ag$wJMXA7}})Lu(xr)qx~`Gp@a z_N?C^533{iB2druonyD2T};@-^1oT8(|~l*g3oa z`{Frx^-K`m!Lsh+0vv2MscgVa%E!t=sR9f{sO9f=Hv|C=>jUW+tjRi}@IAd{6`TLnf$0dt(NZT z0s01g^ zpA0R+f%d0AZPoO}FuU7rOMHrOWa`L61Go8t3Z{tAQHpy+N>PMFGv&5~Oer>r@D62t$u>DPp$V zTzBbcC@Ljyr38MLxx6kE($@E@7`IlnB8yNczQ<^PG2UK`ttB->hn@;Ub6U>c(f4;2 z*KfQHgS6dmDuNfO@QD%6u5M+d1@SSQyq&b4BFNd}MvwVs^JG6z-B&laF-T z4xTRxgCEDLuH?b4<6egk{>0&x_rg5<`Mf6KfORBPazDyLN7;iC1sQ*uL+iinmx%&5 zkt|c?;u)@0$~VOBb`Dms zPWU7NJYCB)%R%tO2&XG^6VW*7@8F1LIe21t?zz*|L^yMF$&UQI9Fq=B#A5cz4pUWg zXl#-&jB6r#st_hO++>98{oAHXL^fg-J?;#Cnv5UjTtt?IZkY`3OvXmN_tty6XCdhH zm`AHKlQD$z5z{g;Gk^cM_Rh(eza^nw{O$}i?r_|+g=q>(7SjYVDIMY6Hq~-(n~LaX z?sjFK=|J4~iXJ0V;X5nDdf)X7Dn{@zXKN<5+PUsbx{BbJChG&g6KF1!ULL%9hsIWieTH}$4Dd^68%^NEG zobS73bOSn9%IBN%p1eKvN=C5c7pNeQ!xNLRGIqzVzV0eWl! z2dYqwXK618ahih3G^RMGm0O*6FG#@#o(uKM1=megM5cmehY!nh@wH&k0!pGh<$Uz@|us`R3 zr%;D%RdC_^9X;2zq_0K=9mzY&!ND|}skPr`!?+!TwTS1=b~BGG`1a!9c@7v~buT#+ zA#y&Q1AlToyv@&m#KPCj!4$?id`hp*BZsT7z_Dt|@=ocH+sYmlrnNaDejP}|n#O$5 zS0U-%`(c%LQ(?j5K7|SgHD(fNB?S>O*P?=Q$3nM}<;j@4bcw!AU^X;gmgMg?O2$}? zt!1g(vmvW5CM02uyg+5cY)_}q?FSRF?Y?8L)T(UE8MS-DUF$@sepi~g(C_85&3p56 zzXTY#uIb)uZ#Hb$&s`pmkFWEsOUtuCVpVHB+XSpw@yTLyR~0^QBf!YkM0}>y_wUsz zyd-tV+&Brblv*ghtwI&|UVI~d>XAjeKPAxjTISA^v6z~(=C{m2Hrrr!rJy>Gn**Yx zlL0Ul<9Hojmjiix8JGsf&c1k;gBVjoiyB2~pl}^cz|~l4go=q*(xKk@`)FowHTcx} z`e+8y88txb*z4CjZFTmMo*?I4YWhBsp>|vrh;rC8++Gbc(qJ-gl ze?K38rbY@agEwx_>*Kg`sTw;dRrHgdPbpK5&4m@k6I(jwK#Tj44${xk@1BEeLG;2S zbKyZ_;Dx~|uv&7XejZZm5V-z?3Y40o0c#%gY5chRlM3;?z>vOhNp=k7ROh%`=HrM~ z;Ke@Qa;kMbQ|j^hDFaTZ?!(N7h0J&4!j@Fiv7|3FW1KPFOGtF4Pb?p%G)`Je++Te@ z&(iwYEsN+J!op5II7X-cUoo2SIv$*lC!9Y> zC`KI)R7VgGn69-eh{q;n*iZRb@Pjx0Qg90Mbn^oEOTDlZYnhX^hR#3o`76Wui6vV7 zbts42eF`Ce%D`%h*5N|%Z{cA`ydsWMbuL0{K3Ix!aGJlOs0fKN2UU&^eEcWQQZJh= zwDP@2AFeUdB97BD{sw=*rl0fX-o93XKfG=~`hbyAA1cKnnFp@Gu&pF~n3duc#V;Q_ zRA4jLL5?bg5$9obDlxilKT6)0K<031j?S^Q)0QPzDf6Y2dE)%w!4mA1x!R97t_*uL zc3?^M_owJ163rtLzd@VJrI;;?4!?k>Tg@jD7fKq)zaoNz zhW2GxM6Iip4Zfir^LFo-A%YhW`)|>tA2r1)yGfr@xR&f66UNg&V!5>?EHm>GIx=Scs@|1 z|HKzLdsY*lk_5_uKT&;u%Yhd!1VKMwPBG*DaVq@gfWoLBV1%P{BNZlvP%~5Fcf8@; z&*UuJAbH(=|96}roW`wL8K};gRUv@w+zr!lmATLdt8m;SGyLGkRQ%meygrR898`6i z*0wAKHCgXI^bHHC{qXRR6!bEGQo~F88x(vnou3RD%YOA0yNKG-IVbVo`}Y-EeTSaQ zqV-&k`%Au{+uVkJ8RqdwB|2v7>t7J2)A!QaW^w4nMIjx(pwgKSKCzJFBqw@4$d3GC zRV;kCN_OmLxJm9{EN)Bg*C%XU7dBPc#Nq_!zIJ|ssaCEku1_4EP57_qtL`UgUSAxj zc8dd3nIAcQ#Gvt~KWZDtRp(t+e}v3645i*4wzp~#zh!{<6W2Z#%SpaJ;8BU)tD7}& z7#ss>$F``%-KXDIF1!qv+%}J+u5ZpD|t7RZuN$t~3`+(Kg)|xf)3;?U5Z%jTw z#*@g8A|-hWKbFGJ|ZorZtnpFyLXb3kS2OjM>vdvDJC-Arb*5*|qAAO2{mWtH$b)7Nu_uWaBKrNQrEyxGw`fV<>N&m3A+dV{ zeFe+=&tuAPX?FJav2|%KA`i-C*niV|)6mTU|8oANn6P4k$Iq5Qq9mFULg!1N_m68< z{ktv`P$|QYY`H*&T;yfN@yu>`(m#L*Gk>)hAZzV4eeOJ~?GPy%} z5ME!;=ONIj=eB9{Lc}wQGY;>U2dR&C3>Lh6|5ucYa8jFc5h!ZOxenz}>LPoa~3`idtqPo+lcXui29g}LOF(0PgXsmtC%%H*5XFwXGk z9e3`Xpr?`er)ub)8AT3M_kTRn97s;$0(TKk>7{W~b8u|Ql>2XQyNWpp%ZyGo&cOn? z4>^g#%N($#ekBJdlE~hOU42a^n%z@Do2so2?z9&(m?yYgh1w5S6in`DBdYysDm;?+ zvB!dIrF!0_^HNQR?rVovZ;|IWc zoN-Y58h9Q|BAtg?{7!WwJxg{`rlUjhl>c5^$KslT2ApQ;@UD71+$bg$YRiS_)p$EG9Ql#BHupmo4)Z*IqQ z3~W^+qlZ~I*4<2-S$l69hMA zJo`n>ROsC9c>S7ZB(@N>e7RW)+Wk8Bu+P;f=v$5}jzcnZ8+ix%?~g`P*?*V}0@43I zx+DgFcABT1pOuWE+{H)xaarZ4A;Ywiu`qQ+o^f(4EKa)*wmzSPhRl=P5r^_)J?eM< zlL*;g`a2FSSN~qKVhH)nX0>&4B2GWoUHWJyqDlCgMiXfNH>c^$B7;PnTzP$X-A%FR zRWoYl+m?xVbJx3Hnq3TbEREVUuP6bUCf?`1wU%vxJG!>LCs^}NK+f3iih6E zk%X|1g63S)3p-9GfLssc0Ej~CiOzYI{plXuNt4PEB0*5Q7xP^baf;XJCsA-5)Td(? zy(Da4EbHt@Sdk0g@;voh%K2CXiYBMc&U}&xZs%$m5r$;Ouyjs>;n|^Cd2XS&Pg9`a z>?GtHF^D(}&WplVncPpt2KSKROEN>To-2=+rl4!w@4H8Tg<|3A=gxl&QW146DO=Yv z3`d&|JF(a{6{Zpl6^?Cm+*6K7Ly_y4o+ih_{`pzbaIX2T=2!oOfIp7-Ch6$O$BY5= zGaK%#==mGfA;tYq5W;1CECXBCe14L1jr{a%BT3Fcwf`UliZ-hQv>#?dpLt(3Lh(fA zWHVt$ZA}_mg3)MPnV!Ni14oa}bl%o07&)U_`&wPgK+1UCitkOfgQU{#=>E^8@1~%WcFWPnbZyKsB6>H<3zeDk-8RHgPrejUFfDV-(y>OjnU*B}%9Pyw8*fQJv@%4gq z%%}9WhJ`OS_qA(~Yo3ABh<3IouJ}RP`_5$GRF*ntcx_+A9G#vOJcH&6geml|P*$&> z!5K*2^Cr;z0kDl1o@E)}7QA8Ky#D1!G9c@Ezj`1qHJ=oVbo|ucNKb?twon@K_lk5} znDVj1p+IM>BD(WRq`{h6F_(65!ttAv>|F1pfha@whxK@a;iQ(eZI_Oxslyr_`~DIq zMsKOCxSx(Mn;KMHUGM_)46P2BS)`+UseYd?XwB-7NOp{e3_TF8^s$9ktB_nIik1&3`M_Z=+yiC4GMTPHv5O z>nJLZy&M0iM+WAv(V>PBM=^LL1q9cp!;*Q&K8}K?Ww(sz^BC?p)%t_IFm#-z5jZ{_ z#@g@Cc57%O7D+!xI-Y&(@y)%{V{x00rS$z8*m&9HrQ)VgP+QRLqv<$&&URX);*5btmGUC<(wj#9w!-DqB`e0{lNom1jm zPKm~8&tzm+ZEZW&aJkUn9REyOpQlegW~w|0SB?YFe0@y&i$gVrpT=gVJ&wJ{Bw?mJ zs9B?%EZd^>?S@Ca?Gbmdp2j(udUQ{uRAcu(k6;t*+Q1|{0c$37$|+p+94Gk=^hrQc z&&?)18ooj;HwqGWN`P67>(6KYwxPKt)pVYShgnV6tgy|`k;I*U@o{+cXWt3!ZyYeCcYv4^_osIjRd@tzRI^)@>6qP&eR0-=K2*K7dUXub!a@~*!!m|(u{P@ADa}9 z9;c5u9IWAizVbdBkG~Qe;0cL~-x~+VTAV2J!kFJZH$@Vs!`N=@m_BdW+k?cwM;B z6_e@Z89fMxq0E`MVPgE++KcCf!9AG`RBkZRtDtu}6lszV;fCyKvj@Ji4aLWp5pyfj zonh8Cx@Y!}FtqU6wcdZZ3j#RLG%Nx&mhH?Nzt$NO)^oWk?USaq$Gd!fi`hJfco&6W znHP3~(c&5>_Wg*4u_HU5$=}7cq2)1X!u(TfS1gFxzp{=-EI!IPvI{z4EKRIq@tN`V ziO!(ggmu~cI2@4qEEjli-k(k={m=7;;e^BNhguSkNVz{a;%|TQ=e~}^K=0y9PE51=<{*it)KTXiRk@DTg1Xp@=I@HLW+Xgp5<^-e9onkQUb~B4fUmha1Tfb*t2F1X<>65Bax;G$WPluyXb7*6e z#b54YqwLCxf>r+H(5@#Qp@a0*MWW4Rjl-i(J;mSP8`FJr!m(&-ch&cY&oMS>*2~up zVX!zcXL(93Tg=P5Pf4s$9DE*jH$K(|L&JLznL7l>w0o~sS=!;8_GDF3W)MsW|1|T? zE6^pq&G@W96b7_T3$FKuxT<5mj=Ua#)Lsu>HFtOeMR22T7Bsh5Otn5%9i5O+XQ0F3 zWIuelmbhbhtP|8fb{{NwW>p;d{9(k|Ndm!`TfDI;QJ@~BGCJLL%oKS3#8MISv_CGK2iU@U=xTZ z=eb}U=p*}ig0Q(>cm3J3y`ka0z-xZDLcGl?S3$IgLvKNSm+}<}a$`4bvcR!vuKLE989Q?eT@@VNN zXEYh+kKo+(-;ccX!?)(IOcECoXF6!Zr@wZ-h}P~&E4(jyTX0~^7pUMUX5~@=Z2^r!CFVo*&>ectbLdokPUNzMTb~x;X zUySoH4!}&Vn@OQu=cJ%MSB3{5Q0mBT*vd6GWdT@z$ebo`E?7yiLSy>(xj~*CR=PMr zlgG&u1F@hmyXViHj%Y1=R0DCM0iQ@5v4$w|FJ1;iXvuvF|96yhnxx3zDB+K(|0Jfg27w~BF1AJM%o~o1mVRbIXa7C?~ zP)b2Q06Ve+=XN*0iql;@Z^d4u+wR%)hNi!+!d>nM0T`e7&OhVnO{~lM8(deAWlKdGZE_7Y*T;FT2_$zar zfU(T;UJ^ItfXoM*y9_o6+jCuL-{AsaFVt>1CC}vfJ(2W?fnx7q@+j7?-r^_XFxLY2 z^TNcycDrhvwGsE(&bq=2SEY~2PUzLXaba+@52Dpo;|vXK#2lh&)x1Ky_2S=`LQcOB zdsYqEu`kOPJ>VNRN}GZJGzKpgVwzsLH^-q@eCAk7u#j<%?A!;6P9gn^iCY-d`w3r zUS^Cjjo0j00|HLMgwE6-v>>-QlKAmNO2r4`;Za zU|eL`5odTF=7qywEJepz$(*<*$bjs%JNL;kjCJKP2x_%S&dBw8EH8k5uL36_}LkG9+wEXKUr@ZKuu zKAuWlm3RsZzFj;~bQjNfEN2xY7I(Q8JuK`B4$-?lL?c+#(eBy*z?&o3d~w#z-X6iC z+Wk)bI-RSY7dwTB11U{*PtQ9p1`>7V;n`5(WT}6>r_&YjbN!PctwO@Yf(Ea0s~%k! z5leK}j!mFGtS5d4&$PNPIvu%|;^G}9=5C$dAnnUzksH%K{rjPC@l(%jbql?F;?B}* z4-pd~LKr)*|3D0(DU{>U2+_vvv;WH5*TitDQFXQr7w7Cwj30LKw76yPGv*X=KO$o< z+YPF{RW zCWC$bfXkeZhIctiJqCU}*Rz?{rnzLT%s;gMV^#<&3U2C`idKj zZ++ENVaGm-?Rts<9Eg3cuNX-*$esIj6i>N6W1WGb1@q<`H&$$udYpkmvj+PK8z^p? z>mE~vHB&QC){^0ct0CkOh7o%HrarW-20w$;j;M0(Ec#~nO=y)aTN zmPIv%$Yap@3?qe#p9`Ol_|S_6j+2ZO1uP@%{th>;S$kxtxWc^8qQA(OIiKc=wXBaO ze8&v-tM+NGXv4bR>2EOM;=wQjg*VquUa3L_+ZcE0E8a3Nu*Xj%?xDotx~7WV`A)+c zzyAx#=kD87q2_+-8XC%#Y{RIdr;z;s1O8%OZvsGT=qe0&%zEQ*wXRfOM+;;nTA$#%SD$VSs;8VGeIo{n2(n9gIn-Cm`*VDAh2*7W zG*s3+P2(J2BgIk9ajwx(9^^d8kQRzzl1tiD8Amb3R<nd9o-RQX)? zVb@f867JG0QBS##eVl`}m6o*{{?>b?t31Pop+U8j>$$G$g^tpY^L0CGEB_MJI4w{| zxsC2MkKdxRQ_1G90;tT%$53vuJRjU zA@?-ZP|A3?i>`7X3m=}p(41NizO2ww%DBMoA82As3+URW$_<4-Z7bux;(u#>rHVk< zUn}uSc8lmMFH5d7`I|YnII5}Ag1>i$1-Nu3GT6mYPdSh4hT?NEifh9ZP**S&OG~2{W>MIBFecrjAqP;hf`uEmT4x&n~VOts~+VSAi zi=Izf7OWaIQrta2Bv{)z%GVP3+E`JWan?g>D|<3Hb8HiZ=lK;L&#bkTVf^1u)lqa} zVD}AeB~v8M30=j=#%!e4QjTPBZV|UPqC7XeUoY^FK?m;#<8j5 zFj40YHTa8Q>HE=FXuT&t^vYlGm;LJoit4%7Hw+{^zwOm#idk%U&`5ZB&g(8PB*SSA%=WPRvSK%!cp)AOq)KDf<4cvl))(XjC z*RQEu#YT_xwu%oN+gPooyeP3^?G%rs@8mbe(Q={DsJ-GE=Me)c@V~Rniih9tdS)iz I|DIp|AB{f-U;qFB literal 0 HcmV?d00001 diff --git a/tests/global/global_stereo.shx b/tests/global/global_stereo.shx new file mode 100644 index 0000000000000000000000000000000000000000..b2f88bf7a5ca0dcba4ed275d458ea5977a247aba GIT binary patch literal 580 zcmZvYJxBvV6ohAUN$zsSM2{qf7=wQ(SOm1Olhek(#!AHMPzYir77_x&l~$5UE$sva z5epFuNfoRtEG$xpT8Nd9&PtufZX@J4?z?&K?dif}MJaQ8hIfWPnre<*0MwyjOV=YMas`O&>U!i} zYs)jN7>eH-5Y1^F-8Q;O(l=)2V=^!Lnhb!r%u~lV$Sv}bq>nwJwSA=f_PeH2Hms35 zfb&I$Ktl2b3@Zb_ZCQvkv?hirzC#}q?bs} zr^W99=#l(f*IOZ90Z;m;B=acxXG90&CHW0x&sleX9QRXspzn$G5Xi6cZh5~2_GP|5 i#Xf(F9zWc#$n}aW(a)*_LodugzsVLo$VHwv8vX%0NJ*jq literal 0 HcmV?d00001 diff --git a/tests/test_README.py b/tests/test_README.py deleted file mode 100644 index c1c7811..0000000 --- a/tests/test_README.py +++ /dev/null @@ -1,8 +0,0 @@ -import os - -current_dir = os.path.abspath(os.path.dirname(__file__)) -parent_dir = os.path.abspath(current_dir + "/../") - - -def test_readme(): - os.system(f"pytest --codeblocks {parent_dir}/README.md") diff --git a/tests/test_bathymetric_gradient_function.py b/tests/test_bathymetric_gradient_function.py index b4c3865..0bb0d72 100644 --- a/tests/test_bathymetric_gradient_function.py +++ b/tests/test_bathymetric_gradient_function.py @@ -1,5 +1,6 @@ import os +import pytest import oceanmesh as om fname = os.path.join(os.path.dirname(__file__), "GSHHS_i_L1.shp") @@ -48,6 +49,7 @@ def mesh_plot(points, cells, plot_title=""): pt.show() +@pytest.mark.skip(reason="not implemented yet") def test_bathymetric_gradient_function(): EPSG = 4326 # EPSG:4326 or WGS84 bbox = (-74.4, -73.4, 40.2, 41.2) @@ -82,7 +84,7 @@ def test_bathymetric_gradient_function(): "Bathymetric Gradient", "Feature Sizing & Bathymetric Gradient", ], - [edge_length1, edge_length2, edge_length3], + [edge_length1, edge_length3], ): print(f"Generating mesh associated with {name_}") edge_length_ = om.enforce_mesh_gradation(edge_length, gradation=0.15) diff --git a/tests/test_edgefx.py b/tests/test_edgefx.py index 3afbb35..7dfb379 100644 --- a/tests/test_edgefx.py +++ b/tests/test_edgefx.py @@ -22,17 +22,17 @@ def test_edgefx(): dis3 = dis2.interpolate_to(dis1) - ax = dis3.plot(hold=True) + fig, ax, pc = dis3.plot(holding=True) shore1.plot(ax) dis4 = dis1.interpolate_to(dis2) - ax = dis4.plot(hold=False) + _, ax, _ = dis4.plot(holding=False) def test_edgefx_elevation_bounds(): region = om.Region(extent=(-95.24, -95.21, 28.95, 29.00), crs=4326) - dem = om.DEM(dfname, bbox=region.bbox, crs=4326) + dem = om.DEM(dfname, bbox=region, crs=4326) sho = om.Shoreline(fname, region.bbox, 0.005) sho.plot() @@ -57,12 +57,12 @@ def test_edgefx_medial_axis(): edge_length = om.feature_sizing_function( shoreline, sdf, max_edge_length=5e3, plot=True ) - ax = edge_length.plot( + fig, ax, pc = edge_length.plot( xlabel="longitude (WGS84 degrees)", ylabel="latitude (WGS84 degrees)", title="Feature sizing function", cbarlabel="mesh size (degrees)", - hold=True, + holding=True, xlim=[-74.3, -73.8], ylim=[40.3, 40.8], ) diff --git a/tests/test_geodata.py b/tests/test_geodata.py index f8f2113..059fa57 100644 --- a/tests/test_geodata.py +++ b/tests/test_geodata.py @@ -43,5 +43,5 @@ def test_geodata(files_bboxes): f, bbox = files_bboxes region = Region(bbox, 4326) - dem = DEM(f, bbox=region.bbox, crs=region.crs) + dem = DEM(f, bbox=region, crs=region.crs) assert isinstance(dem, DEM), "DEM class did not form" diff --git a/tests/test_global_stereo.py b/tests/test_global_stereo.py new file mode 100644 index 0000000..cfd7c8e --- /dev/null +++ b/tests/test_global_stereo.py @@ -0,0 +1,60 @@ +import os + +import oceanmesh as om + +# Note: global_stereo.shp has been generated using global_tag in pyposeidon +# https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452 +fname = os.path.join(os.path.dirname(__file__), "global", "global_latlon.shp") +fname2 = os.path.join(os.path.dirname(__file__), "global", "global_stereo.shp") + + +def test_global_stereo(): + # it is necessary to define all the coastlines at once: + # the Shoreline class will the detect the biggest coastline (antartica and define it + # as the outside boundary) + + EPSG = 4326 # EPSG:4326 or WGS84 + bbox = (-180.00, 180.00, -89.00, 90.00) + extent = om.Region(extent=bbox, crs=4326) + + min_edge_length = 0.5 # minimum mesh size in domain in meters + max_edge_length = 2 # maximum mesh size in domain in meters + shoreline = om.Shoreline(fname, extent.bbox, min_edge_length) + sdf = om.signed_distance_function(shoreline) + edge_length0 = om.distance_sizing_function(shoreline, rate=0.11) + edge_length1 = om.feature_sizing_function( + shoreline, + sdf, + min_edge_length=min_edge_length, + max_edge_length=max_edge_length, + crs=EPSG, + ) + + edge_length = om.compute_minimum([edge_length0, edge_length1]) + edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.09, stereo=True) + + # once the size functions have been defined, wed need to mesh inside domain in + # stereographic projections. This is way we use another coastline which is + # already translated in a sterographic projection + shoreline_stereo = om.Shoreline(fname2, extent.bbox, min_edge_length, stereo=True) + domain = om.signed_distance_function(shoreline_stereo) + + points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=100) + + # remove degenerate mesh faces and other common problems in the mesh + points, cells = om.make_mesh_boundaries_traversable(points, cells) + points, cells = om.delete_faces_connected_to_one_face(points, cells) + + # apply a Laplacian smoother + points, cells = om.laplacian2(points, cells, max_iter=100) + + # plot + fig, ax, _ = edge_length.plot( + holding=True, + plot_colorbar=True, + stereo=True, + vmax=max_edge_length, + ) + + ax.triplot(points[:, 0], points[:, 1], cells, color="gray", linewidth=0.5) + shoreline_stereo.plot(ax=ax) diff --git a/tests/test_grade.py b/tests/test_grade.py index bc9756c..eb50bef 100644 --- a/tests/test_grade.py +++ b/tests/test_grade.py @@ -17,7 +17,7 @@ def test_grade(): edge_length = om.distance_sizing_function(shore, rate=0.35) test_edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.20) - test_edge_length.plot(show=False, filename="test_grade_edge_length.png") + test_edge_length.plot(filename="test_grade_edge_length.png") domain = om.signed_distance_function(shore) diff --git a/tox.ini b/tox.ini index afa2657..f6aac90 100644 --- a/tox.ini +++ b/tox.ini @@ -15,4 +15,4 @@ extras = all setenv = MPLBACKEND = agg commands = - pytest {posargs} --codeblocks + pytest {posargs} -v --codeblocks