From 9e8f1d63f5bff242d00a56d4cfa9b819ec308815 Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 14:05:00 -0400 Subject: [PATCH 1/7] added band1 and band2 to virtual catalog class --- virtual/catalogs.py | 4 +++- virtual/marblecutter_overload_nd.py | 0 2 files changed, 3 insertions(+), 1 deletion(-) create mode 100644 virtual/marblecutter_overload_nd.py diff --git a/virtual/catalogs.py b/virtual/catalogs.py index f291d86..aa93028 100644 --- a/virtual/catalogs.py +++ b/virtual/catalogs.py @@ -13,12 +13,14 @@ class VirtualCatalog(Catalog): - def __init__(self, uri, rgb=None, nodata=None, linear_stretch=None): + def __init__(self, uri, rgb=None, nodata=None, linear_stretch=None, band1=None, band2=None): self._uri = uri self._rgb = rgb self._nodata = nodata self._linear_stretch = linear_stretch self._meta = {} + self.band1 = band1 + self.band2 = band2 with get_source(self._uri) as src: self._bounds = warp.transform_bounds(src.crs, WGS84_CRS, *src.bounds) diff --git a/virtual/marblecutter_overload_nd.py b/virtual/marblecutter_overload_nd.py new file mode 100644 index 0000000..e69de29 From e049f83a9c9fd3a8a526889419a464f5449405eb Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 14:08:31 -0400 Subject: [PATCH 2/7] added render_nd_png path /ndtiles/ and marblecutter_overload_nd.py for testing --- virtual/marblecutter_overload_nd.py | 34 +++++++++++++++++++++++++++++ virtual/web.py | 22 +++++++++++++++++++ 2 files changed, 56 insertions(+) diff --git a/virtual/marblecutter_overload_nd.py b/virtual/marblecutter_overload_nd.py index e69de29..3c93708 100644 --- a/virtual/marblecutter_overload_nd.py +++ b/virtual/marblecutter_overload_nd.py @@ -0,0 +1,34 @@ +# coding=utf-8 +from __future__ import absolute_import + +import logging + +import mercantile +from affine import Affine +from rasterio.crs import CRS + +from marblecutter import Bounds, render + +LOG = logging.getLogger(__name__) +TILE_SHAPE = (256, 256) +WEB_MERCATOR_CRS = CRS.from_epsg(3857) +WGS84_CRS = CRS.from_epsg(4326) + + +def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, data_band_count=3): + """Render a tile into Web Mercator.""" + bounds = Bounds(mercantile.xy_bounds(tile), WEB_MERCATOR_CRS) + shape = tuple(map(int, Affine.scale(scale) * TILE_SHAPE)) + + catalog.validate(tile) + + return render( + bounds, + shape, + WEB_MERCATOR_CRS, + catalog=catalog, + format=format, + data_band_count=data_band_count, + transformation=transformation, + ) + diff --git a/virtual/web.py b/virtual/web.py index eb7fe81..f35c184 100644 --- a/virtual/web.py +++ b/virtual/web.py @@ -9,6 +9,7 @@ from marblecutter.formats.optimal import Optimal from marblecutter.transformations import Image from marblecutter.web import app +from . import marblecutter_overload_nd as marblecutter_nd from mercantile import Tile import urllib @@ -119,3 +120,24 @@ def render_png(z, x, y, scale=1, prefix=None): headers.update(catalog.headers) return data, 200, headers + + +@app.route("/ndtiles///") +@app.route("/ndtiles///@x") +@app.route("//ndtiles////") +@app.route("//ndtiles///") +def render_nd_png(z, x, y, scale=1, prefix=None): + catalog = make_catalog(request.args) + tile = Tile(x, y, z) + + headers, data = marblecutter_nd.render_tile( + tile, + catalog, + format=IMAGE_FORMAT, + transformation=IMAGE_TRANSFORMATION, + scale=scale, + ) + + headers.update(catalog.headers) + + return data, 200, headers From afeee72630e39051ea600d1d550a3e7526f7ba79 Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 14:30:58 -0400 Subject: [PATCH 3/7] added render overload and bandmath function --- virtual/marblecutter_overload_nd.py | 125 +++++++++++++++++++++++++++- 1 file changed, 124 insertions(+), 1 deletion(-) diff --git a/virtual/marblecutter_overload_nd.py b/virtual/marblecutter_overload_nd.py index 3c93708..cbfc6a4 100644 --- a/virtual/marblecutter_overload_nd.py +++ b/virtual/marblecutter_overload_nd.py @@ -15,6 +15,31 @@ WGS84_CRS = CRS.from_epsg(4326) +# coding=utf-8 +from __future__ import absolute_import, division, print_function + +import logging +import math +import unicodedata + +from haversine import haversine +import numpy as np +import rasterio +from rasterio import transform, warp, windows +from rasterio._err import CPLE_OutOfMemoryError +from rasterio.crs import CRS +from rasterio.enums import MaskFlags +from rasterio.transform import Affine +from rasterio.vrt import WarpedVRT +from rasterio.warp import Resampling + +from marblecutter import mosaic +from marblecutter.stats import Timer +from marblecutter.utils import Bounds, PixelCollection +from marblecutter import get_resolution_in_meters, NoDataAvailable + + + def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, data_band_count=3): """Render a tile into Web Mercator.""" bounds = Bounds(mercantile.xy_bounds(tile), WEB_MERCATOR_CRS) @@ -22,7 +47,7 @@ def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, dat catalog.validate(tile) - return render( + return render_nd( bounds, shape, WEB_MERCATOR_CRS, @@ -32,3 +57,101 @@ def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, dat transformation=transformation, ) + +def performBandMath(data, catalog): + + if (catalog.band1 is not None) and (catalog.band2 is not None): + + band1 = int(catalog.band1) if isinstance(catalog.band1, str) else catalog.band1 + band2 = int(catalog.band2) if isinstance(catalog.band2, str) else catalog.band2 + + bandTile1 = data[band1] + bandTile2 = data[band2] + + np.seterr(divide='ignore', invalid='ignore') + data = (bandTile2.astype(float) - bandTile1.astype(float)) / (bandTile2.astype(float) + bandTile1.astype(float)) + + else: + ##TODO introduce new error + raise NoDataAvailable() + + return data + +def render_nd( + bounds, + shape, + target_crs, + format, + data_band_count, + catalog=None, + sources=None, + transformation=None, +): + """Render data intersecting bounds into shape using an optional + transformation. And perform normative diference using bands specified in catalog""" + resolution_m = get_resolution_in_meters(bounds, shape) + stats = [] + + if sources is None and catalog is None: + raise Exception("Either sources or a catalog must be provided.") + + if transformation: + bounds, shape, offsets = transformation.expand(bounds, shape) + + if sources is None and catalog is not None: + with Timer() as t: + sources = catalog.get_sources(bounds, resolution_m) + stats.append(("Get Sources", t.elapsed)) + + # TODO try to avoid materializing the iterator + sources = list(sources) + + if sources is None or len(sources) == 0: + raise NoDataAvailable() + + with Timer() as t: + sources_used, pixels = mosaic.composite( + sources, bounds, shape, target_crs, data_band_count + ) + stats.append(("Composite", t.elapsed)) + + if pixels.data is None: + raise NoDataAvailable() + + data_format = "raw" + + if transformation: + with Timer() as t: + pixels, data_format = transformation.transform(pixels) + stats.append(("Transform", t.elapsed)) + + with Timer() as t: + pixels = transformation.postprocess(pixels, data_format, offsets) + + stats.append(("Post-process", t.elapsed)) + + with Timer() as t: + (content_type, formatted) = format(pixels, data_format) + stats.append(("Format", t.elapsed)) + + headers = { + "Content-Type": content_type, + "Server-Timing": [ + 'op{};desc="{}";dur={:0.2f}'.format(i, name, time) + for (i, (name, time)) in enumerate(stats) + ] + + [ + 'src{};desc="{} - {}"'.format( + i, + unicodedata.normalize("NFKD", unicode(name)).encode( + "ascii", "ignore" + ).replace( + '"', '\\"' + ), + url, + ) + for (i, (name, url)) in enumerate(sources_used) + ], + } + + return (headers, formatted) From 1e6246096c5e806904dd0b21595fa8fecc7dd4c3 Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 14:50:21 -0400 Subject: [PATCH 4/7] render_nd overload --- virtual/marblecutter_overload_nd.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/virtual/marblecutter_overload_nd.py b/virtual/marblecutter_overload_nd.py index cbfc6a4..07d34ba 100644 --- a/virtual/marblecutter_overload_nd.py +++ b/virtual/marblecutter_overload_nd.py @@ -58,24 +58,26 @@ def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, dat ) -def performBandMath(data, catalog): +def performBandMath(pixels, catalog): if (catalog.band1 is not None) and (catalog.band2 is not None): band1 = int(catalog.band1) if isinstance(catalog.band1, str) else catalog.band1 band2 = int(catalog.band2) if isinstance(catalog.band2, str) else catalog.band2 - bandTile1 = data[band1] - bandTile2 = data[band2] + bandTile1 = pixels.data[band1] + bandTile2 = pixels.data[band2] np.seterr(divide='ignore', invalid='ignore') - data = (bandTile2.astype(float) - bandTile1.astype(float)) / (bandTile2.astype(float) + bandTile1.astype(float)) + new_pixel_collection = PixelCollection((bandTile2.astype(float) - bandTile1.astype(float)) / (bandTile2.astype(float) + bandTile1.astype(float)), + pixels.bounds, + pixels.band) else: ##TODO introduce new error raise NoDataAvailable() - return data + return new_pixel_collection def render_nd( bounds, @@ -86,6 +88,8 @@ def render_nd( catalog=None, sources=None, transformation=None, + nd_calc=True + ): """Render data intersecting bounds into shape using an optional transformation. And perform normative diference using bands specified in catalog""" @@ -130,6 +134,9 @@ def render_nd( stats.append(("Post-process", t.elapsed)) + if nd_calc: + pixels = performBandMath(pixels, catalog) + with Timer() as t: (content_type, formatted) = format(pixels, data_format) stats.append(("Format", t.elapsed)) From 58bb81669378882de6c84d8652633051cc8d0e70 Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 15:57:56 -0400 Subject: [PATCH 5/7] bandmath implementation - no errors but no NDVI image returned --- virtual/catalogs.py | 7 +++++-- virtual/nd_formats.py | 36 ++++++++++++++++++++++++++++++++++++ virtual/web.py | 17 +++++++++++++---- 3 files changed, 54 insertions(+), 6 deletions(-) create mode 100644 virtual/nd_formats.py diff --git a/virtual/catalogs.py b/virtual/catalogs.py index aa93028..6f5ed17 100644 --- a/virtual/catalogs.py +++ b/virtual/catalogs.py @@ -19,8 +19,11 @@ def __init__(self, uri, rgb=None, nodata=None, linear_stretch=None, band1=None, self._nodata = nodata self._linear_stretch = linear_stretch self._meta = {} - self.band1 = band1 - self.band2 = band2 + self.band1 = int(band1) + self.band2 = int(band2) + + LOG.info("band1:{}".format(band1)) + LOG.info("band2:{}".format(band2)) with get_source(self._uri) as src: self._bounds = warp.transform_bounds(src.crs, WGS84_CRS, *src.bounds) diff --git a/virtual/nd_formats.py b/virtual/nd_formats.py new file mode 100644 index 0000000..a57b3f0 --- /dev/null +++ b/virtual/nd_formats.py @@ -0,0 +1,36 @@ +# coding=utf-8 +from __future__ import absolute_import + +import logging + +import numpy as np + +from marblecutter import _isimage +from marblecutter.utils import PixelCollection +from marblecutter.formats.jpeg import JPEG +from marblecutter.formats.png import PNG + +LOG = logging.getLogger(__name__) +JPEG_FORMAT = JPEG() +PNG_FORMAT = PNG() + + +def Optimalnd(): + + def _format(pixels, data_format): + if not _isimage(data_format): + raise Exception("Must be an image format") + + alpha = pixels.data[:, :, 3] + + if np.equal(alpha, 255).all(): + # solid + rgb_pixels = PixelCollection( + pixels.data[:, :, 0:3], pixels.bounds, pixels.band + ) + return JPEG_FORMAT(rgb_pixels, "RGB") + else: + # partially transparent + return PNG_FORMAT(pixels, data_format) + + return _format diff --git a/virtual/web.py b/virtual/web.py index f35c184..374c1f7 100644 --- a/virtual/web.py +++ b/virtual/web.py @@ -12,13 +12,21 @@ from . import marblecutter_overload_nd as marblecutter_nd from mercantile import Tile import urllib +from .nd_formats import Optimalnd from .catalogs import VirtualCatalog + + + + + + LOG = logging.getLogger(__name__) IMAGE_TRANSFORMATION = Image() IMAGE_FORMAT = Optimal() +IMAGE_FORMAT_ND = Optimalnd() @lru_cache() @@ -27,10 +35,12 @@ def make_catalog(args): rgb = args.get("rgb") nodata = args.get("nodata") linear_stretch = args.get("linearStretch") + band1 = args.get('band1') + band2 = args.get('band2') try: return VirtualCatalog( - source, rgb=rgb, nodata=nodata, linear_stretch=linear_stretch + source, rgb=rgb, nodata=nodata, linear_stretch=linear_stretch, band1=band1, band2=band2 ) except Exception as e: LOG.warn(e.message) @@ -129,11 +139,10 @@ def render_png(z, x, y, scale=1, prefix=None): def render_nd_png(z, x, y, scale=1, prefix=None): catalog = make_catalog(request.args) tile = Tile(x, y, z) - - headers, data = marblecutter_nd.render_tile( + headers, data = marblecutter_nd.render_nd_tile( tile, catalog, - format=IMAGE_FORMAT, + format=IMAGE_FORMAT_ND, transformation=IMAGE_TRANSFORMATION, scale=scale, ) From e1e55b247d28936a8bb8aa6d5fd159d25875d53b Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 18:00:34 -0400 Subject: [PATCH 6/7] working BandMath --- virtual/catalogs.py | 8 +- virtual/marblecutter_overload_nd.py | 47 ++--- virtual/mosaic.py | 123 +++++++++++++ virtual/recipes.py | 274 ++++++++++++++++++++++++++++ virtual/web.py | 2 +- 5 files changed, 424 insertions(+), 30 deletions(-) create mode 100644 virtual/mosaic.py create mode 100644 virtual/recipes.py diff --git a/virtual/catalogs.py b/virtual/catalogs.py index 6f5ed17..2625c84 100644 --- a/virtual/catalogs.py +++ b/virtual/catalogs.py @@ -19,8 +19,8 @@ def __init__(self, uri, rgb=None, nodata=None, linear_stretch=None, band1=None, self._nodata = nodata self._linear_stretch = linear_stretch self._meta = {} - self.band1 = int(band1) - self.band2 = int(band2) + self._band1 = int(band1) + self._band2 = int(band2) LOG.info("band1:{}".format(band1)) LOG.info("band2:{}".format(band2)) @@ -79,6 +79,10 @@ def get_sources(self, bounds, resolution): if self._linear_stretch is not None: recipes["linear_stretch"] = "per_band" + if (self._band1 is not None) and (self._band2 is not None): + recipes['nd_bands'] = (self._band1, self._band2) + LOG.info("nd_bands: {}".format(recipes['nd_bands'])) + yield Source( url=self._uri, name=self._name, diff --git a/virtual/marblecutter_overload_nd.py b/virtual/marblecutter_overload_nd.py index 07d34ba..2e9006e 100644 --- a/virtual/marblecutter_overload_nd.py +++ b/virtual/marblecutter_overload_nd.py @@ -1,13 +1,11 @@ # coding=utf-8 -from __future__ import absolute_import +from __future__ import absolute_import, division, print_function import logging import mercantile -from affine import Affine from rasterio.crs import CRS -from marblecutter import Bounds, render LOG = logging.getLogger(__name__) TILE_SHAPE = (256, 256) @@ -16,24 +14,12 @@ # coding=utf-8 -from __future__ import absolute_import, division, print_function -import logging -import math import unicodedata -from haversine import haversine import numpy as np -import rasterio -from rasterio import transform, warp, windows -from rasterio._err import CPLE_OutOfMemoryError -from rasterio.crs import CRS -from rasterio.enums import MaskFlags from rasterio.transform import Affine -from rasterio.vrt import WarpedVRT -from rasterio.warp import Resampling - -from marblecutter import mosaic +from virtual import mosaic from marblecutter.stats import Timer from marblecutter.utils import Bounds, PixelCollection from marblecutter import get_resolution_in_meters, NoDataAvailable @@ -59,17 +45,28 @@ def render_nd_tile(tile, catalog, transformation=None, format=None, scale=1, dat def performBandMath(pixels, catalog): - + LOG.info("bandmath") if (catalog.band1 is not None) and (catalog.band2 is not None): - + LOG.info("Perform Bandmath") band1 = int(catalog.band1) if isinstance(catalog.band1, str) else catalog.band1 band2 = int(catalog.band2) if isinstance(catalog.band2, str) else catalog.band2 - bandTile1 = pixels.data[band1] - bandTile2 = pixels.data[band2] + LOG.info("pixels.data shape = {}".format(pixels.data.shape)) + data = np.asarray(pixels.data) + LOG.info("data shape = {}".format(pixels.data.shape)) + LOG.info("pixels = {}".format(pixels)) + + + bandTile1 = data[:, :, band1] + bandTile2 = data[:, :, band2] + + newBand = (bandTile2.astype(float) - bandTile1.astype(float)) / (bandTile2.astype(float) + bandTile1.astype(float)) + newData = data + for idx in range(np.shape(data)[2]): + newData[:,:,idx] = newBand np.seterr(divide='ignore', invalid='ignore') - new_pixel_collection = PixelCollection((bandTile2.astype(float) - bandTile1.astype(float)) / (bandTile2.astype(float) + bandTile1.astype(float)), + new_pixel_collection = PixelCollection(newData, pixels.bounds, pixels.band) @@ -77,6 +74,7 @@ def performBandMath(pixels, catalog): ##TODO introduce new error raise NoDataAvailable() + LOG.info("return Bandmath") return new_pixel_collection def render_nd( @@ -87,9 +85,7 @@ def render_nd( data_band_count, catalog=None, sources=None, - transformation=None, - nd_calc=True - + transformation=None ): """Render data intersecting bounds into shape using an optional transformation. And perform normative diference using bands specified in catalog""" @@ -134,9 +130,6 @@ def render_nd( stats.append(("Post-process", t.elapsed)) - if nd_calc: - pixels = performBandMath(pixels, catalog) - with Timer() as t: (content_type, formatted) = format(pixels, data_format) stats.append(("Format", t.elapsed)) diff --git a/virtual/mosaic.py b/virtual/mosaic.py new file mode 100644 index 0000000..1cc74c8 --- /dev/null +++ b/virtual/mosaic.py @@ -0,0 +1,123 @@ +# noqa +# coding=utf-8 +from __future__ import absolute_import, division, print_function + +import logging +import multiprocessing +from concurrent import futures + +import numpy as np + +from rasterio import warp + +from virtual import recipes +from marblecutter.utils import Bounds, PixelCollection + +LOG = logging.getLogger(__name__) + + +def composite(sources, bounds, shape, target_crs, band_count): + """Composite data from sources into a single raster covering bounds, but in + the target CRS.""" + # avoid circular dependencies + from marblecutter import _nodata, get_resolution_in_meters, get_source, read_window + + # TODO this belongs in render + if bounds.crs == target_crs: + canvas_bounds = bounds + else: + canvas_bounds = Bounds( + warp.transform_bounds(bounds.crs, target_crs, *bounds.bounds), target_crs + ) + + resolution = get_resolution_in_meters(bounds, shape) + sources = recipes.preprocess(sources, resolution=resolution) + + def _read_window(source): + with get_source(source.url) as src: + LOG.info( + "Compositing %s (%s) as band %s", source.url, source.name, source.band + ) + + # read a window from the source data + # TODO ask for a buffer here, get back an updated bounding box + # reflecting it + try: + window_data = read_window(src, canvas_bounds, shape, source.recipes) + except Exception as e: + LOG.exception("Error reading %s: %s", source.url, e) + return + + return source, PixelCollection( + window_data.data, window_data.bounds, source.band + ) + + # iterate over available sources, sorted by decreasing "quality" + with futures.ThreadPoolExecutor( + max_workers=multiprocessing.cpu_count() * 5 + ) as executor: + ws = executor.map(_read_window, sources) + + sources_used = [] + + ws = recipes.postprocess(ws) + + canvas_data = np.ma.zeros( + (band_count,) + shape, dtype=np.float32, fill_value=_nodata(np.float32) + ) + canvas_data.mask = True + + canvas = PixelCollection(canvas_data, canvas_bounds) + + for source, window_data in filter(None, ws): + window_data = recipes.apply(source.recipes, window_data, source=source) + + if window_data.data is None: + continue + + # paste the resulting data onto a canvas + canvas = paste( + PixelCollection( + window_data.data.astype(np.float32), + window_data.bounds, + window_data.band, + ), + canvas, + ) + sources_used.append(source) + + if not canvas.data.mask.any(): + # stop if all pixels are valid + break + + return map(lambda s: (s.name, s.url), sources_used), canvas + + +def paste(window_pixels, canvas_pixels): + """ "Reproject" src data into the correct position within a larger image""" + window_data, (window_bounds, window_crs), band = window_pixels + canvas, (canvas_bounds, canvas_crs), _ = canvas_pixels + if window_crs != canvas_crs: + raise Exception("CRSes must match: {} != {}".format(window_crs, canvas_crs)) + + if window_bounds != canvas_bounds: + raise Exception( + "Bounds must match: {} != {}".format(window_bounds, canvas_bounds) + ) + + if band is None and window_data.shape != canvas.shape: + raise Exception( + "Data shapes must match: {} != {}".format(window_data.shape, canvas.shape) + ) + + if band is None: + merged = np.ma.where(canvas.mask & ~window_data.mask, window_data, canvas) + merged.fill_value = canvas.fill_value + else: + merged_band = np.ma.where( + canvas.mask[band] & ~window_data.mask, window_data, canvas[band] + ) + canvas[band] = merged_band[0] + merged = canvas + + return PixelCollection(merged, canvas_pixels.bounds) diff --git a/virtual/recipes.py b/virtual/recipes.py new file mode 100644 index 0000000..92bee88 --- /dev/null +++ b/virtual/recipes.py @@ -0,0 +1,274 @@ +# coding=utf-8 +from __future__ import absolute_import, division, print_function + +import itertools +import logging +from functools import reduce + +import numpy as np + +from rio_pansharpen.methods import Brovey +from rio_tiler import utils +from rio_toa import reflectance + +from marblecutter.utils import PixelCollection, Source + +BAND_MAPPING = {"r": 0, "g": 1, "b": 2, "pan": 4} +LOG = logging.getLogger(__name__) + + +def apply(recipes, pixels, source=None): + data = pixels.data + dtype_min = np.iinfo(data.dtype).min + dtype_max = np.iinfo(data.dtype).max + + if "landsat8" in recipes: + LOG.info("Applying landsat 8 recipe") + + out = np.ma.empty(shape=(data.shape), dtype=np.float32) + + for bdx, source_band in enumerate((4, 3, 2)): + sun_elev = source.meta["L1_METADATA_FILE"]["IMAGE_ATTRIBUTES"][ + "SUN_ELEVATION" + ] + multi_reflect = source.meta["L1_METADATA_FILE"][ + "RADIOMETRIC_RESCALING" + ].get( + "REFLECTANCE_MULT_BAND_{}".format(source_band) + ) + add_reflect = source.meta["L1_METADATA_FILE"]["RADIOMETRIC_RESCALING"].get( + "REFLECTANCE_ADD_BAND_{}".format(source_band) + ) + + min_val = source.meta.get("values", {}).get(str(source_band), {}).get( + "min", dtype_min + ) + max_val = source.meta.get("values", {}).get(str(source_band), {}).get( + "max", dtype_max + ) + + band_data = 10000 * reflectance.reflectance( + data[bdx], multi_reflect, add_reflect, sun_elev, src_nodata=0 + ) + + # calculate local min/max as fallbacks + if ( + min_val == dtype_min + and max_val == dtype_max + and len(data.compressed()) > 0 + ): + local_min, local_max = np.percentile(band_data.compressed(), (2, 98)) + min_val = max(min_val, local_min) + max_val = min(max_val, local_max) + + out[bdx] = np.ma.where( + band_data > 0, + utils.linear_rescale( + band_data, in_range=[min_val, max_val], out_range=[0, 1] + ), + 0, + ) + + data = out + + if "imagery" in recipes: + LOG.info("Applying imagery recipe") + LOG.info("recipes={}".format(recipes)) + if "rgb_bands" in recipes: + data = np.ma.array([data[i - 1] for i in recipes["rgb_bands"]]) + + if "nd_bands" in recipes: + LOG.info("nd_bands") + band1_tile = np.ma.array(data[recipes['nd_bands'][0]-1]) + band2_tile = np.ma.array(data[recipes['nd_bands'][1]-1]) + + + data = np.ma.array([(band2_tile.astype(float) - band1_tile.astype(float)) / (band2_tile.astype(float) + band1_tile.astype(float))]) + + if data.shape[0] > 3: + # alpha band (and beyond) present; drop it (them) + # TODO use band 4 as a mask instead + data = data[0:3] + + if "linear_stretch" in recipes: + if recipes["linear_stretch"] == "global": + data = utils.linear_rescale( + data, + in_range=(np.min(data), np.max(data)), + out_range=(dtype_min, dtype_max), + ) + elif recipes["linear_stretch"] == "per_band": + for band in xrange(0, data.shape[0]): + min_val = source.meta.get("values", {}).get(band, {}).get( + "min", np.min(data[band]) + ) + max_val = source.meta.get("values", {}).get(band, {}).get( + "max", np.max(data[band]) + ) + data[band] = np.ma.where( + data[band] > 0, + utils.linear_rescale( + data[band], + in_range=(min_val, max_val), + out_range=(dtype_min, dtype_max), + ), + 0, + ) + else: + # rescale after reducing and before increasing dimensionality + if data.dtype != np.uint8 and not np.issubdtype(data.dtype, np.floating): + # rescale non-8-bit sources (assuming that they're raw sensor data) + + for band in xrange(0, data.shape[0]): + min_val = source.meta.get("values", {}).get(band, {}).get( + "min", dtype_min + ) + max_val = source.meta.get("values", {}).get(band, {}).get( + "max", dtype_max + ) + + if ( + min_val == dtype_min + and max_val == dtype_max + and len(data.compressed()) > 0 + ): + local_min, local_max = np.percentile( + data[band].compressed(), (2, 98) + ) + min_val = max(min_val, local_min) + max_val = min(max_val, local_max) + + data[band] = np.ma.where( + data[band] > 0, + utils.linear_rescale( + data[band], + in_range=(min_val, max_val), + out_range=(dtype_min, dtype_max), + ), + 0, + ) + + if data.shape[0] == 1: + # likely greyscale image; use the same band on all channels + data = np.ma.array([data[0], data[0], data[0]]) + + # normalize to 0..1 based on the range of the source type (only + # for int*s) + if not np.issubdtype(data.dtype, np.floating): + data = data.astype(np.float32) / np.iinfo(data.dtype).max + + return PixelCollection(data, pixels.bounds) + + +def preprocess(sources, resolution=None): + for idx, source in enumerate(sources): + # TODO make this configurable + # limit the number of sources used + if idx == 15: + return + + if "landsat8" in source.recipes: + for target_band, source_band in iter(source.band_info.items()): + band = BAND_MAPPING.get(target_band) + s = Source( + source.url.replace("{band}", str(source_band)), + source.name, + source.resolution, + source.band_info, + source.meta, + source.recipes, + source.acquired_at, + band, + source.priority, + ) + if target_band == "pan": + if min(resolution) < source.resolution * 2: + yield s + elif band is not None: + yield s + else: + yield source + + +def is_rgb(band): + return band <= 2 + + +def _reduce_landsat_windows(canvas, window): + from marblecutter.mosaic import paste + + _, pixels = window + + return paste(pixels, canvas) + + +def postprocess(windows): + from marblecutter import _nodata + + windows = list(filter(None, windows)) + + landsat_windows = filter(lambda sw: "LC08_" in sw[0].url, windows) + landsat_windows = dict( + [ + (k, list(v)) + for k, v in itertools.groupby( + landsat_windows, lambda x: x[0].url.split("/")[-2] + ) + ] + ) + + pan_bands = dict( + list( + map( + lambda sw: (sw[0].url.split("/")[-2], sw[1]), + filter(lambda sw: sw[0].band == 4, filter(None, windows)), + ) + ) + ) + + for source, pixels in windows: + if pixels is None or pixels.data is None: + continue + + if "landsat8" in source.recipes: + scene_id = source.url.split("/")[-2] + source = Source( + "/".join(source.url.split("/")[0:-1]), + source.name, + source.resolution, + source.band_info, + source.meta, + source.recipes, + source.acquired_at, + None, + source.priority, + source.coverage, + ) + + # pick out all bands for the same scene the first time it's seen + if scene_id in landsat_windows: + ws = filter( + lambda sw: is_rgb(sw[0].band), landsat_windows.pop(scene_id) + ) + canvas_data = np.ma.zeros( + (3,) + pixels.data.shape[1:], + dtype=pixels.data.dtype, + fill_value=_nodata(np.int16), + ) + canvas_data.mask = True + + canvas = PixelCollection(canvas_data, pixels.bounds) + pixels = reduce(_reduce_landsat_windows, ws, canvas) + + if scene_id in pan_bands: + pan = pan_bands[scene_id] + + pansharpened, _ = Brovey( + pixels.data, pan.data[0], 0.2, pan.data.dtype + ) + + yield source, PixelCollection(pansharpened, pixels.bounds) + else: + yield source, pixels + else: + yield source, pixels diff --git a/virtual/web.py b/virtual/web.py index 374c1f7..6f51985 100644 --- a/virtual/web.py +++ b/virtual/web.py @@ -26,7 +26,7 @@ IMAGE_TRANSFORMATION = Image() IMAGE_FORMAT = Optimal() -IMAGE_FORMAT_ND = Optimalnd() +IMAGE_FORMAT_ND = Optimal() @lru_cache() From a7d31c7b9f90b9b72eb0e35fdbad1f007f12862e Mon Sep 17 00:00:00 2001 From: dlindenbaum Date: Tue, 3 Jul 2018 18:16:23 -0400 Subject: [PATCH 7/7] updated README with nd example using 3 bands --- README.md | 31 +++++++++++++++++++++++++++++++ docs/nd_test.jpeg | Bin 0 -> 33957 bytes 2 files changed, 31 insertions(+) create mode 100644 docs/nd_test.jpeg diff --git a/README.md b/README.md index 76d0e1f..932ff2a 100644 --- a/README.md +++ b/README.md @@ -90,6 +90,37 @@ $ curl "http://localhost:8000/tiles/14/3851/6812@2x?url=https%3A%2F%2Fs3-us-west ![greyscale stretched](docs/greyscale_stretched.png) + +### `/ndtiles/{z}/{x}/{y}` - Tiles + +Calculated normalized diffference index between two bands. The most common example is NDVI + +tile = (Band2 - Band1)/(Band2 _ Band 1) + +#### Parameters + +* `url` - a URL to a valid COG. Required. +* `nodata` - a custom NODATA value. +* `band1` - Band1 for Normalized Difference +* `band2` - Band2 for Normalized Difference + +`@2x` can be added to the filename (after the `{y}` coordinate) to request +retina tiles. The map preview will detect support for retina displays and +request tiles accordingly. + +PNGs or JPEGs will be rendered depending on the presence of NODATA values in the +source image (surfaced as transparency in the output). + +#### Examples + +```bash +$ curl "http://localhost:8000/tiles/14/3851/6812@2x?url=https%3A%2F%2Fs3-us-west-2.amazonaws.com%2Fplanet-disaster-data%2Fhurricane-harvey%2FSkySat_Freeport_s03_20170831T162740Z3.tif" | imgcat +``` + +![RGB](docs/nd_test.jpeg) + + + ### `/tiles` - TileJSON #### Parameters diff --git a/docs/nd_test.jpeg b/docs/nd_test.jpeg new file mode 100644 index 0000000000000000000000000000000000000000..8ebdc85ca3e99c993a320991e235943e78d40fe3 GIT binary patch literal 33957 zcmbSyWmHtr`|hDT1q6o<>FzFRk&tc)iJ^1o5~Nc?8tLxtZjpu|hZF>y0fj*jFTek- zweF{T?|atyboM%HKYQ;tpZ)$@|F;Jat1GK111Kl}fb#SJ|26?d02>n%3ljqy3kwSe z2OAfki~t`G51)$k84(%H3p!ew7u3}B%sj00j9g6A)NG>cTzvdOLPB(`;?iOQQapk} z0{{653JwkqJ{~?L0Rg1|12u!d|FirX1W2&|Gbc0@MgWxr1&svd-zY%$G*1ka{}SMT z2ns41ItC^dHV!V{(||T&02Kud4HX>?0|OoXX>{<@cL1FPgOq_^9`l*bJ1j;IGJ){a zVr(XbrXh0O1vs;ywPyqlF2!?7suwITS=rb*goH&z#l$5PUnwc8sH&;I)zddHG%_}^ zv9)_|@8Ia<)s+oG&G%W#tu>Rn^Tct!?ccon76- zBco&E6O&WZi%ZKZt842UP}u&#;nDHQ_tP`P)%CC6H@C>UKmWn?ADsVF{yVV$0~g5? zE>v`MG<2;0;6g$5dwS7G&@mYJF-hfhu-9fo3@p=;S;xT1pp=s>h778UXhU5$nj|Ypsy|(N8c9bJd%5dE1K%L|P?2J@Vy|1*6N)=0DUwnq z>$f30V7UX>r+=c>V(XlFKVOV(d2IG3G!rCqSu^4}(@VG=;=|R*S1br*#|ub%`MR9~ zWQy(XKea1e8}~bBI1+x_yBM_Up?RZq8=}1I%>66(Gl@Mo#@2n+4-ebu(uyT4OR~9w zo}iPN=4{Luo1RV{nR}4@(WTXxZp5&v(y44&iLAZ3)=BVM@N#4Po7~V>)dQ+}*5Tdf zt^WY5@4Hq@QdIaPJ$H*SztwPi%k+)Ie5{Kus9nWwe^{8>xdhjq-AD-cjW2lkHV&qV zla!`qn2+Kgh?%7bM$x!)+7HQKt@gTADq(!$B~RLC9Y_AzKj6UJ5{)1cUUeHq3! zQOum=Vc%L}XiO3lMV(5{#hz#@`;O-wVVo5Be1mjjW-2~mS zJ`-*8_17{T^nnrdsz&s->mY|ycqjorDw|W@m zes*FrqJ+BNY%fQR$A0c(uAV^k9p6jMO;W#uAG^-#Aq;ia}-J?$3lJd^~m^ ztkk~R`rVDmveS(CGz|xOfZqibj1WwqWkWlekCDZCP_wPXRrUh~dTF4cy5F_sSAR-KlpU2Vh~Zl1srR%N9cW3d=fo=L z*8K@^w~-nmY?WGLWR}XepYJnInHC#<Np)Rroi>V=uCm&`Rzy4BDA0gingTk2J_?dPbCJB^_@oH5=%BSj_bN2LY+wujou zd)%b5L}utYtcoPQ5dWMn?|#_m?=?6YEQb~DFEDnb#;nVz-m`D7i&nqiWvP?#E_I6p z+UIryS*klhgj7{ZAJ4>VXg8U};(k`_m_5V#jfmlso zAuvB#QEje)F_mA{*EL+2k?mC}B9Aa~k9@Knna*WtUtrm~u0e}T_Z*8>**a-{$o2Yr ztg{m)eQg2HMqu=f^>$yjIAV+qn+cfAP(mk})=557rEtEYc40bx&nkY9xDQ zP;P-;So~qX;ECJ${CJF`Ly^mc zJw8MC(u4@X$nW5Ve}DremmMN3?f8!>&;BAk*x?mubuTwokv-_)1N)X3{RI>KfmBws z6q?t*Ea4gSVTq*d^yc|nDn=$hsLktIE>LQU{tPhUIIj4tBk+JVUW$9aklZw0t-hXg zcxQh{BObE9!JMlWQ}Xvjo?*oH0Qhu2+K45O;2dG!zu-8qY)M+BBXv->Uc1nL!9=f* zHETZc30jt#%u9tPF4uOzF};_UU{fw(*a)$fwizy9%zBeE#DuU&wLNAu>IP?WDEf2w zy9p~Kv!~EyN+%S(n39Fa+YHitkFn85qf4qGZbl;Q^N+@cd-&TT9m>+uazNXv+BBgLJ5XA>&NgR?IGAu{t=l|D`?hnDIz zZMb>f33h8GG@QewtokKsL{iXuASN~dD&A~H91iMIYM;?{8#vP9I-AIhzH!2_Us7|LYd%tYncokRrmjDQ`ANm=~ z4bznO2E(Wpc`Gp0aP&kL+n>*Jxr}dhE-VA+#`x9Xoz3;x1wr7^LA^SDN^@y3# zxE3bY>Df`*GwHGDT|qIjtb4`GKIFethPC2L^!OlNqkiK3`;w9ypueG|uAwItC85t1 zuY=)4sBKN1%Ret#?}Y@+%+*td!@_9R!!Whbn3FVv6ZZxPOp#i{?JxBBhRB*m%TJ4^ zNatJKf(qWK9!Q5#0)SNDf#^tA+Lac`X}rT6qlWZFT{l^9^}|=hxZ;zKzEZ5 zoZP=_w4}Ko>Qp7Sgd5+LZ;=z?Mrb$>YGuEhg8C1MiKT0-djADqM84^~ejM|@@uOqn znMwy`r3fY@DXa7WkctbEdP%eIJ2Z)<*p+VO7=~a{b!GV$oSUtdVFaMF{uS%gpzI#&$*}>7Oa^ z^oBpXZ~9gp-5{=IcGONzgPGJKNB+pg6FUAo2>s!~@oysB9i*Fa({kaENt*KF znTCJJj%kZf!t|mY%a9IZmA2gAVUpg;00l%;1<{^QQ&2-4Hhd5342#NS8uv}Z9zsJW zmf!mgB%>n=M_l}1^BZk-vTt+tH$n^s+6+34A81jEfe&GQ98QU-b6N)rwS-Q-W zefg+*5dEkhGbg30QP?8nWLhzkao1QmXbk^22BI+>vA&;}s$ev+U^LpM3>wez)|}mE z!Pj5n(5S?_|Ar~~>yf2pRqEYdywL+e{rW$^eh1puS{K-_B|-t%&obnxiZl}s7Us7^ zYdfGF^;zn?Sc~y7e>;7pbTQi@V{ox4+T_wT_*F1r&(IF8cy7O6pL z8bl?F@cqSWx%&?|mHLG8n;BtpCem&w=(y{Q$MtqNT@Gz@G(oNXHJL5G+V??JU>*Z~ zQBmrTv%sv|3o0JQ9I&!OmREjJbChChiij(ro^3&$-g%owPa-bo*zzxpzFQ9!nsY*f z9SBscki(?or|koIXA>V{Vo55yV+QwV zW4Y`}(~Y6+?F_*3HQwAyFOygp8Jh|zi>7`yj;Eq(?{1lDI|Fn_?T-l7|abRM@S#Z(Pxmhc!i;(qLD_t_CVuVK_?h`G)rqMtAI4TS*~rXQ{UA{K>8l#5Rn?0 zc29BF&gRj2EvQG+j`oXlRVxpO12lj*OT2Sl>GX)Q%_ceuUbrx z4F&cn4R@Icm(t|PGw}C5?GNPqtVblwC(uRhs_R>jp9H)!%jNUKLbvW)oVh0&e}NxD zmm|*bjdw@2mB0Y82u_r-vSFO+WWiyZE$X*mj3W$y)VGywe}A!ZnI4ks>cZks2!ySF zpGm>d<&&SFkK764=s9tUOj!AzWowlb(Xz&~b-=n>U&&SMXR}GIjjnRLL=BhbSOlKa z@aufBR!mx9Kh#F;P%jgDFJMar4S_$ldE{Z-57%hAS64`#v;p!QZggxXP4#PUT4cnM zhBa~eCk|+2*sLjrL%A$faH~rjdq|pX^~`5du+4u&P#QkRFE#|=+tSsVGJ!Sb*3m?_p8$8ytpJ1^oobK>GSP8glI?+EU-c=J7Xph$AG3P+RwoY3 zZV^&rSi$}jnGf-W^;i!63yg?hMD}li3!z2Z&{StGSt@el;axNQ=s*V?m1HfSYK)AyO+UNXakJ zdjF-`ffCIQ!g~BZeOs#{!+wE|PTU6Jzeu-Ou!n zw9rx@WZjo#MA5ynmEgCXKCK=8s2XhYOek072j7~8lZcqz!Gpb2H9vN79(-6o{_d=J zoXPdhV#n7Mh{9Fvil@Fn@f^j!Fu?F*!`pgzgG`3u&Eq)b5cDlNt!3dh)%wovm^r)u zVE+_U%8u(^dKgxU(iw*D^S4Cp2TUI6B>T6zu`M(X>itNBps7!_W68hl3$Rd);8`ar z3V+FbfN!3ujSXBGU5H`Lt{Rsv?zMocsOSrf(Upa9%U1b}L_F2jQP8cN(Vul5xD5x; z`~#AEEh|nth?v-o7Wehc%Q+oV-&X3PU{@EVMvQr_38DP31f*Va3)M6{`&!uBH`mr6 z&1i798U4s=W)-HAl{#YBa>2Uu8J)yra{{AR4WFM)^&nB9nvD`b83a7yW)Y$vYafhZ zKk7&y)k@mCHu{RJ>P5Ochv<{4CABfz`L!*LP0vh%=MoqFGuQZq1XVig53LmOlT^@u zRWXFST$F^rP-v=D5>Wj*ijrg;%N0Dud|Tmn5(nJ5?T^R%6hU&hc_>*s;149BqvbC> z5!H0rfPY`^;TJ0w2D+vS#4z z1dXL~2Qt*J!ra^>U)L~NL&KM$Y(86H&qt@p7Su_j>_2puMn>d0L8lJf{ve7gpgi?TM!4CdYhq=u`}1v<(n2ZHtFCEU6b891zV$s4 zIok00AoJy)QfehiReu}gJHP!@UEXWyvX8rchpES)WxjB(FZt>2BTo@&(9X$N2*XO| zA9f!ad+o*H7U{Z+UrVJ9MjS#P$kA&Ah*g0`3-Pgl^L{3z9YW#qSeAI>A7JH#smt8} z)lWQ502HON*FjD$n=khSP@7njn8rhKg7jhBQm^O!R+Nr$)z>-p^@fMNPE$!v_bi&C zLbJl?zm$_m-ndrl=vw(Ltu{%zLHWUEDzT-{?(Ma}bBn^knSwiN!&b0`z@oP$w${Fv zc>fZ^=JmRO#()6!yzhaF~($C0U~`%mgXTOaf4H<&2&R_>f> z4Mjeo8;?It5j@4!n){U8X-{#%QL4-9Xz9hvM3vc|=F&9C5G5#lGV+nC=COSDq5L0k z0|h3K8I=%+AcgLw-?I2K&No-#UOr`(jGzzcr2dH@#mJY`1x0D9Typ3Evcu+{b0O4o zTsdP$7a36t?Xi~^uRIwqlATp6ubv+b>ZFY}YK!S;BD>&c_v7xymJY|umK)6R%$|vL z?jKOA>6ZD#2*U5j#c3W)Ro5Qg!8=SKBAQNv;%UPOqamIU^Lu%eF8Nrp>kcX=`mjKY zV)(qK)-&***wu=bg*Hx=bi?Kx;%{B=7_JTM@AOkAvc*`i&+9SX|JWgg-zKzb5ju82 z2cW%jt9UYnmD#o(5+Cp)bn-`z;!-AlN@VhOpjV|2T2j>tH^PyHI!#MIO=)%GeIsZ^ z)9+VPhxT}PaHa;g4T3zr{NOu#VMgrOYzfiVXDU4RH?0#Nt(eY1PY5cf-{Z~v&|HHf zfSQ7G2rV4Gkl8dS)M`O|P08|gX}*-U5M7FUcI!*A#f)xy=#f~FGMqC&=#Kdlgk()B zM2P$Y&N-KVwt{zwsnV|BECy*}Z5&?7wqDpR`WSA|yH}6%QISdT9FbJEfMp8H!Ix~gU&Rw!AVee5D2Aj~N9_A8w74$)F9bJbg_A(P9dyePI*DmR5cxuyTTwSEK0`k1zRDTlcWeilWy_lt~_ znAcebqBxANqZHo@GfVpHVVtyk3ZQ(Q+bg8#N$jUwg=1-*-whyhPSS+^Ku6jGChL8R zZG}WyU@!lQT}Q91Ywog~=iwa;SZu0T^cmtb0_=Y*25KEv$9FYXh4I*oJ6e4^bYGP< zqA;_giJ=v~lOq!`y;CCpWS=5m`8dxOZkNVN($WfM&?VoH<3jnaieH+v`LWJh7M_F5 z{s+u1)SWEdtGhjN)M_4Fgaq?&^B*nk3{&iqmpF<`G!e!b5loA zCsNFEy>0%c=Vjm^&#Ren3Z)c2cabWT@F9Js6n=$OM;9K(&W=pg_pyDstVTp;T>0_W8-W?>5)M3?0ryOt zJWoO@$X4V8uh3>6=xw3K@}1wF`5XchkVx4c9Hd#9j}^A0`ra7o@Ll_dAJldIvK2up3f&X90M{5> zGsR!<6@8G$QXA%wDPpVZ1DXra^)6bir0&HiOB@#8JWo`2$F;xl%$w9<$>EHn%(uY%9N z$0@~r6%plya6lx-0|w&Ud{!5D?3*(8pe7@({{W{{MkN$W>%%YNTkGcgWFW2aXcVC`|`h2>r?-=?A zkUpsF9k_rUy>1d;QY9o;%svNW-Va$~q%T_?_%TzbK-LCPKGumfa;JD=0D)#VmP9f< zY~9gkviOjg?s0_LQMea#0uF+rI7|FgkE|Jeaa`r*^87)qwtt>cb&rfCR)zVFFIpK&5`FR-ErRHa z`GY#rj-tcX`4?4N->W+X(4_ZKrf;fyrRc1IL5hnRt=1d0s=a3TTCy-5k<%mLa&7IU z5Y2UcJa?zQINMQB%hlL5Kc=n9T{f|rQvI_rs0K1{KRJR~1Uf3=Dcq5Yak&{w|Ek;O zusGfqUaePs>5H;#Uj)$h?eT&vAkhc93qJSgOTCKa9X~IYVYUhoTiww;cbq9Y0b=(- z!jBg4$BcbV;noZ5Ii5O=kI>O7d;)KVtp2FW z6=-B`!b1|vrYD(#7Uh*8Vs1%d9zJWhMOGWvehtIFZKW{~Azouflf@GNBX$1)f$MFTxvZabDAm6^2T`&rc zEHDYuQ^*w)Y^YDV4wC%+kf8{y5f_RQ+T=<_f&R&(rc}GnIB)a(xs32_pDpG26%p>A zj5a$2lm_5i@F0A{B{e>SMtu#|-pxpiQ!#UijPsqFvS1s~pG05s)CHLrZ5ho39vHh3 zGu0?Pk1FjnPEy-uO&P5Mg;8=44h`M4?5!4(Ja`muZ4KsJmZtyQ5b{f4b=Xitx3(H1 zRGXwcomXF3lHqZp5S{XdVus&DdxLvAxlN#FKM3nL%Mw3w&7@ytuvVvz0Ra?s`+vc2 z5DKTO{{Uidw+oQJh4`-vL)FabO>{}Sa>-<%$czskHLr;^gb{oN$x-tWKSru?kH}21 zEWgVUH;bbD2bdMdj04hQ6Hbw%G9qx^QHTp>pJU~BlU7H@x(D2r3E~lQL#E?Y!N|n0 z2vu~{(Q6#|{;olR@IsIaEH}I29Zh3T{)H>EweQ}bUZZ3)Ff}>!Q;ue}qul9{`fAf? zp`Xln6rskJW~yR#zL}vfE7?Ras4p!(#DR+pxgUG6{W*T5KoXu*uoyZitGraV@GK%| zH%;|n`TP8>sW@wRq}hDZBlVetWvusiZ{xN;CpKS&hQRVUFPa1=;xAd)Q(XSd8JUvxYXa-%bHY4OwcjL~zwc;(%7p&kn4L0QefZSY`W!sm&>)DwU)+U;R zxIAow@{K;Exiw8xK|fO5V| z!etL{E03?;M8=PWuCEy;b)+Wh+IGK6MN@97jRfdbVB`C5a4y?xvy8gnr5!Sa$Y`AB zUK+<>(J6PCxffxVJ9dp}BIx7lS85g&#|(;CL>+XT z6v|kRMYtw^*AihXE?Hww)$-K=JM!L^SI+0It~eoBpK5g#LFtW0DeP#VaB43D3BPFJ z8m4_Cp+3vuw;X1gv}MIxx!;*s$O=WLFpB8KWh|ILo`UNEtg#QAV0#z=r)SOA3B)#&K&>{{vj`FZx|hv3-^O*7gnG7!A?see(Pvc5Kt13eL#gA8*&xOh|vzg z2B?vN7HQG~oBCeWqjBjbyAy|Icj0gL44onlz^>TX=)ROyK3dogMDD_>F$+EvpF&4R zw0t_!RQ$CsLPTczh7+4PZ9C9jHc6?JV=Xe!3bqG!wk``Q?b8T)dztcCym!btPB?MA zlNw|3))N~3QeS#@V?w2L;b2XN9D7;d-o7Z_cj@A>RK`%VC; zBUqD*n1%`;XOq2}Hm`z@FSnewO-wf(PK;Y3)escvsEW?#H;d32 zTi-+fCqq9YC~y~5ivJeWGW5$W{vUw%JKt~g+m8_Hru( zSXl*iM*RfvkMUBV0}V|9uhQ9^wscURK(M;7aY|_+bBzAyxRelE6*`FPkP|Q{!?N5f zd81Ls-m3xY!C@?Lk%9UruGEe->f5MM2$MVfm@DjrqDjOTSQQ zQCO88F)?E7-SgQrcl1T$htd~Pam+3oI6m}sH8U-Y7!q@ojbmq1wHRU1P< z`rw{?kV)Sn0JUyeIfv?b&U#X)`ZZVAi`lVOurkNh!LNU_msU9nV+oTc?n^`mQRs)GRZ03)ls!$7Tx7Mey%CvZM5ZJveU`$wh-e z68IF4frHU~LGLvez`YjNHbAtw>6jV%T)+>(NVmOf$9REy*tF}cfclcxCUU$J?tb%g%8 z{eZ65Ol%Kjt%8OP?Py!(V9091#jErzbu+jc6l*I=Vy?lc(kX-&^$Lnw<2s{OE=(T=T}Fjn!Llg5;z9OV zuX}r!%7ec!Gi%~Q24y*ePUvYHhGZhiCw3p_1GWq~x?Ni;w&XjwZgJs!&IPEIK?yc_ zw^AOH=rU;VJ9w@JZ%i*|$&YU7zr)ZCCyFcRw(#h&JhkKV&uk!~pE$$Om+QE`6Z`ox zS}12r5ZxuSWtWMK#RJ9O$8Y?U=!vmkG9oJRm{d|A`gRJc1WbgiW^k2QDRGcY!8Z}>ZNcW zwe+u5K_@=-Z#%(p;s$Joun228N-A60~rGipTZ4R`rLc)KLDP7Y7lWS zLOfJ03!QCQw^UVZbyhx#)*XCFc<^xx^PC5@npVOx>5(K;)iz6YWg8ZZ5HuZD7g_Am zwcWdTPel_;bM_dEs2F~czqPPszPnDvBLDp618$`ND$GQBNROoRGxqrPS4=WAlWDQH zZ+p`f^cIC{K8pH5rSIrY>Ad}mb=Dc&KKU>FG;6^!^clYLd!%iX__K+$4wf7#VB{cp zjPj&A56H762iS(yh?TOBhzoO%N%}ob(Nvds+lY{=&bxnQ4z%u>O`>{p{{ik30iUq= z3skMt^X`*4lMoFb->tQRq$NjThk_PuWtlZJmCsC6N13@fitTJu>Z!uzvl5&Yrvb>D z3Ro=t09f8lCbuT9@I1aU&;KK{HS~33hMK{kHO`&>>$#~D+XxMRyN=G-{eqY5%%1E| zgIT=#B0Z_u)ReV*DyyH@hH|JjrWz0yFmZlgyrous*3ZPWe@cz|_au~D(rcgrvu_6v;$oGwv%?&XoxK9aJe$Uxn+#Z^Seln?0_0jiIjOq{lY%{<$ip6GJ+=(xb$GQ%Lo9!?P@Mu-8TmYT~!HNCxrh)!o7HEM#Hu0N0A;JxhD z&r#LjXHEuc{%*k>0?=|;eiaokdFI=0cEPbHWVTjzT7UuPN#(k=phlQq{5{DyRI_)i zG@4$T;N-HSV%7n%^OO(_NqS_Cxgl$YcJZZTZNOEKv|&7N5H~t}g4qHL%jL?18G4M6 z`X;}@!~h_Dk8QAEjS*iDrz^;pe}G{-{&c*Mwn_O_SCa7Es7IycCl35uzy-(h@u{|% zFMD1wD0sJHnRU}}F(oy(KK|gm(Ia)>x^*)_ews;#^Z3F^G<9PgG7X;jEZB@KENM&w z#}U*1ECv(~mks>Q|8?7KuoNWk<>uCLtWjf^-(h!@`Rk zqT?6)ioLS*pyTmH#UY2Vs@5^Y;Iu-ok5owtro6}ijU#Ig8;XCZK(@oGw;;S@nX8B$ zE!c$trES)f?oAq6Zw5dj(p39=~NakN>z55X@Ti zU=Ft!phPp&z_}4@l}MuJK(onWea9LBSi1)&iT)-cLVAr+{O}1bcl;AK?@B6J@}h?_ zn^Z*gskq8aE`7y$&=db7zQ+Gu7ekI|q!07Y7gF;Kl+4r0S46C?=k|H_IX&e+vLwv= zzAs+>CXm?oj?$^+E!YvzewU+k*yXA^!lpZcM#l#Hp?f4NjDoFb{!RQ^WgQxg{az)3 zpvTyQQQYob{`6=BWdxsx*~JrQ_1lJJI8ks^R0RG3CG%{<-g-mozn?p4|C1zzwf)f8 zy5D%FftDb&ip|jvGj)#D{a|*q85z|~<;jUzqi5Higg9s314$CwbM=KFq&khFcF(o1 zO2ux_&&IS+cM6p^vc3L8Oi5v64NtUaor0vDgzEss+tveysnw}3Oo|lcXaU5QHQZ7&aOZ0(jj`#J-rV`?nR`FTf zStZMC<~1ijE@GC9AREYZ&t{IiHKhs8_*-7$kmYNry|t;9pDs`D?MzBY>FzF17X?7D zMzrD$cF5sB9R3V{wQkEO@3?wrD4=|D92;VXJG42wH5P9EV`$T26ZQ{yek#MhXtOM2 zoVXxknd@8|0S=z)xQ#P%8K=CDlJxT^;2jUe_~qg&=vB-tp@f{`q_&^TJsPf(9}Nn> zX1g4%A?1hpNw~!Zq12n2X}NfnI}x`dnZ-t^@y|UnJfcyeRSzSV34@w2F+tJ|xCIG>k*K-9OjtM6knIV}8eTD@+|aNJ)lL+I6Mk_~}LU z`{Fl!pxk_Fl}M&H_yUm(xM@kWA`w)FbV+9^tYeI6WDIk0t$EpOg zd*IoViN1}!d`fjzE`-H4ETq+k9&$g}aYk>6S)&4YXszB8CcrTZygzn7QC!LtX0A-Y zX|%GRLw}mp)|V2&5~o@w-xun~Dv9ZCRzqUJQ1382uJ*zA+4`38Lr#7|vztgW1#BrH zn^FM7+9*Qqy<}?PX}*R*^U@_(FMH*WE4+N&9f!2vmtW$}RVPR;Ug-m=$1d(gn29P} zBU1e>MVbdOQawMs7oa0X5#+rzD{&Dk7ql+@to>dmN zW3X)@=jp9n(^aGSOgTyWB;_U4D?yG9p(<=;^$(}%^dph=R_%Fu7mn4u`e6O zTyEaF5mXNmeA1y@MUMx0w7y)##x24Gxd$YWB)rIZ9^>8c%rR1D;V0?Ya+l^=nR??H zl@V%1@z4_^qQEXpgZ5@+Zum@OVO!okg)2ohzOA!)4e#CIh0T$)l4P&T(~gSFUO$aps;zOoMm8B=r1F& zu4Fmkv|TDqJfF6W^m6RmJi!{j0zr6v{{B4CyiK50`q%H)9fdOnX&ZT!*JPP*qv!>J z5qAYXugGg#gFhW30~DG*u)n&iH(;vc_7pi)d!t(_QX>$dHizqbplE_9V-g4*H>Vq= zpcT(KF&=naG`5aC=!~OC64ROa5bwy$%0Cs}d5AiiX3{q?3^l)v?vTZZ2mi5=EV%sKK)6Z@YHGX z?0$3xlR@O9&w~~dvAuzu^dV4IyHJU-yRvV~WNKovee0lCFDl9r9%ln~o!U_a7UBPZ zaOWVG{a`;0gl}KHW`o03JY=Il?BRzJNN0t{xa)rOmyAeA-V8MnCQa3 zo64w9k zp+Dci%L^-C=K3JOqfOgA-NaBhqli{-GoN`=!Mdjagi2UtjE?f{h?YO>;<;9L>X)!P ziNvo(tPG{PryQBrQp?cE11~~CY{~`{s9=(&?=J2GgWUn$i@%ZDsbd*^J?0wJK~xdV z^>f>LqzBAiVRA#_00ozfAmrq&S3P?w^cvO& z9wIucXZxgL|A0cy&__}8r%VnT>i8lkAUeHf)8uqD<(c2h8qQ1Z-ju;F)M=FT+9j9u zq7=d-bJY-gg?5TRBH1}(SG{3p(9%~b5mYK!;LdUYI#jyll~=h35GiEscZ{0KVj6>q zU=WN}NYa6BcWnEXD=NIeJfIT zbuQi)oK|DC2Xj5Xcg@DB51r$_dqjdv;v8=<7j%Y-hzi??HcoVt)@ejXMLGmCfx59C zxng&J{IV(W3Cc_Vg90KGq~(D$#XGVbmhI`5?$UM#SMrDy(V)*9gz}(>^wXhTKehJ+H3!;+QPE5qM{mHsss97mZHmhfMX~kC!GKMQXCTD1(tOx0e`MB%gvV zw{WX3aWciwy2xPx$zu&MP?y8Zj)45Jx?!6j!BhD7X=-W)|Dn75g9B<|bQ`Getx6aC zi;9iA>leF4lS{F6=-1iWAE2qnhBAI6=N&=-|J{+uTyXMs>HQKHExXN&{D{idXh{03U#G#O62n>) ze+5X7P#ubp{UH4QTbrh-X1dnO1$VYeXeHGDIe{nHXLMNfNo+)OZIM+T7B~HSJK$34 z8RGMy$SN&=EcaplqL{F%!r%P)%-UC_{J006uA9T3`KEI<1KN%XKJQDDYM55S`j*_R z(Ze>=UWqYeeB4t+n?R?{xQa*0j(5NL$WzA?6^QW(AtjghcaWTKhOg^gWZg?vACW+>1Jm$0a5cg~|XYM1qnC9VUGL0ZVC%lJ6& zm%V-STS3iCDV%98S3@5>Lj$CT)YG>Bf2e6Ck{^KzjyHhYoiF0SlmBLUNKo16;}o)m zN4EK9Qn@M{azHJoNhe7lR%ODm)pYeufT=aV6FU1Kd<-aqc5*y{Sxk+Yd>o% zRDFM>M*dkEw*Z+YYBD^oEM8ggy;|y@{3Ja5s;0}&KcQ1QAOrn`dgBN(!k^TdJ2LuUx;QX#OmW7x;bQeb+8X#OPJ1_Z+*+T%%p2M6^z$L?yPaoB`-=N~ry4O^QH?tQ3bUgWlgfJmOtv62+q?FOh#{a?UCC-a zi**MSe!g+(Yyg`T6MNGuXMJ3IR%U`yTYW@U7{xz8rQlwnwY?n?q-q!aINFx0#%?E& zSL5TegL$>SJ>#kA0#WD1G2cBiIr-Q>#d36CH#QQ3?ls=bxJS-WU-@X6EgHN(7A;SH z6e88G_3Xi*sd=SW5nSb9faBrshml`HcmCbx*yxZ=ZN^`{^m2I234EFtO&iaT*C*#XgVyC|vV>QW6)=i?nL@%`*aZ)EaRK|g z%)t1})H?SUMz2Gxe3b;bWTjg;@rHlgElSDg6o z{7A8f4<|Yg+rdkM>2lN<9S>&7AKn^mU02OHt#+t+-IvGx;kn}Q%udz6Ik-&CI3d3( zO`UTQ%Dg3!5*$5ZPEvZJ)Ak~W=p&*1dF|EQj1P2%*9so2;mkh}gu~_(n_sg=L2;X}0GR-rS!XMtRG`UPOLB@{8!PPuuotVL z$jmalRRi0Ag;!=jZ{U+;`3_DEC%pk!;SIL%!ogn}S{Ife^WY>(KQ)##^&$>-EbrzB zIbAecUV36bw3Pn_$v`&0Jr6>>@OOqZ88r_NTH4!aF+8h#sh!HF%Ongt0`~;F?OnNF zyaQw^um_lUci=_*w(o!9yQi_Yjak~#$$a~Q0B-6S1HzCLIbd;sYt;M*x(1tKkXdLy z@QyX?uC}(f#ZvO|BO967NiKG?NB|_^Ks*KlsA=C6V%BxP_(rU4CWl#)*5mB=8ik|V z!EtdGRBewe;xfr8JBTagmF2og@ViyiJJIIgX4Ww*7t2!I9L1hB_}l8Et- zn{#eFWR#!}C>Q$A#XTP1L)Nw18x4C*zG#};Q6kh``JXCbSFFg&>9io?LC$lY&%l2K zbfwUI{{R?RNj9Z3AM|gu67U5vJPwOCC9U)PdgJZC(dO zB1xEt)wj7jNm7Uv0+Mmb+`+olg;MF}22IY2gco9fT#>h>~a}Qj$j^hF4NBK?i~e?LZua{u{p4e0Qbzn%i8G z=S3Fb;VRFWl<+eWjucXn_Y&LdJelxjZjQFLS|z=u)5&Y3>8TW(+vP@*>7tX% zNdQ&~0A?RDwgFab{U_p1c=%MZ6B13xrUEnToM5Q7O$A#x+`-%-IEI` zB7Z*JV`%|a-Y~JZBo)pfu#oE1{o%e?Io2l>gS>pR7keh=zwk1QmTP_686=HlQa38F`CN4ZfISWm8fl&u@r)iI z@U6zpD(lHjUrrJES#mcyiWY7ZWwznlW2-K9GH?0|QPcFl8+b{rrFkUNbtn>4k}_j$ zrw)!XG6*|P2v7vv%?> zOe_1?m@)}fTXK-fPaNRnbq60scvHt3PP5@e@m2PtBr6lGt=N+>!EG~nN>l;CLfcq7 z7Ahn;WgCE(;wOS4@z00Fpp_uLw{tD5?tp;U!A3a9M&mf{6yWW_@ZW~o>^=nXUDk`P zN1HuLR@3zxc~&D7EUXn>oRS$qX;p^`6zJi_`duN--WtKd~vC7}PV;SD#7}^M7l`I2E zrQ3KH;iZnTsOh?;sDQjT3R_fp33(DvCL)okK33K&q=^Xyn1ki1pbS45cs@-gYdc#h z;CrRtABShWM&Ai2;Mtk?02mipE!nwNjh7>WeB5}F`@)_iHeMCKp4J^YDJkY^Ro=sK3Nm@l6rRbbH7^`^!6MXrNXg*;01sK-qr&sV=_*VZ?I@sm(n=A@%7sI)pOk}M*Ja|n zUlaIWQq=q*GHNqGlE(tc8^*10s2Nc2+Hsa%atA#{eE$IA4K5$-`Lx--BU)U=rzg>{W0-%P!=k5$$+8SXyU1>47KXKiW}u{-v*ViALo!I#oJ{kD4{uXC}^%!-uf$lunt#09pO}n9x62B!vvKx;cWE-mCAG$WW)9@-0!B#-5C9HB{o3}FxwFyLwEqAOX*1l&b8|G7 zum?$5Wk{u0ncWqqIDD3J!(lCw0ALOQeg{FN_`vw4Epy7Xx4-_;TMK= zP5|0Auh4d{QMJ`2z1A8CEs{+u>``3Ds4T>+%_BGpfG*YnQ*mwI#=LWy=SPq3E_@}S zYW5l?>)l%{x^30ihE<7GjFO<*&oY1k0T8n*0vDWDOKq&#&la8MKiKYd2)x)V?#PaF z@`h#DPIo9jFd4-F62Wb%$9-pYBo`2gRZMXbuK@~<12Alf0U_0f(lUf{0~T6b>eqMo zB7HtP8#pdDMDoZblI}f+EEP^78;anP#DH?DNXD738~Bg$mhV8+ZahDuY42~T#LSai zU8*vomC7_hSb|7k!*@&oROEmG>z9|-dS$Gach^u_+a@>5Jb@8%KpA7e$sm3;Lq_o> zhl#8q)HExYFD#l?j#)P{Goq{H?d;ZQ1A#&YJGE2aEn6+v`3j)+Swg((2;Y z$s)Hbw`dk8A(Rl2urjV(V=6JV$QaT17vib3t$uAL>`0n!pvxM&XUmG~3{Yf{pvbE0 zj22>u2LR9q(g5eJ9qTeb5a};}Et$QJ8>ruFL2?W(a^V9ixCM~A1rE7ZCjgUL7ZEMR zys)`&k|9Ph$jHDEju}WQeZd5Fpa=A>B|a5+TmBNOuMkgbCDgYr@-?l?2`$8r_n&DQ zB%Q|ukP93Eg~R+?@gqm@#+i3?qD?|ecQCSB5wT|Ed2J#2LxYT<-5^uDVNgi<+s5A& zyhGzy;HrE z&l3p=ZSr8j5g>!LQhs-aDBO=B*1TV->ppyK*VfZ~1+>COA*t=x|H&rB2Xr zahz9qrGC#k*Y=&d>bjPimK(D>TFCDyQP~3^89g$3XE->)`Y0!JN9#ZwHlOf=!1`^o zM!F(gv)h|_CApR4W&{kd$-+6tcOX(Saf;wPTkx_CI?~p{D~M%^Yxx^ZF6IdkG2Aru zk=eoCH?lh}1ET|zHG9Yz>0HN;JRfi4UkhC6QvA0onXKe4R3za@QrHJ8$pm2VPC3B< zA1Zi{!;$Fm>GyV8j2ecYe+YOdjHq>RsLvFrtlK1A>yeZ9hdn^S1fQ64!?ErwsnxXE zb!{GfV91H9>8<4I*0Pe94tMQaf~*Ek=2)Bzw&e`jG7>p_?;F6w4x)fQnfOid5z|@tFH{oI6R%eZ_8z|oAI~63#i%4;=2=x?V?q@g(At30l<_M z8Qh@!#PP9ujkrDurfRpkCZh(YrQAKW%q|%vPN?6=f;b@XK_KK;*IoO$vCh!L;aqZf~Br_H?@L>o7QotS1uIrHy} zegPg9)&`?wbEaxCyF{xN(QN%pTn(X4NkPHFk&l-q@4~%OZF0qYFebEvSkj!|czWh# zE3t#_o(5%L%#I1lfPPZs=W2ifLT`m$FYx}0;!DeW&-jx^Y9q-|{*g-mV0W{GzYRh5}Q4I;a8Nhj96CU~0b zRrozU?VZYLQ;XlRYIaR>jVNU-lH02;cOBSnF}g5MLP7I;e+udvw}I^7lJ8ErzlPl1 zX!iC~Lf02I?03r>E<}t`jHN&fLq~9@s!#{i*u{A+ygn3Zo-Ml4uk8HIYX0l}HYK%| z;n}5)*kF}nyNrrQkV9@M*c^J`8T>1WmTf~&yuTVop?wwY#@3%l))Gs2mM{`Ie$O%} zbzlJfV3g6j_M^ zp)z@D6m1>QD&>`GhsB=**!U~O5L;Nj6?Sk*BO@f{zM0m% zGo|bInoZWB2bnIHb}sFha6^D_3G25zSb#=KE^-Dd&%P`8P7fB@Sm}Nnn!-(QMHpQs zJD8(nw{U>rINOXxjFlTn>Pm_Lt>Rr)4~gFjY;kX!$>Nqs%o%o%QGpBj0wrR9c^I6j zW5Vqu7G)Hl4gMor>OK;<({1dpqMJ(6q7r?YQ*(7JilP-Z&=p~~2#iU>t3S*LCC+!k zsdTRey~l_AO|08Cr>S1TA<|GAmA08ADy5k>MzWU+#z9?(N$fel7He1Y_@_wIyc2O_ zq3N3RR+CzdEf(3!J2ZH>c!qZ?nG~ocMh6>#pbuWZ@m`6g>lRAtntrjXUcqU1WeVFp z_VUZ+G{t0LRiS+sH?_^p30O1?HJA)#uYL^qZX(g3W zR53eqk(K#Du6M>i2<NnqS&rc}r`O{iHzmKQu}k0BsC&B#p6Gb{~+;wb^NU z+g^CH?C{MU=9_nAEzPE(71GX=h~@#11uC(a%)5A2E;jCbp#%dj#@8Chg0zci{2Sr# zwtcHr`$nT3v~MKRZet{e0bB3jsdM)MPzO!CcgMdCjY~qX@r9B&dsSGj1c=DqOPrx9 z7oi?nFu^9jnRQg}fg!ZHCY&&eBc(|jGGF}+Ldrv8rK_W9NA|1$}vyHn;uT}u#_`}6lo;BAlHEWwoOGRcTHI5VvtA$V* zzzre7j#%-|0)vXfeNItowyNaBxXLPp3_)_N*dH?wF~DPz0l+xofIWZ5cmDtk@4RJm z;u&uz)+G^n66toUBxyYA3aqU;APDfPzq`*xVYp--ZtvqiieU0>bUESDwCiY=3u{PX zS%OBqCg+edGX=mWfN`7*n&Qkx=R%THiP&2*yweqqHI1CI`C$F&NeYakBW~Gngl>A> zYdAbBad7M5%e`LfO=z0(T~5|mt!>K@uo#BdA&3Q9Mo`Mc7=h-18G7C<_=hZ-_OCaa zZDkOKeM)%N7$88v?IlzOU}etL4Ci+M@GI0jC-6r})O39}P;E|0V7a&RCY(S+Nx5WF zu_UMjXag$=8C{zJm&)qg<-dS#JZ-Dl>iVCGV3$LE@L4c`Lv?j+Y`9@-Xn}&LA%IcH z`^7SFeWSwuD289L=(gS&)-?NgpX_jnV-78&j@lVl&moVNO2`#KEJ@CF_`{|6J6XKZ-qw9bQH#rT ziX`%*j54X%@)AjG64_M78CR$#yeq~&KiB+f$hxh}V%jI)9i%C;ByFG;#z|w>mM0sV z7zY#q^=Q{^2Dthx@+H0DReN~lo*>sMK_te3%i(sgC6shH$OP9rckxTZk=b0_X`gS_ zzRPtLr0J#Fw547`%7hl$tfiA^Do#r-KJc%b1>N6;^cn5i%6og=des6>^+s0zKmLD(6^5kV;uDdY2u5dA)TOc0S zFUMaCYFe^s_u6I2(s$d=mx||lNUa=$d|(0sI62zQoty>@UH$Kft~3t=$>Iya{jOuh z%ZmaSc@jzw?$`hwV~x0C7jR9-Zb=^_G~!t$)vc@&&_g^SVTJ=WvSnD|jIQL8Qz2V$ zVZWeU9=g-we~7#>tvq^Vx@q>YDnj}+%I^eDiLvmkccCx2$OIC{<}7kZ2d#L{D_6I) z_?h9Dt#ud@22Dx@SKFytI0(VO<0vzb>&D-k{a9M?UlV9n`b5@u+GO+GX&Q`1dnu=g zqu?=DjriJK7+~!rAC-m!0bfOYBk_Kf;ypI{C|XN>S*{-P>c&ZnqBG-r!Ml+gK3J|! z)^a)(Cy?=Hg{^!itz2oBt-r!vC6UBgK(XmJryC{y6$+r70LzpN4&a|K2ZLuEuiasS zE3xojjqbc9;&^o{dx;tgrVXiDrMe_;o7*cIDGQZTl|Lp4D~-9`h2k9w^TQhL-jAp( zV&>U`N#(ZXAb`ri`V}DagU=uyDwUnj*l!Y6f(g}gnCDH%;DQJ}#z8!S7p1u7)>(T%W!SD5uJk$MkfRyK7%w1aARc*he;3402;6w0-U#ia zfkYlE)2-wIuIx{PQ2Hq!P8&AE2A zxpk5(ZysNoravr~3<@(6fE7M>;oWlo0Ki@^gH*SH_dx9x&7`HG+^2_)`-7f&2MnN) zafwcIioP<@ zwOmtkANM z!=_Xa7!OiO(eN*Zbxj*n)i1T%4LvnGq>RgdrcWS>@bVbR3Brj%UB$wHeAA44*#NWP z?-Jd3p6gP%O;1X@(;=5_xUsaH5R4V$-z$)&K!@jzSLGyRj+M8sLu+@Y8ymm1s$MvO zB(;U)cw{`JfU38aRe)h3z)}Io;;Lz$4AAti6Y4%9w6yZ|JwqaC0}!D|O~4?+^2ivH z3l?VQoK{bdb&V;08eBfDsKI9*t83&zJQB++CTEPY0P!y7jN`7{GVn$Y2m~50jc4%( zhoHB-(BQw+Z#P@b1%#sN$wtGt&o~As)ZmjCX&fu33|ke=>&L`D4m>*7@Mu~Do*G+m zZ+)a&H^>?1+vYPoVWlP?E+R?gn{iSGE2r>9``COcA7|3g>KZM)_S5OMvMak>!MHM= zfZ^n2-M}aJfc)9X)nih|2^0mYf+}u3NactDQo1c@Ut^pU8?L>F}meE*TH5`FiuTs>Q_JTkoa;t zi(A?M0JNZp>voeC^4TsFfqxRWo3&eVCEBW}Az0up;a7pZ-!RS;oSz#>rvHVw^x?lH!jhp2@gEY@e2o0X5tPQ2~UQfrmC)FqNA;NO_z`-rL)6zx?k6ku*47jt!N z0!cj32d;RZ<0h+pbuGG)Yi%AugHF5J3xOLh>yS4{p~Q;9As$EFk%0!h!&~tc&a-jm zPqPZn=^SabRhazDd%K({F6EuJjm+4ov`+`?z9qDp`pVZfFiMckbt45tb|3%;-ge}C zs2Pss020a%x%eyK?t$U??sbiH{{U!O!55hfx`ayvI~fBrybY4OMn+^|&g`P?R$vba z(L5`xc!oE)f;5LuiNeo+cOZt^46=`!M#AUkQ?zA>JOf=%gtYGo9~5bNJ*3cU8lBXr zziiZ1q0`}PBjIAo$vb37fH{p=l^afTUZ<@1ZqLVF8PczO7_e$PWfpBdtsIGLEs`Ry zvL;Zm?r7s(_UNP>;I11C*8U^-d*Ppm8%s8KdZx5e+8s{yN0!eD2J+d<6Xwn%b;-aD zpS-vVpcF2C9C%wyvb(T0G0EZ`Mdh`!zt!i8$_XZPDDMmkbFm!!sv;+I$F!Bhtk*vc zd{uekXmtC{V_BcV(>=s(cP^$H{uqp5T^#+^idRU^+%^?PLK49F2gGj>{{U}z&sDg! z)qdZnPbxf8CAZnbs)MjH;1+1Zv|-q8AbK2ExY+zxPYu{x+4$);YK>sGv0T_~o*UTM z4ZBJlEQIb;(K4il3d=wp`&j<~VR&xp-@|%*3v&&r7LYedWJrQU9f<_+h`BAC0nZiA z&G8l=4|p#});wirCYyhJtQyKJGwxMp&eC!NVMY(!GAgkgnzi8x^$jn<6Kg&vzOiF3 zo2kFpB3~xj7an0%G3U(wCN1+4Sn}U^{Pp868fyC2jbQPGwVlPCx%o7SZf%j*WM>48 z%2k*9s~;*)9xj)1wTd07oE}dko0QRO*8dfegxK4}|1_bu=0RAcf zZT|p-?(Z)j!58ZV?!Tkr=SjJ_CM48sfCub=?pF091-9*xLm*Ps!rIwwB)VzB z;(s|5A#4Mb!6z9daDTWPKR(3z*Tx8Z1J8cZapa^v zbcQI|ogyrDg1`nKbOeQ5p8H3jucZDD{74=kuw4^azVkJBvTfc`hPj9V!EF4}E*m`w zTpyT*HS=39%bfh{)PGvk)4WfrX`0@ft3`2}iw%n#Fve9nBbCE4?fE%Ba|@gR00*{w zVfc4$E#~n)k!btb%W5!0psvyhGvhf2INgAAoyQnreDmS&j<0+peJ+b8oBge%>T;}i z7F$e?;DUrknc7*m06_!`%i%3)<+HbtXtyo?hbeOYTVpK~9Apr4jOQzkNezGxs_=8iO2P3SrK~51 zBE0adHaGWI62P#$?C%ekBdYnIy9Ba;z>)|7a@jb+!11DZKKk!Ry3{;LaiU*dMB>j> z)2{x`XMh0gjW7iac`Czk>y4+Ap4!4Mv)Efh_Gm4_7*Xr2N80ECCd_CMLzR?%ke+TDcTKJAJV1dp7RAQQ&y6TvmikB8Qp)rO;>Tqgej z3eR@Bwv&5ld^bYU0hpqUCczlYq;*}gt^sV7>y3_&X`)S~>QYY)rll*{`C68u@;v48 zff3$DZ=FU}Re|TAITe?8rg)3Oeh|`he+8s^obX&M4JNZL?<(ByC77u~LY;sG1cev? z;(#OYg}wc@hi2NMY>Py-`w_JJ1di)G8{;<+NK`Ydksbn}>wt5(9=i7O(@e9H?*8rz zTWImVY2*>2)H9$1aiYV6gYbK$!ys2|53A+^)y zXOC`}ZWcK((ImspB<_e{x)mW#K*Msm0QGMj_(x0lgW!EO&qbEX3y~6C=_QAhzD5u* z+qAOBvD1JC7)P9RuQE*!!JiGhai{27zmcw9TuT|A_URRE>?3ouD?(g`Vha+^aDbzQ zY!G{!{d-5%{1<294+OGJCxh+;g5E?6DYQjsl$EWG~&&FL* zwGn)pXNvTwTTQxxMu~3WRE6Y>;dc$o84Z%8E=E3eAOeSsG~H)P@pgr)-gwtai^8Vt zOQuAd0^8dpZzkyA^4d0JLQxzBC7Ldo0D5+)JU6Czb`4X+ntT@)i>Jj6&HdGjIoZBF zvARjRZe4=~ZP+SH9oQ#^)Vy`@>%e*jqu~8E((3&guHm_~i3(xZkjW&iw2{uh>%4_Q z$-vKt>sPvFt*EulhL=8{r`uW;OBd9{zELEB<+H9z?{MZY2|1AELCS!3TA#(eO5Xnf z#3I`EOjYr2e`i`{}uoc|1 zvrJ_L0<$rZNdZF+zMOh;tRraR^3m0Te)O}34BL+?FgA_5Mh z44mL&C3DL7E8%vr;!Sz(~{{Xz$=%qknx&HuI2|44BTIXi)*NCBM z{?VZ7k(eSAO}@@W^5Z+O?UBO#kICcZK8B#N*0oItO4?SdePU+Em@O_IFEEVbJ4ntl ze=|TAFZCUI$4a^I_|+$$L>EmE5KPZ$GsM|GTOJDK!{B2(k%-RRLvwp~tHq~j?7xyvywa0@vs%uYUvpb8hV8Rd?Dv%wD7VvM8-upNl0yZz?dxsP@v_peU) zQ8d?Z>pl!yE4ypWZ7s}~i4plM@w#2*vJ^BwYLhBPjvIfL~$LYuwBD@g>3n2 zmn0F4f(K*a4+rSFt;yGJu4R+M7JH<=hAh3srT`%(RROj)%sj}|vJjFZ`L=TMFZOly z<%Pbd3s_w3w`;whyOskYNRG;(nSSaqQS&Yek(`PE&sgjl$;&;{-Lfv!3ac76k5O0K_Wdw|4xa(S%s~OOAndQ`@x^ZtUymp90Z)D3klMXhYF&SlMQV8XU z&+hPTYw&1M&Z}W?ge8@@iLDi#E~MLu-a{jYi0==9uw-pJR1m}h8;|Wdt#oY~_AO3$ znkSaV)#5Q7^PDM%A&a9Fjko|4ybKe(`sct84Z~xm>RvIPcXlZyqFcZOq(-|Pn5t#V zk^l*Sa5IM9wc{GLw}0X*=fBmjr@XjgM){w121(kZqOM88i~*hnPCQrQYk4MnO?y$d zGe4Oqx4C#kkF+Uv2ta2cM&>-RQm;=dR}rD-pB1)O$41-_?nqDfhI638|?NZ4rgLw3Q~hIIXOFIMeGE zEVwMUI}2CK1I@R z20(0NL}M!C5O7K6+ccV?J!F?O7mKTxQNrb{PIFXfsVjBcUmHBW%$s>Tk z8y*Ak>v&tn!&I3bEjH;0ytS~6Ld}FbNJb%4@56Dl46slX?XRbGD-AEk-V-*mypn0U zlnEk988NoghesJ11Rbq{J79yzuZ$s)t)*Z*yt1-_@vDKkN|S;J&JPEuI3olDUcK;J z;#8Vt{2mn69!9gB)*A*?TtO>45~l-k4$st+w1L=?03)?F{{R;LBMlb&Q;yTa_6c$Q zm2YPXKukkosziVx*1j#% z^((CxP`A`{>&Y#zE!;eEq2wRrK^X)d2qYd2b=RL8{BL)n{{UwFx(M9v4|@awCv`a( z$oW^O!BtipM>)6_4i_hYSb#`8@NwdPG5Gta z>oyYDcz(`(Cv>EJv(HwM;lmBPc7E-SROMSLM+d!mUZH(|tIoGNoy*+Z1qf!B=8?up z8{4aqoMnk106G8zY|XS6_S0Qj{jSY3(lU&a&&D_r< zvblLAk&832IV6+GB=zD zD+>8gBw10=w$L{>t2S8oCaZ1T+JGm5_Tno`n^@)94bIUp3O6`Ykfin7oD#XgDx`u{ zaF204&D678Jd;f`nIw`%EY8H_l20U)-l0+y;EJ-Nv=G21R;j!yCwSUliJ#JG7TYx^~kMpHu~+^B3(_-fUe2zckT&gf9Ky z0Kl0l{4Vfsh95#-8Eco4UdL~7_Gi<+!exymW96!{?PDo#HYQ!5V4nCm^-lolap~)G z6n+){%e1xM7M*)M%Ov*}@=GMUNGc?aK`M4XBX>VACV)O(@szP#c;8;Sf>`E;T}I+a zq>Y0ssgZ~vjtC>YVuw<>X6C)S#l8#C{B5Clm&RH=pJ2bSy0{lMjSI$%A<1&;z%Z59 z5@Y~AV5~Q8bK_bhnq9uTWj3{IJ=U0#gL4cLYq3r?K}HA#jKx0wrQpF8Dmg<$#zSYM0DB_z}!e|x0Z89!QU36 zw$ztO({wARvUhV8jkI8e`|?Qe#){5`j^ZL_kwI`g$khp84;m0nCY$#5wV&)pUnlJb zGv)c1-LrR1*ykC}Ij-AZ@qP8Tf#C36z@F~uCTqLtr@1&1NFB^2R4X$h$gbXF5+z0i z5Td(J7+rW)*GlmJ0K^Xu$u*?!YbBy#QIhZ8|&KsI-Xg9?|6wZ!PLb(WKu#YA#zCQGRiaR zK<_{o^<4(~!uIp~LJLa*uqEs;iMJ%tNa_?o;hW4KK6AM6kCI3=n$10(Xl9}3+oRgmZ@bS)L?LuX|AO% zz$zFpg@9)S5`I9$1to%v65U6`I+ccvXRFB8vRYi*tl~(nV`8AYAXWh6IV!{`$-=*Q ziX0JT=4$C6kcuisc4URYjra+{PIPS z8p!Ns2&@)0-cE8CZ=C8^ns!JC8aY@38AM278_58vc=^MOfP1Vl>3$>d_x4|(Yjxqx zS_Vxz<^ddLF@i)QG63HgQc#7;h=_0*H$RhVS`F@-eEu2KA$>zx)AZL{FD8C9Z8e>< z7RO>%Xo+bf05%L?b+#cq*Teq+iIR9@P=igmv9TTkC{~bp|<3Z zpzU&sE7!GKTOS&DQfup7DdX2}FGM<>m94bWww-o2X5-3VnG!&oWE)#QDfylvAILR* zGV4jyZY-}6BfPeME@VkkK&(2A%1FrQMtC^~kZ1y?$8)+yV%QSI>XJ#8fk7z$O06gx) zMfjWH&k8)&F=|?Vo2XBTqF3@{gr+}w`N4q>2rrcwPnhl7ov())KCU&i@$8;0(`U1_ zfu+@fxSo5JhhqYg$AEG+ags(*GKT;I>36bdx`RvL$JabF28kM5SV5=#s84qzB#It3 z48vjCzFOmRhC9Os09Iz9;Xj6c4e>^SXEujCS9-Ls1+=R?(nq@rr6dMURs6u1AgK&P ze7T?sek5I1$Kn-@r-9#2(wD_n?RMI#Aedau4066OOCUlTGFek3xoo#0yz5K&58{i< zYkRBNt#tTaJ zWI0ebF9#&{Gy!(|;KrY28?e{4%ZU~tD#ZtdZkWDTa8Pr!F~}GsVDNB3$9>_QPf+lG z!@KJd6^zQU+o9LqFwKbD3Tjg@3}W!ed0f(g!5TRX#l0Q@O)r;EFf5KpYf zFPI{p#gaHk!CZ`Q%PVB(b}|4UfI(d!!8)du@TbLo6SuO5SiCx{k%)!Q*_lSjj8f$M zrG_GxX6$V^Ew~k+4=3^WiXX*3GM`msk~rpq-((Lu${-$6mOxbF8;cxpa!x|FK(pX2 zE?M;`uWWT$to0O>v@+Z5D`grO%-0ts*ae72>;UNO3X(B~T|YzAd@JIKTf|xf(@Uk= z{k`BCn%V7L$`(;F*ENX=jSs!%>4ou)5NtSq;^(SCZx=X@oZ= zNf0oEMDk>)W>Dc6kCX!Dqu}j3TGVZBn#xN}7wl0aT0+4jdd8mu%F&ld7jTgpRx&(} z9g4GK<^w)^u6RaWGsDBhR+o@mczW{wTU#8ndF+bv1!n~S84(lsjt6Zk{{RzR==b`* zlj7|MLnwX?Ji{46?1P zK5L{)%P>+vCmjYx;zEu{uR*o=RbX`eKF;sNNvCVqD=XWJjk4M5LwFFiyvjhs23)D# zC(WIJV?5V{ykdMMYpCir4V`M{AFUVaUuC#{H zVY9H0ZMC~R0x}|Bm1LXC1iXZN&cHDwoN<5?Q1I}IIj(fgM%PoeiDa3KTan|0P9aB; z7!8HXDF_QN3<&^pSrN(Q@WUhZt1P9!k}NkJJtzZ3t^U_?kf-*1w|at^yyJBpyfEa3 zz*FU{PNgM68|5K}#m2^F+e>+(%W7|+yp2q6d1Q|ALA9hzM{ybUWsR~v_m&26S=M(q z7Oy;Z>{Vukp_K}T4u_~9r{&>{# zRU=)iy}{ohb#39;09w479~9}<_Z}s5vA&6Xl-IX;%Xf%_#~kRHW_1NsrCrRMkz`Ux zBA&D1ON#@fT3hM*l3w0HmRpC4NG@C`U~Q&DwouB;gc%`-P+YLATay*F)b~1E^T%Z* zlS6LOOrmoU7{q~2Mp3{#ZEu(i;Nz@c9i`p;mhpIVN=+*FQ8G`Zv#MJZXwoDKSOElz zHnow(;p*=JSGS9( z-9vP3uVOLkQ?4DH+qnczgh&fWpteH0MyfgGuU6VvqlZe4=1XW`N$#3kbPW`(u_XD^ zB*njO(%8oQV(laYP9X9sf{{Tgqyz?=S%RwU+MOdZF2WIWNZwY ze9f53=(@;m>EFhdXfeJ8C-*x@i)V< zY1dkQxuV=X_l9q6)okvt2#W=Wkk9v+s2E-U05DOrpDg3uekXiA)4Xk=$#h5A^@}K( zt?rU0a>_DHY^n<~o^Svx3Ce?#t>AtV)NO5&{{X@MA-374>N1F#?lP^WT@0enwV6Wg zCk`?RRlzs~NCe;Tvp^Yd19&Ho`h@=g_I1^G6XHeTztJHA-X0_B>#a6}x8f7)3-Uzz4?@_z0& zW9IwW$8PmU#k%u&U&FSVj-dzIt>e13O+w-cTHHW|z%vB}nUp@!<$`h>1d*Ki&9}!Z zohI8#lV9-s^4MO)?-r$PFk>iEK^qDpnH7LX!BL#DjjW*ee-Jbu5qPt~V@}Yn?$D zh(1Z<{{V%WpMtOKAc9G?OKZ@v!*Lz1+|EfJPSDJE5rRLxlZ@>(zwrCSz9RT#ZKy}8 zcvD-F%ItZO`H@KqMcr=o?}K(X+vL%Z<4i>!El?(^+m+4`%;3$AvE+8CrPW?H#X6^+i`1xvngHbrIksMOJ@v#54>u(k3217;wv_r#5$Ca_WyYlQUVu~^5Lna9kZ5+Piy@BN9 zmE#MZOH24`<2$WKPFr0)?JlK3bsNtw*`&7%DNW=yA1YN?>%{EXcuRD9BW?RphB6{{V_2x7R!!aW<1Sk8mw5#nr=|B_Y+s8&Ja;jZ`m`iqQs{ z2vVH}6o~*VA&Twpyfvs=X%|f$nl_!L++T~3@~}mS+1#|M(M2H$?F&Brin{}EJ9%w0 zSMe8vbsbv!!up2YBF<4I&7Hv!d2GfY==QQm(haW55pYT4hw5uKZ36GY4X0htXQS(} zYBp2J98sBWB8^pXBX5kVmp~soUAADiA1hZub!7*NJSVL9y7JS^(qWie@0Q@cp=C2k z7|O}CeXRnz045=T`Bj&0_5f4jJtJ1s^!+o!UKpBtZANB&OF*~?_VP+pjqo;0Lcv-= z<}SbjA;1bnK3@cAaCn7HowS#xXckNF4VhXi{Wwbqn>Os_MZ?B7e54c5GLm-&Y99>r zyK5bHTD$RWjQV7IMWy|Px0wOCa!Vws2Ja!1q8ywe~2aByWdxsj4v(*fy z;Uz5{f60G(bf&roeUHyEATe6L0^LIL(*D&y8U<{d&SUJeZax+zB@PyWyRM*$K zt(Ck&A#q^rS)s!rK+%1)>%`LsWmon@eBuKlN ziRb_xy-5IMkSbpf>i+=Py2ZWR_OR+Ui71kL3!%Ch@T|FIIU9g?+DJP|&JQ5i8kOdQ zV{@X~-cJCzKk+hDNg-n3GWiY(1gZvLM%BwRgS;v=VU=|x#rJ*rpbHk-r1~bbnx2U? z_ENzqdv-1zkaK{?05=>Tlm_P`2NlxlddG%5OQl1%VG{TUo4H3fx?np558i)^tGYAwYrl$%b62FCjn)f8rJKjFL;@`@4qJ;)q3Q zsK*inPnZysMu289mphPe8;(SuDz0s&r4_D*Y$CTyxNad~I$uT@Th3hKTXs?z4#kKb zaKsS7fLvxW1w*ufo^X0l2HZXd@deJIronG&XSda55!zqE8cS(1SZ8vAR{t2=MPmFqAw{K~t>umO#!#p;zu|p!>NWd=f zu|SNL&QyQ~MhWKt9;KT3`{KufqM8f6Uq;j9x4F_-Lwjcy8BM7S@!Y0x+WuQu?NTNSMh(5Ere^~&CM6u(7Q{+9kXD8@)gEYZsC%0NZNgt-sT&aBDjKNnh8`q Taz??G)QpA%a6ugAfIt7)+#d(c literal 0 HcmV?d00001