From 5546ca784520d9c2bbdf373bd1608e53b9cd530e Mon Sep 17 00:00:00 2001 From: patel-zeel Date: Thu, 26 Sep 2024 08:26:10 +0530 Subject: [PATCH] major changes. Checkpoint --- .gitignore | 3 + README.md | 11 +- garuda/annotate.py | 247 +++++++++ garuda/box.py | 684 ++++++++++++++++++++++++ garuda/bulk_ops.py | 5 +- garuda/{ops.py => core.py} | 669 ++++++++++++----------- garuda/od.py | 4 +- garuda/plot.py | 2 +- garuda/utils.py | 27 + lab/test.tif | Bin 64902 -> 0 bytes lab/testing.ipynb | 773 ++------------------------- requirements.txt | 2 +- tests/{test_ops.py => test_utils.py} | 2 +- 13 files changed, 1375 insertions(+), 1054 deletions(-) create mode 100644 garuda/annotate.py create mode 100644 garuda/box.py rename garuda/{ops.py => core.py} (53%) create mode 100644 garuda/utils.py delete mode 100644 lab/test.tif rename tests/{test_ops.py => test_utils.py} (98%) diff --git a/.gitignore b/.gitignore index 08141be..dc00b5f 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,6 @@ __pycache__/ *.pyc *.egg-info/ garuda/_version.py +lab/* +!lab/*.ipynb +build \ No newline at end of file diff --git a/README.md b/README.md index c3df1d0..4187879 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,11 @@ Latest version: pip install git+https://github.com/patel-zeel/garuda ``` +## Non-circular imports +* core -> box +* box -> utils +* core + box -> annotate + ## Terminology | Term | Description | @@ -157,4 +162,8 @@ img = plt.imread('data/images/22.32,87.93.png') fig, ax = plt.subplots() ax = plot_webm_pixel_to_geo(img, img_center_lat, img_center_lon, zoom, ax) -``` \ No newline at end of file +``` + +### Why `OBBLabel` and not `OBBLabels` +* Non-vectorized operations are easier to understand and debug. +* It'd be easy for us to separate out exact false positives and false negatives at a single object level. \ No newline at end of file diff --git a/garuda/annotate.py b/garuda/annotate.py new file mode 100644 index 0000000..b6a2e48 --- /dev/null +++ b/garuda/annotate.py @@ -0,0 +1,247 @@ +import os +import cv2 +from glob import glob +import numpy as np +from ipyleaflet import GeomanDrawControl +import leafmap +from copy import deepcopy +import geojson +from ipywidgets import Button, Label, HBox, Dropdown +from IPython.display import display +from garuda.core import geo_to_webm_pixel, webm_pixel_to_geo, xywhr2xyxyxyxy +from garuda.box import OBBLabel +from shapely.geometry import Polygon + +class AnnotationTool: + def __init__(self, labels, classes, zoom, cache_dir, clear_cache=False): + self.original_labels = deepcopy(labels) + self.labels = deepcopy(labels) + self.classes = classes + self.zoom = zoom + self.cache_dir = cache_dir + os.makedirs(self.cache_dir, exist_ok=True) + if clear_cache: + label_files = glob(os.path.join(self.cache_dir, "label_*.geojson")) + for label_file in label_files: + os.remove(label_file) + + # initialize + self.index = 0 + + # Map + self.m = leafmap.Map(center=(27, 77), zoom=self.zoom) + self.m.add_basemap("Esri.WorldImagery") + self.m.remove_control(self.m.draw_control) + self.draw_control = GeomanDrawControl(position='topright') + + def on_draw(*args, **kwargs): + self.status_label.value = "Submit the label to update it." + self.disable_buttons() + + self.draw_control.on_draw(on_draw) + + self.draw_control.rectangle = { + "pathOptions": { + "fillColor": "#fca45d", + "color": "#fca45d", + "fillOpacity": 0.0 + } + } + + self.m.add_control(self.draw_control) + + # Interface elements + # current label + self.show_label = Label(f"Label {self.index+1}/{len(labels)}") + + # status label + self.status_label = Label("") + + # next_button + self.next_button = Button(description="next") + self.next_button.on_click(self.next_button_clicked) + + # previous_button + self.previous_button = Button(description="previous") + self.previous_button.on_click(self.previous_button_clicked) + + # submit button + self.submit_button = Button(description="submit") + self.submit_button.on_click(self.submit_button_clicked) + + # reset button + self.reset_button = Button(description="reset_current_label") + self.reset_button.on_click(self.reset_button_clicked) + + # classes dropdown + self.classes_dropdown = Dropdown(options=self.classes) + self.classes_dropdown.on_trait_change(self.on_dropdown_change, 'value') + + display(self.show_label) + display(self.status_label) + display(HBox([self.submit_button, self.previous_button, self.next_button, self.reset_button, self.classes_dropdown])) + display(self.m) + + # initialize + loaded_from_cache = self.show_current_label() + while loaded_from_cache: + self.next_button_clicked() + loaded_from_cache = self.show_current_label() + + def show_current_label(self): + self.disable_buttons() + loaded_from_cache = False + if os.path.exists(f"{self.cache_dir}/label_{self.index}.geojson"): + with open(f"{self.cache_dir}/label_{self.index}.geojson", "r") as f: + data = f.read().strip() + if data == "Empty_label": + self.labels[self.index] = None + else: + feature = geojson.loads(data)['features'][0] + self.labels[self.index] = OBBLabel.from_geojson(feature) + self.enable_buttons() # allow to move around if label is already present + loaded_from_cache = True + + label = self.labels[self.index] + if label is None: + self.draw_control.data = [] + self.status_label.value = "No label available." + original_label = self.original_labels[self.index] + self.m.set_center(original_label.properties['center_lon'], original_label.properties['center_lat'], zoom=self.zoom) + self.enable_buttons() # allow to move around if label is empty + else: + if loaded_from_cache: + self.status_label.value = "Label loaded from cache. Submit only to make changes." + else: + self.status_label.value = "Label is a valid polygon. Make changes if needed and submit." + # set pov + self.m.set_center(label.properties['center_lon'], label.properties['center_lat'], zoom=self.zoom) + + # show current label + feature = label.to_geojson(source=None, task_name=None) + # show the boundary in red color + feature['properties']['style'] = {'color': self.get_color(feature), 'fillColor': self.get_color(feature), 'fillOpacity': 0.0} + self.draw_control.data = [] # first clear the existing data to trigger the changes in GUI + self.draw_control.data = [feature] + self.classes_dropdown.value = feature['properties']['class_name'] + + # update label + self.show_label.value = f"Label {self.index+1}/{len(self.labels)}" + return loaded_from_cache + + def disable_buttons(self): + self.next_button.disabled = True + # self.previous_button.disabled = True + + def enable_buttons(self): + self.next_button.disabled = False + # self.previous_button.disabled = False + + @staticmethod + def get_color(feature): + class_name = feature['properties']['class_name'] + if class_name == "CFCBK": + # red + return "#ff0000" + elif class_name == "FCBK": + # orange + return "#ffa500" + elif class_name == "Zigzag": + # green + return "#00ff00" + else: + # blue + return "#0000ff" + + def submit_button_clicked(self, *args, **kwargs): + if len(self.draw_control.data) == 0: + self.labels[self.index] = None + # remove label from cache + cache_path = f"{self.cache_dir}/label_{self.index}.geojson" + with open(cache_path, "w") as f: + f.write("Empty_label") + self.enable_buttons() # allow to move around if label is empty + return + + feature = self.draw_control.data[-1] + try: + assert feature['geometry']['type'] == 'Polygon' + except AssertionError: + if feature['geometry']['type'] != 'Polygon': + self.status_label.value = "Invalid label. Please correct it or delete it." + return + + coords = [] + for lon, lat in feature['geometry']['coordinates'][0]: + x, y = geo_to_webm_pixel(lat, lon, self.zoom) + coords.append([x, y]) + coords = np.array(coords, dtype=np.float32) + + (x, y), (w, h), r = cv2.minAreaRect(coords) + r = np.deg2rad(r) + rect = xywhr2xyxyxyxy(np.array([x, y, w, h, r])) + + coords = [] + for pair in rect: + lat, lon = webm_pixel_to_geo(pair[0], pair[1], self.zoom) + coords.append([lon, lat]) + poly = Polygon(coords) + + feature['geometry']['coordinates'] = [list(poly.exterior.coords)] + feature['properties']['source'] = 'hand_validated' + feature['properties']['task_name'] = 'hand_validation' + feature['properties']['class_name'] = self.classes_dropdown.value + self.labels[self.index] = OBBLabel.from_geojson(feature) + self.cache_label() + self.show_current_label() + self.enable_buttons() + + def cache_label(self): + cache_path = f"{self.cache_dir}/label_{self.index}.geojson" + feature = self.labels[self.index].to_geojson(source=None, task_name=None) + collection = geojson.FeatureCollection([feature]) + with open(cache_path, "w") as f: + geojson.dump(collection, f) + + def reset_button_clicked(self, *args, **kwargs): + original_label = self.original_labels[self.index] + self.labels[self.index] = original_label + self.cache_label() + self.show_current_label() + + def on_dropdown_change(self, old, new): + if new != self.labels[self.index].properties['class_name']: + self.labels[self.index].properties['class_name'] = new + self.cache_label() + self.show_current_label() + + def next_button_clicked(self, *args, **kwargs): + # show next label + if self.index >= (len(self.labels) - 1): + pass # do nothing + else: + self.index += 1 + self.show_current_label() + + def previous_button_clicked(self, *args, **kwargs): + # show next label + if self.index <= 0: + pass # do nothing + else: + self.index -= 1 + self.show_current_label() + + def to_geojson(self): + features = [] + labels = glob(f"{self.cache_dir}/label_*.geojson") + for label in labels: + with open(label, "r") as f: + feature = geojson.load(f)['features'][0] + features.append(feature) + collection = geojson.FeatureCollection(features) + return collection + + def save_to_geojson(self, save_path): + collection = self.to_geojson() + with open(save_path, "w") as f: + geojson.dump(collection, f) \ No newline at end of file diff --git a/garuda/box.py b/garuda/box.py new file mode 100644 index 0000000..275d372 --- /dev/null +++ b/garuda/box.py @@ -0,0 +1,684 @@ +import numpy as np +from numpy import ndarray + +from beartype import beartype +from beartype.typing import Sequence, Tuple +from jaxtyping import Float, Int, jaxtyped +import warnings +from einops import rearrange +from geopy import distance +from copy import deepcopy +from geemap import geemap +from hashlib import sha256 + +from garuda.core import xywh2xyxy, local_to_geo, geo_to_local + + +# This class enforces type checking using beartype and jaxtyping on all the methods. +class AutoTypeChecker: + def __init_subclass__(cls, **kwargs): + super().__init_subclass__(**kwargs) + for attr, value in cls.__dict__.items(): + if callable(value): + setattr(cls, attr, jaxtyped(typechecker=beartype)(value)) + + +class BB(AutoTypeChecker): + def __init__( + self, + geo_box: Float[ndarray, "4 2"] | Float[ndarray, "2 2"], + class_name: str, + confidence: float | None, + ): + """ + Base class for bounding boxes. + + Parameters + ---------- + geo_box: ndarray + Bounding box in geo coordinates. + For OBB: [[lon1, lat1], [lon2, lat2], [lon3, lat3], [lon4, lat4]] + For AABB: [[lon1, lat1], [lon2, lat2]] + + class_name: str + Class name of the bounding box. + + confidence: float + Confidence of the bounding box. This will be ignored if the child class represents a label. + """ + + if type(self) is BB: + raise TypeError( + f"Instantiation of '{self.__class__.__name__}' is not allowed. Use one of the children classes instead." + ) + self.geo_box = geo_box + self.class_name = class_name + self.confidence = confidence + self.properties = { + "geo_box": self.geo_box.tolist(), + "class_name": self.class_name, + "confidence": self.confidence, + } + + self.add_additional_properties() + + def add_additional_properties(self): + self.properties["max_lon"] = self.geo_box[:, 0].max().item() + self.properties["min_lon"] = self.geo_box[:, 0].min().item() + self.properties["max_lat"] = self.geo_box[:, 1].max().item() + self.properties["min_lat"] = self.geo_box[:, 1].min().item() + self.properties["center_lat"] = ( + self.properties["max_lat"] + self.properties["min_lat"] + ) / 2 + self.properties["center_lon"] = ( + self.properties["max_lon"] + self.properties["min_lon"] + ) / 2 + self.properties["width_of_box"] = distance.geodesic( + (self.properties["max_lat"], self.properties["min_lon"]), + (self.properties["max_lat"], self.properties["max_lon"]), + ).meters + self.properties["height_of_box"] = distance.geodesic( + (self.properties["min_lat"], self.properties["max_lon"]), + (self.properties["max_lat"], self.properties["max_lon"]), + ).meters + + if self.geo_box.shape == (4, 2): # OBB + side1_length = distance.geodesic( + (self.geo_box[0][1], self.geo_box[0][0]), + (self.geo_box[1][1], self.geo_box[1][0]), + ).meters + side2_length = distance.geodesic( + (self.geo_box[1][1], self.geo_box[1][0]), + (self.geo_box[2][1], self.geo_box[2][0]), + ).meters + self.properties["length_of_object"] = max(side1_length, side2_length) + self.properties["width_of_object"] = min(side1_length, side2_length) + elif self.geo_box.shape == (2, 2): # AABB + self.properties["length_of_object"] = max( + self.properties["width_of_box"], self.properties["height_of_box"] + ) + self.properties["width_of_object"] = min( + self.properties["width_of_box"], self.properties["height_of_box"] + ) + else: + raise ValueError( + f"Invalid {self.geo_box=}. Box should be either OBB with (4, 2) shape or AABB with (2, 2) shape." + ) + + self.properties["area"] = ( + self.properties["length_of_object"] * self.properties["width_of_object"] + ) + + # unique ID + self.properties["id"] = sha256(f"{self.properties}".encode()).hexdigest() + + @staticmethod + def _from_ultralytics( + label: ( + Float[ndarray, "10"] + | Float[ndarray, "9"] + | Float[ndarray, "6"] + | Float[ndarray, "5"] + ), + classes: Sequence, + image_width: int, + image_height: int, + ) -> Tuple[Float[ndarray, "8"] | Int[ndarray, "4"], str, float | None]: + """ + Common operations to convert Ultralytics format to OBB or AABB. + + Parameters + ---------- + bb_type: Literal["obb", "aabb"] + Type of bounding box. "obb" for oriented bounding box and "aabb" for axis-aligned bounding box. + + label: ndarray + Label in Ultralytics format. + For OBB: [class_id, x1, y1, x2, y2, x3, y3, x4] or [class_id, x1, y1, x2, y2, x3, y3, x4, confidence] + For AABB: [class_id, x, y, w, h] or [class_id, x, y, w, h, confidence] + + classes: sequence + Sequence of class names. We will use this to get the class name from the class_id. + + image_width, image_height: int + Original image size. + + Returns + ------- + box: ndarray + OBB: [x1, y1, x2, y2, x3, y3, x4, y4] + AABB: [x1, y1, x2, y2] + + class_name: str + Class name + + confidence: float + Confidence of the bounding box. 1.0 if confidence is not provided. + """ + class_name = classes[int(label[0])] + if len(label) in (5, 6): # AABB + box = label[1:5] + box = xywh2xyxy(box[None, ...]).ravel() + confidence = label[5] if len(label) == 6 else 1.0 + elif len(label) in (9, 10): # OBB + box = label[1:9] + confidence = label[9] if len(label) == 10 else 1.0 + + x = box[::2] + y = box[1::2] + + # scale to original image size + # x = x * image_width + # y = y * image_height + + box = np.stack((x, y)).T.ravel()#.round().astype(int) + return box, class_name, confidence + + @classmethod + def from_ultralytics_gms( + cls, + label: ( + Float[ndarray, "10"] + | Float[ndarray, "9"] + | Float[ndarray, "6"] + | Float[ndarray, "5"] + ), + classes: Sequence, + zoom: int, + image_center_lat: float, + image_center_lon: float, + image_width: int, + image_height: int, + ) -> "BBLabel": + """ + Load label from Ultralytics format. + + Parameters + ---------- + label: ndarray + Label in Ultralytics format. + For OBB: [class_id, x1, y1, x2, y2, x3, y3, x4] or [class_id, x1, y1, x2, y2, x3, y3, x4, confidence] + For AABB: [class_id, x, y, w, h] or [class_id, x, y, w, h, confidence] + + classes: sequence + Sequence of class names. We will use this to get the class name from the class_id. + + image_width, image_height: int + Original image size. + + zoom: int + Zoom level of the image. + + image_center_lat, image_center_lon: float + Latitude and longitude of the center of the image + """ + box, class_name, confidence = cls._from_ultralytics( + label, classes, image_width, image_height + ) + geo_box = cls._local_box_to_geo( + box, zoom, image_center_lat, image_center_lon, image_width, image_height + ) + + instance = cls(geo_box, class_name, confidence) + instance.properties["image_width"] = image_width + instance.properties["image_height"] = image_height + instance.properties["image_center_lat"] = image_center_lat + instance.properties["image_center_lon"] = image_center_lon + instance.properties["zoom"] = zoom + return instance + + @staticmethod + def _local_box_to_geo( + box: Float[ndarray, "8"] | Float[ndarray, "4"], + zoom: int, + image_center_lat: float, + image_center_lon: float, + image_width: int, + image_height: int, + ) -> Float[ndarray, "4 2"] | Float[ndarray, "2 2"]: + x = box[::2].tolist() + y = box[1::2].tolist() + if len(box) == 8: # OBB + lat1, lon1 = local_to_geo( + x[0], + y[0], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + lat2, lon2 = local_to_geo( + x[1], + y[1], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + lat3, lon3 = local_to_geo( + x[2], + y[2], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + lat4, lon4 = local_to_geo( + x[3], + y[3], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + return np.array([[lon1, lat1], [lon2, lat2], [lon3, lat3], [lon4, lat4]]) + elif len(box) == 4: # AABB + lat1, lon1 = local_to_geo( + x[0], + y[0], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + lat2, lon2 = local_to_geo( + x[1], + y[1], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + return np.array([[lon1, lat1], [lon2, lat2]]) + else: + raise ValueError( + f"Invalid {box=}. Box should be either OBB with 8 elements or AABB with 4 elements." + ) + + @staticmethod + def _geo_box_to_local( + geo_box: Float[ndarray, "4 2"] | Float[ndarray, "2 2"], + zoom: int, + image_center_lat: float, + image_center_lon: float, + image_width: int, + image_height: int, + ) -> Float[ndarray, "8"] | Float[ndarray, "4"]: + lons = geo_box[:, 0] + lats = geo_box[:, 1] + if geo_box.shape == (4, 2): + x1, y1 = geo_to_local( + lats[0], + lons[0], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + x2, y2 = geo_to_local( + lats[1], + lons[1], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + x3, y3 = geo_to_local( + lats[2], + lons[2], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + x4, y4 = geo_to_local( + lats[3], + lons[3], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + return np.array([x1, y1, x2, y2, x3, y3, x4, y4]) + elif geo_box.shape == (2, 2): + x1, y1 = geo_to_local( + lons[0], + lats[0], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + x2, y2 = geo_to_local( + lons[1], + lats[1], + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + return np.array([x1, y1, x2, y2]) + + def to_geojson(self, source, task_name): + if source is None: + if "source" in self.properties: + source = self.properties["source"] + else: + source = "Drawn|Azure Maps Satellite" + if task_name is None: + if "task_name" in self.properties: + task_name = self.properties["task_name"] + else: + task_name = "" + + if self.geo_box.shape == (4, 2): # OBB + box = self.geo_box + elif self.geo_box.shape == (2, 2): # AABB + # convert AABB to OBB format + lon1, lat1 = self.geo_box[0] + lon2, lat2 = self.geo_box[1] + box = np.array([[lon1, lat1], [lon2, lat1], [lon2, lat2], [lon1, lat2]]) + else: + raise ValueError( + f"Invalid {self.geo_box=}. Box should be either OBB with (4, 2) shape or AABB with (2, 2) shape." + ) + + points = rearrange(box, "n d -> 1 n d").tolist() + # append first point to close the polygon + points[0].append(points[0][0]) + + feature = { + "type": "Feature", + "properties": self.properties, + "geometry": { + "type": "Polygon", + "coordinates": points, + }, + } + return feature + + def to_ultralytics_obb( + self, + classes: Sequence, + zoom: int | None, + image_center_lat: float | None, + image_center_lon: float | None, + image_width: int | None, + image_height: int | None, + ) -> Float[ndarray, "9"] | Float[ndarray, "10"]: + """ + To Ultralytics oriented bounding box (OBB) format. + + Parameters + ---------- + classes: sequence + Sequence of class names. We will use this to get the class_id from the class name. + + zoom: int + Zoom level of the image. + + image_center_lat, image_center_lon: float + Latitude and longitude of the center of the image + + image_width, image_height: int + Width and height of the image. + + Returns + ------- + label: ndarray + Label/Detection in Ultralytics format. + Label: [class_id, x1, y1, x2, y2, x3, y3, x4] + Detection: [class_id, x1, y1, x2, y2, x3, y3, x4, confidence] + """ + assert isinstance( + self, (OBBLabel, OBBDetection) + ), f"this method should not be called on {self.__class__.__name__}. It is only for 'OBBLabel' and 'OBBDetection' instances." + + def auto_set(key): + nonlocal zoom, image_center_lat, image_center_lon, image_width, image_height + if eval(key) is None: + if key in self.properties: + return self.properties[key] + else: + raise ValueError( + f"'{key}' is neither provided nor found in the properties." + ) + else: + return key + + zoom = auto_set("zoom") + image_width = auto_set("image_width") + image_height = auto_set("image_height") + image_center_lat = auto_set("image_center_lat") + image_center_lon = auto_set("image_center_lon") + + class_id = classes.index(self.class_name) + + box = self._geo_box_to_local( + self.geo_box, + zoom, + image_center_lat, + image_center_lon, + image_width, + image_height, + ) + + label = np.zeros(10) * np.nan + label[0] = class_id + label[1:9] = box + if isinstance(self, OBBLabel): + # OBBLabel doesn't have confidence attribute + return label[:9] + elif isinstance(self, OBBDetection): + label[9] = self.confidence + return label + else: + raise ValueError( + f"Invalid instance type: {self.__class__.__name__}. This error should not have occurred because it should have been caught by the assert statement at the beginning of the method." + ) + + def to_ultralytics_aabb(self, classes: Sequence) -> Float[ndarray, "6"]: + """ + To Ultralytics axis-aligned bounding box format. + """ + class_id = classes.index(self.class_name) + x = self.geo_box[::2] / self.original_width + y = self.geo_box[1::2] / self.original_height + x_min, x_max = np.min(x), np.max(x) + y_min, y_max = np.min(y), np.max(y) + box = np.array([x_min, y_min, x_max, y_max]) + + label = np.zeros(6) * np.nan + label[0] = class_id + label[1:5] = box + label[5] = self.confidence + raise NotImplementedError("This method is not implemented yet.") + + def visualize(self, zoom): + geojson_box = self.to_geojson(None, None) + Map = geemap.Map( + center=[self.properties["center_lat"], self.properties["center_lon"]], + zoom=zoom, + ) + Map.add_basemap("SATELLITE") + Map.add_geojson(geojson_box, layer_name=self.class_name) + return Map + + +class BBLabel(BB): + def __init__( + self, + geo_box: Float[ndarray, "4 2"] | Float[ndarray, "2 2"], + class_name: str, + confidence: float | None = None, + ): + super().__init__(geo_box, class_name, confidence=confidence) + del self.confidence # Remove confidence from the label + + # We didn't implement because a path can have multiple labels/predictions. + # @classmethod + # def from_ultralytics_gms_path(cls, path: str, classes: Sequence, image_width: int, image_height: int, zoom: int): + # assert path.endswith(".txt"), "Only .txt files are supported." + # path = path.replace("%2C", ",") + # base_name = basename(path) + # base_name = base_name.replace(".txt", "") + # lat_str, lon_str = base_name.split(",") + # lat, lon = float(lat_str), float(lon_str) + # label = + # return cls.from_ultralytics_gms() + + @staticmethod + def _from_geojson(feature: dict): + class_name = feature["properties"]["class_name"] + points = feature["geometry"]["coordinates"] + points = deepcopy(points) # avoid modifying the original points + points = points[0] # ignore additional axis + # delete extra point which was added to close the polygon + points.pop() + + lon1, lat1 = points[0] + lon2, lat2 = points[1] + lon3, lat3 = points[2] + lon4, lat4 = points[3] + + geo_box = np.array([[lon1, lat1], [lon2, lat2], [lon3, lat3], [lon4, lat4]]) + + return geo_box, class_name + + def __repr__(self): + return f"""BBLabel(box={self.geo_box}, +class_name={self.class_name}, +length_of_object={self.properties['length_of_object']:.2f} m, +width_of_object={self.properties['width_of_object']:.2f} m, +""" + + +class BBDetection(BB): + pass + + +class AABBLabel(BBLabel): + @classmethod + def from_geojson(cls, label: dict) -> "AABBLabel": + geo_box, class_name = cls._from_geojson(label) + min_lon = geo_box[:, 0].min() + max_lon = geo_box[:, 0].max() + min_lat = geo_box[:, 1].min() + max_lat = geo_box[:, 1].max() + geo_box = np.array([[min_lon, min_lat], [max_lon, max_lat]]) + instance = cls(geo_box, class_name) + instance.properties = label["properties"] + return instance + +class AABBDetection(BBDetection): + pass + + +class OBBDetection(BBDetection): + # TODO: nms supress method + @classmethod + def from_ultralytics_gms( + cls, + label: Float[ndarray, "9"] | Float[ndarray, "10"], + classes: Sequence, + image_width: int, + image_height: int, + zoom: int, + image_center_lat: float, + image_center_lon: float, + ) -> "OBBDetection": + box, class_name, confidence = cls._from_ultralytics( + label, classes, image_width, image_height + ) + if confidence is None: + warnings.warn("Confidence is not provided. Setting it to 1.0.") + confidence = 1.0 + + geo_box = cls._local_box_to_geo( + box, zoom, image_center_lat, image_center_lon, image_width, image_height + ) + + return cls(geo_box, class_name, confidence, image_width, image_height) + + +class OBBLabel(BBLabel): + @classmethod + def from_geojson(cls, label: dict) -> "OBBLabel": + geo_box, class_name = cls._from_geojson(label) + instance = cls(geo_box, class_name) + + if "source" in label["properties"]: + instance.properties['source'] = label["properties"]["source"] + if "task_name" in label["properties"]: + instance.properties['task_name'] = label["properties"]["task_name"] + return instance + + @classmethod + def from_label_studio_csv( + cls, label: dict, zoom: int, image_center_lat: float, image_center_lon: float + ) -> "OBBLabel": + """ + Label Studio CSV label keys: {"x", "y", "width", "height", "rotation", "rectanglelabels", "original_width", "original_height"} + + x, y: top left corner of the bounding box. Scale: [0, 100] + width, height: width and height of the bounding box. Scale: [0, 100] + rotation: rotation of the bounding box in degrees. Scale: [0, 90] + rectanglelabels: class label + original_width, original_height: original image size + """ + x = label["x"] + y = label["y"] + width = label["width"] + height = label["height"] + rotation = label["rotation"] + original_height = label["original_height"] + original_width = label["original_width"] + + rotation_rad = np.radians(rotation) + cos_rot = np.cos(rotation_rad) + sin_rot = np.sin(rotation_rad) + + x_1 = x + width * cos_rot - height * sin_rot + y_1 = y + width * sin_rot + height * cos_rot + x_2 = x + width * cos_rot + y_2 = y + width * sin_rot + x_3 = x + y_3 = y + x_4 = x - height * sin_rot + y_4 = y + height * cos_rot + + # scale to original image size + (x_1, x_2, x_3, x_4) = map( + lambda x: round(x * label["original_width"] / 100), (x_1, x_2, x_3, x_4) + ) + (y_1, y_2, y_3, y_4) = map( + lambda y: round(y * label["original_height"] / 100), (y_1, y_2, y_3, y_4) + ) + box = np.array([x_1, y_1, x_2, y_2, x_3, y_3, x_4, y_4]) + + geo_box = cls._local_box_to_geo( + box, + zoom, + image_center_lat, + image_center_lon, + original_width, + original_height, + ) + class_name = label["rectanglelabels"][0] + return cls(geo_box, class_name) + + +class AABB: + pass diff --git a/garuda/bulk_ops.py b/garuda/bulk_ops.py index 2b3f0c2..97d490d 100644 --- a/garuda/bulk_ops.py +++ b/garuda/bulk_ops.py @@ -4,9 +4,12 @@ import pandas as pd from garuda.od import add_obb_to_label_studio_df from beartype import beartype +from beartype.typing import Sequence, Tuple from jaxtyping import jaxtyped -from garuda.ops import obb_to_aa +from garuda.utils import obb_to_aa +from garuda.core import BB + @jaxtyped(typechecker=beartype) def write_obb_labels_from_label_studio_csv(save_dir: str, df:pd.DataFrame, label_map: dict): diff --git a/garuda/ops.py b/garuda/core.py similarity index 53% rename from garuda/ops.py rename to garuda/core.py index 586163c..d9ef172 100644 --- a/garuda/ops.py +++ b/garuda/core.py @@ -1,278 +1,60 @@ -import numpy as np import cv2 -from numpy import ndarray import pyproj -from beartype.typing import Tuple, Union -import warnings -from beartype import beartype -from jaxtyping import Float, jaxtyped, Int +import geojson +import numpy as np +from numpy import ndarray import planetary_computer as pc from pystac_client import Client from pystac.extensions.eo import EOExtension as eo -from shapely.geometry import box +from shapely.geometry import box, Polygon import xarray as xr import rioxarray +from os.path import basename -# @jaxtyped(typechecker=beartype) -def local_to_geo(x: float, y: float, zoom: int, img_center_lat: float, img_center_lon: float, img_width: int, img_height: int) -> Union[Float[ndarray, "2"], Float[ndarray, "n 2"]]: - """ - Convert local coordinates to Web Mercator projection at a given zoom level for given image center and image dimensions. - - Parameters - ---------- - x : X-coordinate in local coordinates. Generally mentioned in the YOLO label. - Range: [0, 1] - Example: 0.5 - - y : Y-coordinate in local coordinates. Generally mentioned in the YOLO label. - Range: [0, 1] - Example: 0.5 - - zoom : Zoom level. - Range: [0, 20] - Example: 17 - - img_center_lat : Latitude of the center of the image. - Range: approx [-85, 85] (valid range for Web Mercator projection) - Example: 37.7749 - - img_center_lon : Longitude of the center of the image. - Range: [-180, 180] - Example: -122.4194 - - img_width : Width of the image in pixels. - Range: [0, inf] - Example: 640 - - img_height : Height of the image in pixels. - Range: [0, inf] - Example: 480 - - Returns - ------- - (x_webm, y_webm) : X and Y coordinates in Web Mercator projection. - Range: [0, 2^zoom * 256] - Example: (1000, 1000) - """ - - # Get image center in Web Mercator projection - image_center_webm_x, image_center_webm_y = geo_to_webm_pixel(img_center_lat, img_center_lon, zoom) - - # get delta_x_c: (0, 1) -> (-img_width/2, img_width/2) - # get delta_y_c: (0, 1) -> (-img_height/2, img_height/2) - delta_x = x * img_width - img_width/2 - delta_y = y * img_height - img_height/2 - - # Get bbox center in Web Mercator projection - bbox_center_webm_x = image_center_webm_x + delta_x - bbox_center_webm_y = image_center_webm_y + delta_y - - # Convert bbox center to geographic coordinates - bbox_geo = webm_pixel_to_geo(bbox_center_webm_x, bbox_center_webm_y, zoom) - bbox_geo = np.array(bbox_geo).T - - return bbox_geo +from beartype import beartype +from beartype.typing import Union, Sequence, Optional, Tuple, Literal +from jaxtyping import Float, Int, jaxtyped +import warnings +from einops import rearrange +from geopy import distance +from copy import deepcopy +from geemap import geemap +from hashlib import sha256 -@jaxtyped(typechecker=beartype) -def geo_to_webm_pixel(lat:Union[float, Float[ndarray, "n"]], lon:Union[float, Float[ndarray, "n"]], zoom:int) -> Tuple[Union[float, Float[ndarray, "n"]], Union[float, Float[ndarray, "n"]]]: - """ - Convert latitude and longitude to Web Mercator projection at a given zoom level. - - Parameters - ---------- - lat : Latitude in decimal degrees. - Range: approximately [-85, 85] - Example: 37.7749 - Beyond the specified range, the projection becomes distorted. - - lon : Longitude in decimal degrees. - Range: [-180, 180] - Example: -122.4194 - - zoom : Zoom level. - Range: [0, 20] - Example: 17 - - Returns - ------- - x : X-coordinate in Web Mercator projection. - Range: [0, 2^zoom * 128] - Example: 1000 +############################################################ +# Common utils +############################################################ - y : Y-coordinate in Web Mercator projection. - Range: [0, 2^zoom * 128] - Example: 1000 - """ - # Convert latitude and longitude to radians - lat_rad = np.radians(lat) - lon_rad = np.radians(lon) - - # Project latitude and longitude to Web Mercator - x = lon_rad + np.pi - y = np.pi - np.log(np.tan(np.pi/4 + lat_rad/2)) - - if np.any(y < 0): - warnings.warn(f"y-coordinate is negative. Latitude='{lat}' might be beyond the valid range of laitude for Web Mercator projection (approx [-85, 85]).") - elif np.any(y > 2*np.pi): - warnings.warn(f"y-coordinate is greater than 256*2^zoom. Latitude='{lat}' might be beyond the valid range of latitude for Web Mercator projection (approx [-85, 85]).") - - if np.any(x < 0): - warnings.warn(f"x-coordinate is negative. Longitude='{lon}' might be beyond the valid range of longitude for Web Mercator projection ([-180, 180]).") - elif np.any(x > 2*np.pi): - warnings.warn(f"x-coordinate is greater than 256*2^zoom. Longitude='{lon}' might be beyond the valid range of longitude for Web Mercator projection ([-180, 180]).") - - # Scale Web Mercator to zoom level - x = (128/np.pi)*(2**zoom) * x - y = (128/np.pi)*(2**zoom) * y - - return x, y - -def webm_pixel_to_geo(x:float, y:float, zoom:int) -> Tuple[Union[float, Float[ndarray, "n"]], Union[float, Float[ndarray, "n"]]]: +@jaxtyped(typechecker=beartype) +def get_latlon_from_gms_path(path: str) -> dict: """ - Convert Web Mercator projection to latitude and longitude at a given zoom level. + Get latitude and longitude from Google Maps Static Image Label path. Parameters ---------- - x : X-coordinate in Web Mercator projection. - Range: [0, 2^zoom * 256] - Example: 1000 - - y : Y-coordinate in Web Mercator projection. - Range: [0, 2^zoom * 256] - Example: 1000 - - zoom : Zoom level. - Range: [0, 20] - Example: 17 + path: Path of the Google Maps Static Image Label. + Example: "37.7749,-122.4194.txt" Returns ------- - lat : Latitude in decimal degrees. - Range: approximately [-85, 85] - Example: 37.7749 - - lon : Longitude in decimal degrees. - Range: [-180, 180] - Example: -122.4194 + dict: Dictionary containing latitude and longitude in string and float format. + Example: {"str": ("37.7749", "-122.4194"), "float": (37.7749, -122.4194)} """ - # Scale Web Mercator to radians - x_rad = x / (128/np.pi) / (2**zoom) - y_rad = y / (128/np.pi) / (2**zoom) - - if np.any(x_rad<0): - warnings.warn(f"x-coordinate is negative. x='{x}' might be beyond the valid range of x-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") - elif np.any(x_rad>2*np.pi): - warnings.warn(f"x-coordinate is greater than 2*pi. x='{x}' might be beyond the valid range of x-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") - - if np.any(y_rad<0): - warnings.warn(f"y-coordinate is negative. y='{y}' might be beyond the valid range of y-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") - elif np.any(y_rad>2*np.pi): - warnings.warn(f"y-coordinate is greater than 2*pi. y='{y}' might be beyond the valid range of y-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") - - # Inverse project Web Mercator to latitude and longitude - lon_rad = x_rad - np.pi - lat_rad = 2*np.arctan(np.exp(np.pi - y_rad)) - np.pi/2 - - # Convert latitude and longitude to degrees - lat = np.degrees(lat_rad) - lon = np.degrees(lon_rad) - return lat, lon - - -def obb_to_aa(yolo_label: Union[str, Float[ndarray, "n 9"], Float[ndarray, "n 10"]]) -> Union[str, Float[ndarray, "n 5"], Float[ndarray, "n 6"]]: - """ - Convert YOLO OBB labels with format [class_id, x1, y1, x2, y2, x3, y3, x4, y4] to YOLO axis-aligned format [class_id, x_c, y_c, width, height]. - - Parameters - ---------- - yolo_label: YOLO label (or path) in OBB format. - - Returns - ---------- - yolo_label: YOLO label in axis-aligned format. - """ - - if isinstance(yolo_label, str): - yolo_label = np.loadtxt(yolo_label, ndmin=2) - return obb_to_aa(yolo_label) # to trigger type/shape checking - - # Split the label into various components - class_id = yolo_label[:, 0:1] - confidence_scores = yolo_label[:, 9:] # will be (n, 0) array if confidence scores are not present - xyxyxyxy = yolo_label[:, 1:9] - - # Get the x and y coordinates - x = xyxyxyxy[:, ::2] - y = xyxyxyxy[:, 1::2] - - # Find the bounds - x_max = np.max(x, axis=1) - x_min = np.min(x, axis=1) - y_max = np.max(y, axis=1) - y_min = np.min(y, axis=1) - - # Convert to axis-aligned format - x_c = (x_max + x_min) / 2 - y_c = (y_max + y_min) / 2 - width = x_max - x_min - height = y_max - y_min - xywh = np.stack([x_c, y_c, width, height], axis=1) - - # Concatenate the class_id and confidence scores - yolo_label = np.concatenate([class_id, xywh, confidence_scores], axis=1) - - return yolo_label - - -def label_studio_csv_to_obb(x1, y1, width, height, rotation, label, label_map) -> Float[ndarray, "9"]: - """ - Convert from Label studio CSV (x1, y1, w, h, r) to YOLO OBB (x1, y1, x2, y2, x3, y3, x4, y4) format - - Parameters - ---------- - - x1: x-cordinate of one corner of the box, range(0, 100) - y1: y-cordinate of one corner of the box, range(0, 100) - w: width of the box, range(0, 100) - h: height of the box, range(0, 100) - r: rotation angle of the box in degrees - label: label of the box - label_map: dictionary mapping label to class_id - - Returns - ---------- - yolo_label: YOLO OBB label in the format [class_id, x1, y1, x2, y2, x3, y3, x4, y4] - """ - - rotation_rad = np.radians(rotation) - - cos_rot = np.cos(rotation_rad) - sin_rot = np.sin(rotation_rad) - - x_1 = x1 + width * cos_rot - height * sin_rot - y_1 = y1 + width * sin_rot + height * cos_rot - x_2 = x1 + width * cos_rot - y_2 = y1 + width * sin_rot - x_3 = x1 - y_3 = y1 - x_4 = x1 - height * sin_rot - y_4 = y1 + height * cos_rot - xyxyxyxy = np.array([x_1, y_1, x_2, y_2, x_3, y_3, x_4, y_4]) - - # Normalize to range(0, 1) - xyxyxyxy = xyxyxyxy / 100 - - label_id = np.array([label_map[label]]) - yolo_label = np.concatenate([label_id, xyxyxyxy]) - return yolo_label + assert path.endswith(".txt"), "Path should be a .txt file." + path = path.replace("%2C", ",") + base_name = basename(path) + base_name = base_name.replace(".txt", "") + lat_str, lon_str = base_name.split(",") + lat, lon = float(lat_str), float(lon_str) + return {"str": (lat_str, lon_str), "float": (lat, lon)} @jaxtyped(typechecker=beartype) -def get_sentinel2_visual(lat_c: float, lon_c: float, img_height: int, img_width: int, time_of_interest: str, max_cloud_cover: float, max_items: int = 10) -> xr.DataArray: +def get_sentinel2_visual(lat_c: float, lon_c: float, img_height: int, img_width: int, time_of_interest: str, max_cloud_cover: float, max_items: int = 10, nodata_window_size: int = 2) -> xr.DataArray: """ - Get Sentinel-2 image as a PNG file from Microsoft Planetary Computer API. + Get Sentinel-2 image as a xarray file from Microsoft Planetary Computer API. Image is centered at the given latitude and longitude with the given width and height. @@ -307,6 +89,11 @@ def get_sentinel2_visual(lat_c: float, lon_c: float, img_height: int, img_width: Range: [1, inf] Example: 10 + nodata_window_size: Size of the window to check for nodata values. + Range: [0, inf] + Example: 2 + We will check if a (nodata_window_size x nodata_window_size) window has all zeros. If yes, we will discard that item and move to the next one. + Returns ------- tif: GeoTIFF file of the Sentinel-2 image @@ -352,6 +139,21 @@ def get_least_cloudy_raster(sorted_items): try: assert cropped_raster.shape == (3, img_height, img_width) except AssertionError as e: + print(f"Shape of the raster is not as expected. Shape: {cropped_raster.shape}") + sorted_items.pop(0) + return get_least_cloudy_raster(sorted_items) + + try: + np_img = cropped_raster.values + np_img = rearrange(np_img, "c h w -> h w c") + # pad image to make it even sized + # padded_img = np.pad(np_img, ((0, 1), (0, 1), (0, 0)), mode='edge') + mask = (np_img == 0).all(axis=(2)) + n = nodata_window_size + result = (mask.reshape(np_img.shape[0]//n, n, np_img.shape[1]//n, n).all(axis=(1, 3))) + assert not result.any() + except AssertionError as e: + print("Nodata values found in the image. Discarding this item and moving to the next one.") sorted_items.pop(0) return get_least_cloudy_raster(sorted_items) @@ -369,7 +171,130 @@ def get_least_cloudy_raster(sorted_items): return get_least_cloudy_raster(sorted_items) @jaxtyped(typechecker=beartype) -def xyxyxyxy2xywhr(xyxyxyxy: Float[ndarray, "n 8"]) -> Float[ndarray, "n 5"]: +def local_to_geo(x: float, y: float, zoom: int, img_center_lat: float, img_center_lon: float, image_width: int, image_height: int) -> Tuple[float, float]: + """ + Convert local pixel coordinates to geographic coordinates. + + Parameters + ---------- + x: Normalized X-coordinate in local pixel coordinates. + Range: [0, 1] + Example: 0.5 + + y: Y-coordinate in local pixel coordinates. + Range: [0, 1] + Example: 0.3 + + zoom: Zoom level. + Range: [0, 20] + + img_center_lat: Latitude of the center of the image. + Range: [-90, 90] + Example: 37.7749 + + img_center_lon: Longitude of the center of the image. + Range: [-180, 180] + Example: -122.4194 + + image_width: Width of the image in pixels. + Range: [0, inf] + Example: 640 + + image_height: Height of the image in pixels. + Range: [0, inf] + Example: 480 + + Returns + ------- + lat: Latitude in decimal degrees. + Range: [-85, 85] for Web Mercator projection + Example: 37.7749 + + lon: Longitude in decimal degrees. + Range: [-180, 180] + Example: -122.4194 + """ + # Get image center in Web Mercator projection + image_center_webm_x, image_center_webm_y = geo_to_webm_pixel(img_center_lat, img_center_lon, zoom) + + delta_x = x*image_width - image_width/2 # (4,) + delta_y = y*image_height - image_height/2 # (4,) + + # Get bbox center in Web Mercator projection + x = image_center_webm_x + delta_x # (4,) = () + (4,) + y = image_center_webm_y + delta_y # (4,) = () + (4,) + + # Convert bbox center to geographic coordinates + lat, lon = webm_pixel_to_geo(x, y, zoom) + + return lat, lon + +@jaxtyped(typechecker=beartype) +def geo_to_local(lat: float, lon: float, zoom: int, img_center_lat: float, img_center_lon: float, image_width: int, image_height: int) -> Tuple[float, float]: + """ + Convert geographic coordinates to local pixel coordinates. + + Parameters + ---------- + lat: Latitude in decimal degrees. + Range: [-85, 85] for Web Mercator projection + Example: 37.7749 + + lon: Longitude in decimal degrees. + Range: [-180, 180] + Example: -122.4194 + + zoom: Zoom level. + Range: [0, 20] + + img_center_lat: Latitude of the center of the image. + Range: [-90, 90] + Example: 37.7749 + + img_center_lon: Longitude of the center of the image. + Range: [-180, 180] + Example: -122.4194 + + image_width: Width of the image in pixels. + Range: [0, inf] + Example: 640 + + image_height: Height of the image in pixels. + Range: [0, inf] + Example: 480 + + Returns + ------- + x: X-coordinate in local pixel coordinates. + Range: [0, img_width] + Example: 100 + + y: Y-coordinate in local pixel coordinates. + Range: [0, img_height] + Example: 100 + """ + # Get image center in Web Mercator projection + image_center_webm_x, image_center_webm_y = geo_to_webm_pixel(img_center_lat, img_center_lon, zoom) + + # Convert a point to Web Mercator projection + x, y = geo_to_webm_pixel(lat, lon, zoom) + + # Get delta from image center + delta_x = x - image_center_webm_x + delta_y = y - image_center_webm_y + + # Get local pixel coordinates + x = image_width/2 + delta_x + y = image_height/2 + delta_y + + # Normalize + x = x / image_width + y = y / image_height + + return x, y + +@jaxtyped(typechecker=beartype) +def xyxyxyxy2xywhr(xyxyxyxy: Int[ndarray, "8"]) -> Float[ndarray, "5"]: """ Convert Oriented Bounding Boxes (OBB) from [x1, y1, x2, y2, x3, y3, x4, y4] format to [x_c, y_c, w, h, r] format. `r` will be returned in radians. Modified from `xyxyxyxy2xywhr` function in Ultralytics library. @@ -381,20 +306,51 @@ def xyxyxyxy2xywhr(xyxyxyxy: Float[ndarray, "n 8"]) -> Float[ndarray, "n 5"]: xywhr: Oriented Bounding Boxes in [x_c, y_c, w, h, r] format. """ - if xyxyxyxy.shape[0] == 0: - return np.zeros((0, 5)) - - points = xyxyxyxy.reshape(len(xyxyxyxy), -1, 2) - rboxes = [] - for pts in points: - # NOTE: Use cv2.minAreaRect to get accurate xywhr, - # especially some objects are cut off by augmentations in dataloader. - (cx, cy), (w, h), angle = cv2.minAreaRect(pts) - rboxes.append([cx, cy, w, h, angle / 180 * np.pi]) - return np.asarray(rboxes) + points = rearrange(xyxyxyxy, "8", "4 2") + (cx, cy), (w, h), angle = cv2.minAreaRect(points) + rbox = [cx, cy, w, h, angle / 180 * np.pi] + return np.asarray(rbox) + + +@jaxtyped(typechecker=beartype) +def xyxytoxywh(xyxy: Float[ndarray, "4"]) -> Float[ndarray, "4"]: + """ + Convert bounding boxes from [x1, y1, x2, y2] format to [x_center, y_center, width, height] format. + + Args: + xyxy: Bounding boxes in [x1, y1, x2, y2] format. + + Returns: + xywh: Bounding boxes in [x_center, y_center, width, height] format. + """ + x1, y1, x2, y2 = xyxy.tolist() + x_center = (x1 + x2) / 2 + y_center = (y1 + y2) / 2 + width = x2 - x1 + height = y2 - y1 + return np.array([x_center, y_center, width, height]) + + +@jaxtyped(typechecker=beartype) +def xywh2xyxy(xywh: Float[ndarray, "4"]) -> Float[ndarray, "4"]: + """ + Convert bounding boxes from [x_center, y_center, width, height] format to [x1, y1, x2, y2] format. + + Args: + xywh: Bounding boxes in [x_center, y_center, width, height] format. + + Returns: + xyxy: Bounding boxes in [x1, y1, x2, y2] format. + """ + x_center, y_center, width, height = xywh.tolist() + x1 = x_center - width / 2 + y1 = y_center - height / 2 + x2 = x_center + width / 2 + y2 = y_center + height / 2 + return np.array([x1, y1, x2, y2]) @jaxtyped(typechecker=beartype) -def xywhr2xyxyxyxy(xywhr: Float[ndarray, "n 5"]) -> Float[ndarray, "n 8"]: +def xywhr2xyxyxyxy(xywhr: Float[ndarray, "5"]) -> Float[ndarray, "4 2"]: """ Convert Oriented Bounding Boxes (OBB) from [x_c, y_c, w, h, r] format to [x1, y1, x2, y2, x3, y3, x4, y4] format. `r` should be in radians. @@ -405,76 +361,159 @@ def xywhr2xyxyxyxy(xywhr: Float[ndarray, "n 5"]) -> Float[ndarray, "n 8"]: xyxyxyxy: Oriented Bounding Boxes in [x1, y1, x2, y2, x3, y3, x4, y4] format. """ - ctr = xywhr[:, :2] - w, h, angle = (xywhr[:, i : i + 1] for i in range(2, 5)) + ctr = xywhr[:2] + w, h, angle = xywhr[2:] cos_value, sin_value = np.cos(angle), np.sin(angle) - vec1 = [w / 2 * cos_value, w / 2 * sin_value] - vec2 = [-h / 2 * sin_value, h / 2 * cos_value] - vec1 = np.concatenate(vec1, -1) - vec2 = np.concatenate(vec2, -1) + vec1 = np.array([w / 2 * cos_value, w / 2 * sin_value]) + vec2 = np.array([-h / 2 * sin_value, h / 2 * cos_value]) pt1 = ctr + vec1 + vec2 pt2 = ctr + vec1 - vec2 pt3 = ctr - vec1 - vec2 pt4 = ctr - vec1 + vec2 - return np.concatenate([pt1, pt2, pt3, pt4], axis=1) + pts = np.concatenate([pt1, pt2, pt3, pt4]) + return pts.reshape(4, 2) @jaxtyped(typechecker=beartype) -def _get_covariance_matrix(boxes: Float[ndarray, "n 5"]): +def geo_to_webm_pixel(lat: float, lon: float, zoom: int) -> Tuple[float, float]: """ - Generating covariance matrix from obbs. - - Args: - boxes: A tensor of shape (N, 5) representing rotated bounding boxes, with xywhr format. + Convert latitude and longitude to Web Mercator projection at a given zoom level. + + Parameters + ---------- + lat : Latitude in decimal degrees. + Range: approximately [-85, 85] + Example: 37.7749 + Beyond the specified range, the projection becomes distorted. + + lon : Longitude in decimal degrees. + Range: [-180, 180] + Example: -122.4194 + + zoom : Zoom level. + Range: [0, 20] + Example: 17 + + Returns + ------- + x : X-coordinate in Web Mercator projection. + Range: [0, 2^zoom * 128] + Example: 1000 - Returns: - (torch.Tensor): Covariance metrixs corresponding to original rotated bounding boxes. + y : Y-coordinate in Web Mercator projection. + Range: [0, 2^zoom * 128] + Example: 1000 """ - # Gaussian bounding boxes, ignore the center points (the first two columns) because they are not needed here. + # Convert latitude and longitude to radians + lat_rad = np.radians(lat) + lon_rad = np.radians(lon) + + # Project latitude and longitude to Web Mercator + x = lon_rad + np.pi + y = np.pi - np.log(np.tan(np.pi/4 + lat_rad/2)) - gbbs = np.concatenate((boxes[:, 2:4] ** 2 / 12, boxes[:, 4:]), axis=-1) - a, b, c = np.split(gbbs, [1, 2], axis=-1) - cos = np.cos(c) - sin = np.sin(c) - cos2 = cos ** 2 - sin2 = sin ** 2 - return a * cos2 + b * sin2, a * sin2 + b * cos2, (a - b) * cos * sin + if np.any(y < 0): + warnings.warn(f"y-coordinate is negative. Latitude='{lat}' might be beyond the valid range of laitude for Web Mercator projection (approx [-85, 85]).") + elif np.any(y > 2*np.pi): + warnings.warn(f"y-coordinate is greater than 256*2^zoom. Latitude='{lat}' might be beyond the valid range of latitude for Web Mercator projection (approx [-85, 85]).") + + if np.any(x < 0): + warnings.warn(f"x-coordinate is negative. Longitude='{lon}' might be beyond the valid range of longitude for Web Mercator projection ([-180, 180]).") + elif np.any(x > 2*np.pi): + warnings.warn(f"x-coordinate is greater than 256*2^zoom. Longitude='{lon}' might be beyond the valid range of longitude for Web Mercator projection ([-180, 180]).") + + # Scale Web Mercator to zoom level + x = (128/np.pi)*(2**zoom) * x + y = (128/np.pi)*(2**zoom) * y + return x, y + @jaxtyped(typechecker=beartype) -def obb_iou(true_obb: Float[ndarray, "n 8"], pred_obb: Float[ndarray, "m 8"], eps=1e-7) -> Float[ndarray, "n m"]: +def webm_pixel_to_geo(x:float, y:float, zoom:int) -> Tuple[float, float]: """ - Probalistic IOU for Oriented Bounding Boxes (OBB). - Inspired from `probiou` function in Ultralytics library. + Convert Web Mercator projection to latitude and longitude at a given zoom level. Parameters ---------- - true_obb: True OBB in [x1, y1, x2, y2, x3, y3, x4, y4] format. - pred_obb: Predicted OBB in [x1, y1, x2, y2, x3, y3, x4, y4] format. - eps: Small value to avoid division by zero. + x : X-coordinate in Web Mercator projection. + Range: [0, 2^zoom * 256] + Example: 1000 + + y : Y-coordinate in Web Mercator projection. + Range: [0, 2^zoom * 256] + Example: 1000 + + zoom : Zoom level. + Range: [0, 20] + Example: 17 Returns ------- - iou: Intersection over Union (IOU) matrix of the two OBBs in [0, 1] range. + lat : Latitude in decimal degrees. + Range: approximately [-85, 85] + Example: 37.7749 + + lon : Longitude in decimal degrees. + Range: [-180, 180] + Example: -122.4194 """ + # Scale Web Mercator to radians + x_rad = x / (128/np.pi) / (2**zoom) + y_rad = y / (128/np.pi) / (2**zoom) + + if np.any(x_rad<0): + warnings.warn(f"x-coordinate is negative. x='{x}' might be beyond the valid range of x-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") + elif np.any(x_rad>2*np.pi): + warnings.warn(f"x-coordinate is greater than 2*pi. x='{x}' might be beyond the valid range of x-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") + + if np.any(y_rad<0): + warnings.warn(f"y-coordinate is negative. y='{y}' might be beyond the valid range of y-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") + elif np.any(y_rad>2*np.pi): + warnings.warn(f"y-coordinate is greater than 2*pi. y='{y}' might be beyond the valid range of y-coordinate for Web Mercator projection ([0, 2^zoom * 256]).") - xywhr_1 = xyxyxyxy2xywhr(true_obb) - xywhr_2 = xyxyxyxy2xywhr(pred_obb) - x1, y1 = xywhr_1[:, 0:1], xywhr_1[:, 1:2] - x2, y2 = xywhr_2[:, 0:1], xywhr_2[:, 1:2] - a1, b1, c1 = _get_covariance_matrix(xywhr_1) - a2, b2, c2 = _get_covariance_matrix(xywhr_2) - - t1 = ( - ((a1 + a2.T) * (y1 - y2.T)**2 + (b1 + b2.T) * (x1 - x2.T)**2) / ((a1 + a2.T) * (b1 + b2.T) - (c1 + c2.T)**2 + eps) - ) * 0.25 - t2 = (((c1 + c2.T) * (x2.T - x1) * (y1 - y2.T)) / ((a1 + a2.T) * (b1 + b2.T) - (c1 + c2.T)**2 + eps)) * 0.5 - t3 = np.log( - ((a1 + a2.T) * (b1 + b2.T) - (c1 + c2.T)**2) - / (4 * (np.clip(a1 * b1 - c1**2, 0, np.inf) * np.clip(a2.T * b2.T - c2.T**2, 0, np.inf))**0.5 + eps) - + eps - ) * 0.5 - - bd = np.clip(t1 + t2 + t3, eps, 100.0) - hd = (1.0 - np.exp(-bd) + eps) ** 0.5 - iou = 1 - hd - - return iou \ No newline at end of file + # Inverse project Web Mercator to latitude and longitude + lon_rad = x_rad - np.pi + lat_rad = 2*np.arctan(np.exp(np.pi - y_rad)) - np.pi/2 + + # Convert latitude and longitude to degrees + lat = np.degrees(lat_rad) + lon = np.degrees(lon_rad) + + return lat, lon + +@jaxtyped(typechecker=beartype) +def obb_iou_shapely(obb1: Float[ndarray, "4 2"], obb2: Float[ndarray, "4 2"]) -> float: + """ + Compute Intersection over Union (IoU) of two Oriented Bounding Boxes (OBB) using Shapely library. + + Args: + obb1: Oriented Bounding Box in [x1, y1, x2, y2, x3, y3, x4, y4] format. + obb2: Oriented Bounding Box in [x1, y1, x2, y2, x3, y3, x4, y4] format. + + Returns: + iou: Intersection over Union (IoU) of the two OBBs in [0, 1] range. + """ + p1 = Polygon(obb1) + p2 = Polygon(obb2) + iou = p1.intersection(p2).area / p1.union(p2).area + return iou + +@jaxtyped(typechecker=beartype) +def obb_smaller_box_ioa(obb1: Float[ndarray, "4 2"], obb2: Float[ndarray, "4 2"]) -> float: + """ + Compute intersection over area of smaller box with the larger box. Here, intersection is intersection between the two boxes and area is area of the smaller box. + IoA = Intersection / Area of smaller box + + Args: + obb1: Oriented Bounding Box in [x1, y1, x2, y2, x3, y3, x4, y4] format. + obb2: Oriented Bounding Box in [x1, y1, x2, y2, x3, y3, x4, y4] format. + + Returns: + ioa: Intersection over area of smaller box with the larger box. + """ + p1 = Polygon(obb1) + p2 = Polygon(obb2) + intersection = p1.intersection(p2).area + area1 = p1.area + area2 = p2.area + ioa = intersection / min(area1, area2) + return ioa \ No newline at end of file diff --git a/garuda/od.py b/garuda/od.py index 01b1b76..43c15f8 100644 --- a/garuda/od.py +++ b/garuda/od.py @@ -4,9 +4,9 @@ from dataclasses import dataclass -from supervision.metrics.detection import ConfusionMatrix as SVConfusionMatrix, MeanAveragePrecision as SVMeanAveragePrecision +# from supervision.metrics.detection import ConfusionMatrix as SVConfusionMatrix, MeanAveragePrecision as SVMeanAveragePrecision -from garuda.ops import webm_pixel_to_geo, geo_to_webm_pixel, local_to_geo, label_studio_csv_to_obb, obb_iou +from garuda.utils import webm_pixel_to_geo, geo_to_webm_pixel, local_to_geo, label_studio_csv_to_obb, obb_iou from beartype import beartype from jaxtyping import Float, Int, jaxtyped from beartype.typing import Union, List diff --git a/garuda/plot.py b/garuda/plot.py index ffcf89e..68618a9 100644 --- a/garuda/plot.py +++ b/garuda/plot.py @@ -4,7 +4,7 @@ from jaxtyping import Float, jaxtyped from beartype import beartype -from garuda.ops import geo_to_webm_pixel, webm_pixel_to_geo +from garuda.utils import geo_to_webm_pixel, webm_pixel_to_geo @jaxtyped(typechecker=beartype) def plot_webm_pixel_to_geo(img: ndarray, img_center_lat: float, img_center_lon: float, zoom: int, ax:Axes) -> Axes: diff --git a/garuda/utils.py b/garuda/utils.py new file mode 100644 index 0000000..4a5d78f --- /dev/null +++ b/garuda/utils.py @@ -0,0 +1,27 @@ +from jaxtyping import jaxtyped +from beartype import beartype +from beartype.typing import Sequence, Tuple + +from garuda.box import BB + +@jaxtyped(typechecker=beartype) +def possible_intersection(obb1: BB, obb2: BB) -> bool: + max_lon1 = obb1.max_lon + min_lon1 = obb1.min_lon + max_lat1 = obb1.max_lat + min_lat1 = obb1.min_lat + + max_lon2 = obb2.max_lon + min_lon2 = obb2.min_lon + max_lat2 = obb2.max_lat + min_lat2 = obb2.min_lat + + if (max_lat1 < min_lat2) or (min_lat1 > max_lat2): + return False + if (max_lon1 < min_lon2) or (min_lon1 > max_lon2): + return False + return True + +@jaxtyped(typechecker=beartype) +def deduplicate(labels: Sequence[BB], ioa_threshold: float, iou_threshold: float) -> Tuple[Sequence[BB], float, float]: + pass \ No newline at end of file diff --git a/lab/test.tif b/lab/test.tif deleted file mode 100644 index eea968d13d320cc5dd164abfc4b009bd41b80427..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 64902 zcmeFZXOvslbv8`f63J;8OhNCx_udD+_c{Ou7_`9v7)>=er|H6yrEvD?h%qxEuis_wIXe@4-3eK6jse_I@_( zcIJl78#Zic+OXjp-`emEK$-_W-?&NsCLn(cDEkojeCvzy_477-7WydHvj4-~hCx{&V@q>J1wfm+Cic1j_#n zkRMq6)&@Kv(={75Jhb}F4Ol?#2IToWzq#Qp)rJjV;MxU1UT+H@_^iJ+fBv4wef8e( ztrg&!jU5{{tV{#@lfeEeu>T>j|97CzpTE^sliuYn3}^D0;f%Hxk5}7=3srn?rdrSt z_1+Ls+dVwoKcwXGh6=sIrQSl1ut%8T3%J>yzAQJ_H(20~m3s4iV?*5D!throae9P2 zfs_~IS&cj)Unt`71w5gc7ZmE#?x5WT?AA;90--o0kP0OdLE6I?^ZDrjuuTI+L96l0 zGo&GbgfHXExk0_f%(Djcp0LZzTd(4a^Yna?LMoPjdH&}$Dg+{)-EB>q^dUVjAmCSd z|GCl|KMV9R4SdZVE_LaKhBZK|#2mhyv)(Q~MWh!=fsI5Kr)H}Cs`bxM;j42bA%Ot+ zTDZRXPv`V&27A7`V)TnEBpkl@tI|PD?yG}JUn#Fq5^gK5}B||A(v;nIC5#1fFl-Ux;P3!Cd(1#Ww|_GAjoIMvYRit zx$S)6dfP?7Mj*VYf0ZvgsNnDgq5tlrAx-wCdjE0{OGQm)BwUJmD|s)^+vq6Q80u;)=M&T`vLiyhY3jgXFS^l+Mc#ES@@wb7( zzv?ekeC;p1wNWVe8}VBIwbhA0@b!#*i=*&w#B2RmQ7HTxg|{{e<@&#k8TqfFP`)10 z{@PP>>w0R|84g3kgQb2>^$WhkEeFbVwn8uzu={wXVA>oE*a1j=&XTy!fe`T1!JBCl z>m&zv$don(x#QKIujvcx^o(21L0BBj^p5ltZe7s*4;VG;U<3 zVX)afL7zQf{u)BJEQI<83!|mN*ey3V{twiY-+rC<_uo+dd=dNCM5|jEe|`Cp-0-bC zU#?Gr*Ey_TYq47vMT3Rjd||MVPyZKn*58u>n97;KiCkZ`e`L5Y$gP%ggMCANUBdvM zxysA*mw4+-1>UHDR|4p$yf13t_V*TbxxSvBLTQ-x&V72(_88)3=Q9MTK~(2ao3yq^$G`g-fq1WPSEJO<#6*quok{b z6Zk9c9iWe2yLWDBC*yFbr%)d1>-`Ff|6@)6Wfk1@dj1x5{cx!YD0#p|N7SzSk$?mb zAkhUckSc`|r9{95a$op-K8w#+^7(&xvba7aeKo{iPJ(}9NkJi_{p%qAE&a4#2V%N= zu+a528vi`~?8^);AdT(o6t(MFdzJrRQ`Djn-`_k6fj)k%6K^q{sOJ6W^xf;3z#BG@ z+pbGr?!PSm%ljW3rRvweIkesN`F>qq&ko+Op1*wC7n#Ce?7w;2L#xpB?>9fRx@iM+ z!$u&P`WqW|lQz_^@4M%2`+UzJMXf=+T#V=8!yz|*vah!QSQ(z_>!s5h+I?+39#c7H z&!pW+mxjZ^F-ZukN|5!K6K;LLCgf32917gmCJnSncmxCv-oz!s`4kiz3HAFdRu#Re zz6M;k84Rw6Hr0!1w1{01vPxrSwS-ACE9fRMQ!Qgh&7!Q!urS`YYxmgJa(g*!vIr?o zEw3E4?x;jpd;KdF|Bg~`&mPJA0 zk*x$$3!=Fm32#BcTU(&@PF zc5!xoVQOi1c7A?tYHD_BaeU9h{8De>uD-}4XO~YFgLAp4Sw^GdS}7D1jfvvp8x>@f zgba~UkxC{;L_!k~bzY0CGokazXc`&5*x~VbGX3)s)8;>xnP?t`B(e1lbs2rh{+@% zIRa9`qloCFE(wK(T_5dIK$UUH*hB~m+pJ;WWfU|S1GAYqW{I=~RM%8n3%zYKtP#`z z-++K`Q1GBxugR!om_%$KO9Uc;xi#bW-2r7e}Sc&NkMTu<}&5^<1y# zjse%n9^27^{cy=TpD+g6gggd9$ir!sRIi=mwM!UO3>4l##k61%EhtzMnbg81Bg6s% zkB?*VP+SRCE+sh)T%nMHLN?;bO*AHwK|_%-t!O5)O@e*r&Ck~Uo7acGdE;!}^VnTW z+Xv&Bgj&R4+H|VXT-R)6bYO5aRZMt&Z8{N!L2Aar8yL8zV7qK6sV_L%yjo7u>u|a_ z0v@?1<{Hm-l#4cxSQhd~=p-Bp-3)K4GqP>xe|hEDmYiPEMr86)NL&*zN_Dq4)`OZF zo5*zhoW;1Uz0jZ1m!g_V#M05@j9NsYHsx^KooLrORYW|tRmMa6<{_*G5-ucC)9mg$vTEM0bndwP2D{hn-l?(t?38!-J zQhE2zido7K5t|)6Qqd;d71y3B7>^gM$BWi8J*Fdh{oYL5REM$MDw5MML7Q;@X!qIO z(+i_zi%Nz=HKS2Y2nYxX0TEHnVg_8v$MP6x0T*T0@C$BfkB2K25D3IpDh+ZyJn{+qG?qWJ0ayRKc>a(pRgmeTDe@BI&)f9Jz@-u=O|FF*Of)qBq0fAP%J z^q5Z0Ct;8xA-)h7cbknqqoUpHNT+RaGmpbW14B1nPL1@Yl0m>BYm*BZEsfx&#zs7e zF}u`vVR1pIv>~ym)@E2;T}|EQjWwHV8tXQ-Hi18H!B4i9eZ%eXh})KnI-)LJw@wl>_(!oeddQqsmLI9P=?R}Buo;k95;?Tm8y`wv3`-gjyUYBJ$;W;?z>C_pmcD3JWt7i3;aGOcPsAi)JJ*la3 zM_(aS?1;H_h9*RF1GupjTq_l*qs{<~glnwbSW~m9cGJe1x|)XC`i7c~_*UrWjrrMf zB|BV6Ce!}e9{*fd_tt9X{a41jgKa96ghp({fJh7l@H! z0bA7u|F1n7@Npv=Nl?exHWZvr#w}Vw&a67PW3qRwi->J$ti25ZsiPu6m{u^ZwGIv0 zOorC8V2utgaW1BvPuqv1#*m3qbhT}p?LW3-X>oZt9Z^~p3?>DNh14L>O%yzeM?+EZ z(B`HZaN}kuqLz+rC1WvM3PMGO6H#C!q7IJ)Ba!t4D%=@1to`Hn{%a#|l>FP9yYAh- zeP|e%02vp~Uw+}+XYM$(vUSVC_}sRIU5B?U?K!lxa$?8$fyMFTyXTLrBw`v(JY|VE zbbZOr>7HDFR6ALWO;tOGD)DNEyC)qBcziex4uGs<<)NK3)0!82$<^TT6hhMvN=&sr6xzS8-B|Os?o-4&xW-4!9o!nL~$72=_8N|im zDoOJ|*yA=RI%3*l)EqR3H4=uBhqK4blWxV0j@%Gt!qsF#mD8}g!{(8)FJY7v9r`0X z1_r7gDiUK?E5as^)S`g5)PXl`1l4T>9>&Io&4`wbZbnPhqn%CL$9g?S7psfogNsAu zV_OHyonD=Y${`>~=w>Xg5siS*h;RbF6#}UPLu;B4HE?(n5s#Mhuu3+Tjs>A0^*CfJ zp3+3+!3`REIB2M5+ouZIT~iYW4zBD!y71UNJ8rc0jeED>e|kqx#XA;pJn`6*51e1# ze`x2V+&{)kH}q$0qrpxP!Bi2=os8bF(CZ>u+JRbgY?m+a2OBD;FFKl9x8 z_AU(U+1h*J@Z7>kda^UPb7^k(YH@YC@W`&p^k7$i&Pm6i^b*SU@#wDE)Ji^aVzhhj zT)B|e$L+#^TR%AxJd>Td(UBWMAf2cLq;pAMPcArCjP)mj{TcJAT?2EKs8Fx|VCNzRwW+D09?{qYg&@fUY*25Dm16^0=W=&&JZ)X7WC{W2(n|FoKY8-v-buSo zsh5&DBn*`XBVZ5|0*pWaV~~w7M12$ZwuZ*rLCv)cObeF+;2(s7fe3)+q_j}!5QBi~ z)!H4QxVPP(2_(~zN-Eu3YyI-IpGJ@aSC+ zKJoNxKlW+>$!!VrP1!>pmyy~ zYiC}2ZlyO-Hb(-nczf72ohm)K=lt4#tgZd)+WD=++p5Vft2km)#%*oeraC^jdgrNq z10&VWxJStZEJGrcgojer^*n@5fKo86(AL_fmfPzZz6pWcP9THOFc2^)kO{4H3Jf?$ z&8J$FI+xAkcLv?wXe8Di3nfw>^IcP!pI@4K@9xQOpBj1c?3S?}%N>ime)osJf8|GS zfB4BSZuI4b@WC4|zw*Mlv-`#)ezQMp@jG>4bK7)ZdZuI4lHJeEh$} z*Yz4MUpllh(0OpYIG(V~bSHO?RE}-yJ2}^Pc zadzv(ojV4f-8b~~@#U>OUB}OyxVa;LCcpCBi?3ciw0*G`7!|8t=+>*8HXZP^&ZhhW z-I2ahva2&-Y~$n5XjpR#5&=cPn(H>#fEt>>jm@>|H{dpWOA`v-h=gtub0IncPRVDO zBs47_?~yahUQ4wY92qGO_jGTsb{twMr5*ZsN^TM%7;JdKCO>&_@!?DJPh8${?}4Gi zL#5g2>_WvoT`Y`M977S)1N)ZWf9vcAuN=5?ajF{E`AjmWS`g5PeHum9=a}yguXF`Z ztX8(qBvMx0fqPay{r1s6>)q#vzv*qyj`SoBj%Fvb#?^A^_{oJMTYGm+R1a<&9PWsC ztV$shNyWA>aVRIf~J-n;;v2*)&&(0jWcjYg&>Wi}bcbz(MN6{u@ zo8-*2$CmV4Lv~|1>6z+_4s=F~`D`H`kZ{?Jjo`-GI#5k5xB*z3)@`n>tJ}P}rlzI| zv>5?z#-O0sRtUJcp{a2bxuL_4i!7C>99u0qbr0M6%{UK zBkdxZT`dopBttoOE@?>F4KqWriI@TSOQnQGXjTh?E_H`l-s#kLdGraVD&y7$^!yH^ z*ka_mM9iv9x+P{kyi}YVPMMUvkXPOtvMp~NJbh+hWhJ_`8krwR_Vk&eUU}Rl8H+jg zl|$8#xTnK=e5pE;^cCZ_;jG_k7E2j892$TSh)+Z?NN^mafzr}|Z?10zZEULFfWbC! zs8A}r4gwfO4YxHl)X)j77KOmAquOLzhcV_4CxelA$Q21klO6GNN3<9WM}}j*lQX3Q z+lO{9SB~xK{rNrj{>5AHonKsjd+Y3{|GHKf32F7>;l9K`cQ|2lmD0ZXUT=SAxX{_2 z3OU&{T0>)fT?1g{*VNWEGyyg~6b6950R)0HgK;e|wScPEsO&1;{_Wci?%VFxutiLy zgpQ$On^|NelLD-EF(L+;On_lv&2$`6gh$dz5IoSS=9=c_1_G>!0NJdd;r2{NdaXK_ zfv1GjVuI+ z(+cKbp=^8$2UjmgVC4*^oJUOs4O=_?d%K-IZh16d?1>vgO1hQ{_e)qJ9tN0IAx)c* z%{2r#7!9vOfB~Ti3vGmeHaBkmR$cu!V2zvM4I3e#O)ZU^u`OVK8@&>;2aH@`nQ9BR zbsCKYb1IT5`ZAuJ+n4bbZ0+qqk9=91exR zTJ>rRfiL1va9*9YyBge<9RQ}ffBxZn4y{m1LrbVo5@a(9R)cB;!CRYwZW4)53KVWt zu)}hOkc=@gC<;1RN<_8Mh-N;^WfB~pNcsF~Dj7tD! z-alQ=t#rD^3wE`TOhzLqNI=y$;n7VjJQz?ra7Y~t3Pz$@;Ajw{r4H8Ih=qZ$kOpK6 z2;8u_X47q;y4zZtHWQ#=U|#lV1eXtwzxwdOi`)8Ny7<8RAO7^Uk6!%Yt53c6#52eD z@15U1v1`1vZE$|~($cy8iwo6`Y$n(fm77JZm67p#o;dUL<7W?FoVf4o*vChAFZTQY z`2DZ$M(>obHx!E7Re7r|n*a>F_UZ0eKH^9Roo<^Ni)^XiSO?fUzzqP*f>tnuBxH^+ z?~}yr!0Z)EMf&=jm00-7eWRx;DTe|0=QpyD(9}jJ>GVqAU*L>NslzNU4<=$R|7^v7 zs4HrgQbNE4;kS$@T%(=#at< zSmNPeEudPRn3DAy;y%6Ktg$QjRx2|mqwkm;eD#sNFJC?O`ZGsA{@ycdYrntIw>`%v z$KCQuR6|6y0Be@UP2Z~Dcsrv+Gq92)!XdbpMpSD< zbL-|tU?J98$6z46CVnBP+dbL&%XiMN{q}{Qe0cAF{QQBn|M&jd+E3Qje)-3>7hgTM z^1z9?2kx9&UC33FmTi66!LEecC=Rv>WO`=6tVn5D=XWnXe17}S-g@fC??2y@H{EFA z`@eW?Ht!fHrLt~w$Zik#R7H<@bf7a6btPgp;1Q6}X?6AWb@dxJ*VNZGfWTlli$yP# ziyQ-)fkJ!ZhHBKaFkXD&{2r~GY7#TNCYefrGfUXx-R*Njk%f|Hp%6SY+P`bQ|I9+_ zJm7+I@;WS{*+Tg8%GiCo=eMpFCzIBx-ee`hIk>Xz_2*8#`S`ilo;m--Q%B#pcJldq z554f<@qNn`i=2hPgUllON~wLmI~8$i!**WG#J92p+jmsHb7jZtkL`c|nX@1N=-Ow$ zdh+{E-1+ej-vIh}_UvSrMz|2u(*bW+eNFx5%?+D3*Mn-oO|_u<&7g+caV=nYOFgu1 zV^c$IYyIu0rcDs=W+=1{hXWCrV745dOKHyU%szdz`tDPE{_yd+KYeum-#<9>FCX9W z+aI6!-P?!P-afwe(WTEmJoNUnyFa;h=Fs7ih=B|nbR%PhfL4&!OCmnYj^*mPt-F5z z$wzBzzq&E(zxnvpwO_rxx>%X%NlsQWMQXf}%ihvkT%7DExa1uUd8egKtCOBN zSpME)XI}mG@$Wr*?1QK8diUvzmo6Us!H?FRNO$bolF|u|_lCo5VsHbf9<-^UzOETu z*HCkNYjYD81t?PxhEi|R$&grN-Nude4I67VehbmKS;(f@Ly}xbasEj0#kWb;lsN=eRBUFUcd7np5J}n{)MY|pFX`(~ zNx`y6DONUS)Gg>Sb0|bGo!D$s(kl_gsE_Fov)m3NM!S1+8M%43##~1qFyLRg6_Rgr? z>~$$zKHg(@9sT|P`PteZe|K{S0r}Uz{&b|AOxpy9W)j7)0f%m-(8&x2mW{{p=)80? zaC~^Cn95BTW5*|sDx3xZnMS8jnKYhA%r)9|)n5C#>GYMmw_Vu2a&m6q-mQguX5x1( zv|riL`P7c;{W~hJ?C*W;$naYy2A(-ExPLr(|K9THj#|95Teo10#>4d!_Xw%C?%}_{fEvOCzYHV(; z2Z7fGI9evAR~_DKj~%;z_rk}&dg9upv6)W5RHW;LM3;uwCSx1;SPc!?0^isS2J8or zoDElUU}$6m9t9<$AVf462LTzCvd2$P9UbgyGq8_s$sAl7xHvcP^qIrYK6uBAubqGN zvAqvn+IejI;M8bb&%`J=qSFWa-nj3;NY>_aYLY&gSHa!6IRD$9e)#Wy_{FFH^1-h^ zdH1)Ue(>36AHDa3rw<+&NqM9*QTx7uaMHoEve=6=W2;-H(y`FY*vQKdT>AM3uYd3T zXaDJk&wuyn=O)K{#``LnNVt-zZkd~%n;cl3&c1YV;-0(DojrW#@uBh4!~Q!*{P!$H z?;Q(2xm0*c-GebZS94d@4~r#2QyJ`N7NKCTO(1MPp1hd`~#^#Ip^t2Svs6Hugyur z;Ttvq&u_!J3IYM<0k{Q%;?ihVxA*Mni5Cy79$OgOv636_?ufY*AvM!3V>onN83#v1 zG{Yf{K#YluhcQqHBbN&ha_Be=4U1(l5gZCZC!j1AVhaOtkAb{0H1p&$XWsbXgCD(h z;?$0z%R6WGZts8Kz@aA|J-kv08iZ5{lQ0yr4doCuCIob1v_^(Y?z2>y;71-EkAvB>AQDbedwu2-+blj&)@&x_y6+4 zwV!|dhd=$`?_PW0!1n3$`*s{Ze)`S}yPvw_jy+2~_w4I`_saIG51xMF^2N*NcRsQ! zd*4j^+3~>LBhhQinU{8TUE5uJZcpXfj_#*d3(xG#J~)$D9kufr2nC`c$yB6mG485VoWFhR2f!`()7sikfARLtT{H2hUJtO8_$-x>>CsAa zamRsp`*eRWmvnS^?3s|qBvLgsH#Guo=DJ$oK>#L$hPnn=V?7Or)hO6f6^%!S%JG;| zIAAjf<8ADC8^0Y09oiLaIs$O`Ta^4rP#dvo=Zfj6-eM|h%S4QYfFl?3^%g=qhx_(# zOIP&TD-S$+0~m8-Q*+Cmdj^W%d*QKv`RM}_(P+XcZsRkH0sTa~zZA7ib~+BuRu3;! z_Kqd{@*V9SN!2dz^R^9!wPlwmXQa2QIAJAmJYbxRsG>@0yNZ5vVCbc%9{SG9_q_A? z!*9NP?cEnH{nK-gz4+>h@7#Cr+?hjrchB$Jv$%iP_@&c(w@mciv90&*SNA@0-?@|f zPhZ}>?eX3H$HtOV*|q}%q5D@lA6?B{S;{`K(*5YR>?2E^4=g1w%_sLv*o1hbl1K7+ z+gwVq&#fzU`t!MD$gJp$*c~cvvcomh7feOX#Y$qJ$5jrR$MeC3;bg+2!y%DScpVDb zq!3W6asT$d$lU%=Z!%UXCI?5;xs)4^Mgl7_z>28d2-rg3s@w2Qz@DgYs)5~p8yZvt zK{PiaY7}&mjD<8wncg;rTgQvowKgSNO2<1@?1_YNxI@3Uw`zeG*u$c$9)dSU()J&So6jE1}7#D`S+FtcsmofxQ{) zSlCc;Nz-P2-mUEn7(xMYHLM;@*mn#k_bjF^@0{+hCiQBePRmg#$@z%I8L}7R);AvA z@|`RD&+eExI9u+|dTl0&l7$Y+8U23s_M+?bSoprF*rQ9yhZfWK&m|t-mAPXoLMJ0> zOgLbgv&ew^nxJNq94d}Q$u{d)3MSdA9GLl!0h?Xe8OcP^cu^hbqx6)$~LmC=)a4TUr|Hz>wOFTpG*owkPyfFjUQxv9~8llcSliPMvvOTel@&Gu;7J}p0C6Zc0< z^F8k4^FtSo%zX0nl`XR!RykY6BRM4$lNc9Kkb3Nja>$srsm7zeW0R4UtffC-3F?J% z5rRo*Cc+!>kOnpZuHq4#N=C#i&4kR}kde=%({V%`2}8g`2}FdDOI{pFJaOOh1N*k* zI&Ca6gouE{5ilqUOhkhe9Q;JudS&Ox%lEE6dvNgTw(?t-MlWycVhdm#TC+pOC6XXK zGD1Yh>G(vCjF2};eJZX+!6-Uy7Q0kHhcoGT3Zj`$K-pz@+bFO^0ssh~@RKw(i08&*b+Uwi(c=ih$!^z%m({?5_9^lvUdcU^A) zPY@hl3oLknKpqNNk3iH(S@=XiH`rmHOj>%Q>fV5QFl6Y7sRt9rp@emy)7lf1jV1&0 z6?3myFzA$y`hkZ+eYnfN*x?%Xs%KNqm1-oPv+B4+o0K2ai|uNLPsZ%?o2R>yr?wBj z_4H%8tQge*CSbvUf0)H+*3yuADpteAYT1~AiN8N>^J{oQ7L<)^=6>!HApi?u0)mHz zNvWtdE?&hUsM#c&M=9kn`E0sXDc31v3@Qr?gLn-ZNh)Bkh=x8nqB103! zJONWq$m@jn?H+ve;-NeDEuB9y^2C|S(S?v%g;#KK1{s@%hw^DSHH#YXI2V?R7w7vY z+tq#r7w{h@gHk00PR2F~$#}C-IP7a1cROMxl|?{uOIbmkaJcAS$aO9BBxWk^{&aY> zr)zS!D;c+v>11&;I*m4(@*bp2KV3Tc1tV)_xN)($fw>&M)co ziSiC{&_D+Q{ZbLJm@v;~f^mbur6L4Xw6vPlsbdu^(!qdnCL1_7H?U_aHI_BaX07ub z;e~?t=ydnCv1o5caI!PCt-Gt#X;RSeCNae=r3Yjji-;LAOS1v(Y|%nI zN08IeZ3>b>B#b6J9qoab(-8)u*A_!(*mPuNcHf2h(UfnxG=KGxYw!K+qqRT%@;can z?PqJ>9V$BG2~9_?t*fLeW~5RTo=t*?c?2L}NTr~Yap&Rj%1g)QUfer;sv0bqapr%ue0D)FNw|@P-SKof}kq1w{{>a6Zq1djh9hkLOx@@Ng1NV-GE{umSjsy>v zTz(USfUh?yXamVWF=>oznC)$>tdpO$0C8j)5J{ZLcrUJOJv~3M+7+5FI7d>(iFU(s z(KDD-<)V>vyRE~cPgz-N9I}mx)v+;l9?r@pgxZ9iQCT5i*tOK9u?(w}eZ+d?zeAk@4$0y1;L?J0h!N%h$u%t!u<}3F+ zHCFOi#D$P%zT_E78k~B@SSC0+F*wsdHQ!$uD#XUB?aKp&nA?seq9|~@$*NA-EQ{s* z>_E3KVALxZKvb)#X(Ot!ZfCdUo|)joi-|{&^N-5x>3~b*$Rb29j;(C5U!%Lc#Hm%67B%5jL zG!-+4q+#K;bVAB0?@hUEda+I;q_UfT_~dJ!kIBDYADMk~al3^R*0bzVQeRZnu3=zc zjdCWRsv>`RZJX+Sq7r|m|2_;8M4`^vqOGi<&J&l_s;Jqoj9;? z?9k#p`=+Zs<5yq)%bsLjSi-Vl@oEG~4aMtVC=M8dZ?ee%D$RQAV5aPRW!UoRtopM> z`yIoIh?VKp(6l_NR!SnVAs;Uvx!z!-iX3+<3+;A~ho@x|I|J5{a%rHeINlrY%?AcE z;n|+N*)DErt_N61O1o60W&&W-a7ch61vmqsrUt-8(b(E@pklv!BKYWB^2(Nud&i>} z$HI4xhNsHyS{*Uwmmgm0KfF{O%-Aw9P1tGf@agj&NmR*+YU!g1?IZigA3MLeFyzVD zC4*Ue#@!x|>D&D(HIL#D5=?BooP$&G0oDXmfI((s%DFsni`Ya8P6}_m8IXYd&Yk;` zZk0tv_G<(kR*^+OWs*^qWMKd5_RdNy<8@9IO2tI5O+>K@=)(#1LSlWfkhBT7RE&y) zy?gb--V>8St2!Sx4`w|F$9fZ8kyoF(|IdEq`k^>aWE%-Y4Fao#qSdW1I+#F&m{dYK z0}6PKsT@ZwE*`-%?9v> z4ZFJRXC@*Srb1^Y1E?(G!6z-IiEt$rt9TxLt~k-V%(|{Mh(2NkZRXxPj}V+=9QZr_+Q-$Qvyad5J_4j zPTh)DK_N6S4%Vnvaq!qCqmnq^(YBoA@6QQVvV6N8PZYp#47frkn4eF_41%j2!L{%0 zSUWS8)^lWY4R$%fpyVdQ z`jK*=obZ=i>YPUsXd_xhB%g{tn|7ZV@|FAso0g%IP_-PCoQ-tJX-);jDWkLLa15*! zi-s~W$Tl&hppZ|GClh`RmqJULB3}*1=TbOiA;4RKrm0!*Zi`98Mgsn41&eG{h&57r zee`LRbCnx9UBLTeD0Z;WB=B|Skm#* z)wlobSFRtLEgYhSK-D4fS}0o93Zd7-nl~$z6b23@CL?lAfyqEt>d-Pdn#^xSU_o$X zBa=ex&AArb!L=uMeO}2YVXLYu>g8I$sk z;!%S(-@a|jZ7MXbxrHetwMhgffgqPkojf{!V0j{BQvp^cVLkqc?9g$6upJc*HyC;2 zrTmU9(+8JEyOQCsQNC*>KN>4M`{uJZe&S|9YLro}Xo?Pr)j(15Rw%O(3ajI>aX2`V zPJx7s6sa6X56Wb1yt z1{AL`aL~MycY4Hkc+hjW-#i=FrJU@Fj*~J=IxSos7bPOJppam#gt0PQTI%nNSrj@E zD`8S4ZHmRXW1?ygDoAn~O3xqy8dOY2X}LI$hMf0G%v?H^Me`|@`*+>(=7$e`_ns4d znMA~`=5YyXGUn&6{=_l4Ht8maYcDg{*a|E`^JD3D2$!l9DtHGB!5-e%^G7Q2hdP&w`COZ$PDYoe0rZ1>Wr zu#3Cz1B{HH|L$Mk`q{gm{N!ifdHpqsRE@*JFKl1_s))T3W>+%jhU95`s->jt7lXqlxjHzuhd(8AQDSf1CRnXdm zI3}uIL2hv==%t{^(x290DK4e)09M97i_~&gb%* z4ho4#hJ|&ygWI+q+_RkNbQv^!9u>`|QR!4N6_4hU$PyllM?~mUO20){Om~b{Gd7)6 z$fOxXxRT#h=yASx?Xmw-^Lhb4iJ@zW7&#oLfMI12gs>Gy0kMc^Af^U|HSj5Ls{rc& z0@+j)2?nCjkVX-khDIfo(rCm`aET60TKAST79rj$Wh>tlM)B>MDXD|Tn3YeIKt+dB+wH05x@X4*3lT33$J&4H5OD3!{_xq~z4P(Mzx>%NKmC#3%7=hpSOOk~g77Jb zl0~(rY~Po+4h4lB7KTScwo66XsQL7G{UHkfj5v?%V#)#@gN)f># zVh>l-Pu;oaz~snaD%_h(4Rv>R=c0v>)2&yA9gYL56RTVLXR4X$(aO|NcTd7MShDXJ zExq#e(_anRSJEHc1|ds}V~Nlj6NM*4Qp9j38R4~YM06~MR8OLUjbcK{Z;3k8GC4&^ z$C!+=7w_HqRRwEnM-I&0e|&s88`Y_4G#r|OgA3_M4i$|? z0<2dMlYrIBIDs~)-=YclwPA<8!(kflaPIF5?wu|k+SA|Le3Dbi%YnG9J;DT$-VjW;EvW>A5Buzm)S{IKAu0g)Li`E34J~j=^$oU#u_h zF8l1AG3(-ZvZo_FHjvmdULNVpq$8efWA1HZ#XFA;{`22Gy7u<9-+uDSZ$G=b_N)8W zK7R7qGm}N1E$NndwZc-xp7tqyaE)IFhxbfe*j*VbdQuUG#VE4NMLID-M#p9Sy5;fiiE^^PGdftzkN0PmCdw|a zkxYOt4$b}SUw;3|&p&w44eY-#L70@BHbV^T)Q%_UBv?yEtxC|HPyJ4X0p3It;%Yps# z-P=ctnTW1Gs%y8aA{K?msk52+c0Fq$>mM&A=X>o-v+=&7xD?{f+9g9#WqW`fwz6Gj zw$(t47&zsaB^xpu%?h83H=2rcB|SVLiA9C6XhbQSr4-VEKNXRSND2|dXV8SKYQQru zqalDiESp&&V-qB7f{?{fE7%$iGG-Ic$HMbDUnOD~F2?7FG7FQ1P*6#SqKX4UAO8H~ z7yj;hkAC}!ty@;eEG`(_*w|3p*jNi~Y3cC@hC-6LxaL5=^!g^}3}}4qm{=i9`qiwnc`3b9opTy^c=0y~`&V%KH!REAJi|9v$fDbj!VJo<+k{ z2{6Fa$)iA(ENsH59?toz3HEP~JoTGvuWuhNj7Nb%Bh#fgiUb4LmK1K2REY4ZX>p6v z?b7I&gksz^UGdpO0xTYZ0&mdZunrC1X;mvE6ebR)=kl}~p-n?maxi8Y+o9#Li7+Ms z!5|=&EFvF|bSi{nU7o$CvXlAXbboG3rE7U0(`nNYV8rhJ;`cxL=-C%uyzmI+t_Afbj!<}t8TOp}NT?R3fd112pMH`$*3 z{-uXs`RKW4Uq1TxXHNe5v6p`EgLk*>Or~u7?IV$Pf7op1JCtI58`mOd7-VEI7YXoM zg$%S&MysZLyQgB;mG(b<^yKb7Jzs{$h@f~e6v=Icvm0UDMyVJXQZfLa5*D#*NYp8hS+y38kW7J7@JJP%*q`gnC9E-@4kh`oqPZXO;%SbLHY>HPLR@<8Tpf$*hS2;_XVMcfvI zsMUNX5zgZg0Tacl;Q0(JyM%0!Q*=U-LCEokW&Lv3^{amQ*~3GfGLj636+;lbW+WS6 zL^Z+LAQY#WBSP9FK#moVFsNJX>()8Ar;9^tzx&PE`RQ0x7dPo*4xe1r>6N`s6B7_i z1r9~;z z@W^TfEu-e8i~uwkH;wuFQt`{zuTS;oR%Y^YCdy?o9GabSDdidwCG3=jt@4U! zn*8X_;$@^n=L&zc=%m=5k!Ie_*&{+QJ`LCG#J6Ai)Zd%YrQ|eGSK6D-bTrh+S zZe=zR7+{3}OT$8KfCHgq?l6eO6tq>sF$>wiLL*^SmEFdmO{x$um}G*0Ll@GCF0DEk zX#=vHv?3gfhyfmUz|GjkCCvAQ4@{J2dU~dciQWDA9kabJ%#bZ* zY}XIA#0|5Zx`B`=XJvCptqnEbsH^!V46&I(MR|R~uv6C^0}?F-Ast19QeVI1oxcNm zz5crUT+t&Wm+aDQS(8P9rxKf$e27Czu<&t#Hd(vZKG2aGDFDFa+YRhe&~7rxT3{Pn z5Ou)R4*}QU;SdrY4&)cB`4}@>T8TS8y>#XH{Aekt==O-&WUSM}TRBj@Yj5%D@Wk3% z-v`$HYim#6zjuCe#FlZ295Ttbzw_kvmOT69ecdVZmi~@?^8=5aUigQX9{crsK!Ek< zYioc1=v_-dVCZy7a40R@5oNgaI1IN1#%_Wz8-ZB>pezDYJQWG`h19!8a?6EirDz+8 zdty$T(=3i#)nNrYpw`P|Vu^y$Ve@&N=2W}6oB%Q~D1dcg7Yi&x9zZSj%4j=E-qQ9RS0y^%MzBL}w(cj>MlW}TwET@cf=+5{+E3qacaMt$Q( z9Ij3-rX}OLppiEl)JJS`3ztrj(7&kRr$4(cTze!KP%EN#_G*`JI`7cym@F17p=B4G z>O#<%ir9ua3x$%y#3zI;!pS6{IM%a#H3(j4c5cLkfk_8X&X=Knr2HEpi^50BHhZ z`r9kH!}HN|+q?EJ6mp&BxJBW2m_-Z%k4Bh8Uw^S98MgV1%9zRh+5dgZt&rR3tVV-LXtVRI5?qhZ9*l6LRJuKgxn9Yy-h1b| z@W#<)p9Rp~^ng({oA-=nj2bbi7?xk!w&UPxxwkW#2n9lRnTA0ImM7astBH&S18o8R z9{^ml8Qaw8Q`;sd%RZN`O(JrescSF1e7)MwUOV3%&>LCGkAHGKNPTnLUwQS&@!93L zOhCt@Ay_mz=W@UQ>N;^3c~?D zI9I5GBLJT>9fV-E;&>1h8-YcE0l(DzaOcv1^Y$9T5b6K@cog(Pa_}ipw>5({HrL;FGj|rCcB2t zN5eIIyj6`?NXQ*+%ys(q+S>m|*IU51b>(Tp`|ZBFJ5vgqZGkO|nVFfHS+ZD`WJ_k| zICkPBj>B=dNtz~Un>0z8n=(${ay#u5rZ~*w@3=`j)7f{=@0aN6TDm&-+;h+U!}C1< z|NHihyx#2WD3#Je#x!Cun0B9PA!hGN>dHy|zKOnrlbuuj-NBGuFJvjGgi>5L)t7fV ztu5F(u=uyaAu2g17t(iSf@U3$NkMfw7%LzC?uSk;48|q0e_b}gU#}d4k(yza6|DZ5wq%rC!hJ-AHQDt+pqumV;;A{qF@ZQ`!ZgoRmYF{O`}D7 zNJ3CD;3^>`Z9qiiP#(R8%5O%J0AvUN9yUnmP(B(Xmk?quEenTf+)zWsBN4EA6bxcg zs8t#so7!@`Z;RS&R!V71E>$NL`u!qE^+rG~X>O`+g*9=xXab}bTek++SdD?!;vqG_ zBcsA0AfUv;!_|7B9rdGgIeGW2UBeoQI8#x5ZQaeAYX7UM=D(ZjZ)s>+*96_rOx##1 z+0XAD*;BDP#Apev*{VWq%^AIR^FN;{*uHQ+NQAb~X(*5=40<&ReyfU4v*~1ZKE*;yKRch?Ju+YrO27Hw{Pp8s zd1Q9Pt8v(9mv`Stc3noEOO@N<5Cbt~hwRv^tx1as_kKf&s*FUkcvTvcHWYaow zmSUUErDj{zypqoq@Qa-WS=1vL3t49eySGksTw3bOhrDR;?80Q3+YN zW;@hp}eaCh)zL?hv)sigR8HPJT@UCPyxv<6IF2d5Ed8861DJ^ zSiS@+<)dv*?UDVn&)sw6iQ7k98WyIZmQRO{_#GqTH{6*Y4z;ygUSHkBrQ^VZ-f7oR z(Dih5i&6q`{z9{a7WG()VSm7^G|5%R8sJ>7MZ;<|{8q1&lX2;?1!KY^ zSMxC|AKZY14dw&`4nxi>2~&kXd+)|h@S~}e&rXGy#3H>Ddw$`&@b%w*yR!1Dl~vmC zjrQ=xJ5)ks(jgu#dLvGOSq7cQT^{nEKZ&p*CreiQ`G(ID!q zqM?EAPSK?ayVN`?hKfb<0M`uBq5+PlM#U2myk^Db(WFI!0vS4Rg=}h45x|{oT$7p( z5~Wf;eKZqrnOZw_tn<^!D_fFRrjwT^5|=lpAJ~yUFzy(1>#Z)scgEX|PVVS*U{TE? z8cIY$%4JxQuvNxI0nDM>!grXs!&&>YcOClpsr@e=?$k>V4I9_hRd1vso3eh%)vb$f zeglmE6Dz-Z^bb$$_bWxf9Fg7TcQ?SB52Tv={k@Z@i;4S-(V)> z)5yrpReW3%uDOP_04Zz{p+tkTKV%Y|M~^+&Tf3@t@rOuX<12| ztS1Z5$dR6yHDHk1Ozce9)b7=EMoqmHPqAoP=n3x_OQd{S5ZUVpTgHdm4o!A-w3+x= zBm)QIW8qw4i&e?88#r>%{(QKD&#UWLUyIi*cnxg zWvr*R6z-fW+_Aa+&YA4l`HnlbcHTLa-aAwPWb(XQadx78PcaE%*Ag07h^t^t)sV(o z4h_j+w{kcwIx)eihX- zvYG*A6#?5or`FN2>#ZbIDID(V^9GH=>4@oI-nOUBA25pjGTvNDePXz{YhNX*mzlNn zR8&y%DkpM=gM+D)-QfohbUb-t`hnBaFJ3)x$Jx7^MI!6d(W|HSORBO?;P?*PS$XragWuE8%tace?gZ4ra0uA(P5y zlO$hSuuF+S^(;t%q&gR998KapO;GZ~h61 z-Av`QQu(b6DT*t`=#BI%i+wL#J@w&(3qL=RkyGF{EhA!2Y#BQE`Dd@becwZm9qj)4 z?#+KY)AiNy{);;@J1WWHgv+SqWrKE?O-ZMs1T;cGC)(Xziuuj3=Jhf;cf69{H`X^h zQE&-)81fniTC<9m{;vJu)t||gNDy?QAgViqmfmQbO+&au^rfU>)U2?m*law4gKbrF z@K!F*W~L4n!`r6XA3rzt(&_OB7u!yaIJ0KSKxg9o^69^>fVF7l>o=~s?fk4otPl> zh%Hh&fk%UoS~n4zYZ)jA0z6m|O(=LR8;$Bq1^Nm>A&p2P!dhGET4B{Z5|q@^NToH( z*%YmiWDrm+N`^s9GK)w)5vvr`@0#oA@SAdSXTlqDd9*%@MoLF}&6Y^OgyX=lycQg{ znIwRcB`q=w;jvxQFFkzZgY&a1mxlvp{4dX+fA;Z}zkW|qylb=j;$rdbGu_`Ft-O1r zu)9a@H!~G-d@f)QcnoSOEpD~MomP{TZ?{Vr0z8%3>a(c#jQ}l|NWmdh3*a$R`TF2r z|KDP6YFgW{VdJd;#ab`~K%&#dBaXNvHZ~T~R7Yv8;d8krhsX(%H6G^5+t-)+PhZ+U zzT9gZbc)9qE4WNd8QZF{z~jB!8RD5_X!MRvAFWge__% zhoNHwQ;mtoH{wyvm0;Fp5!u983y&J>H-3Mx)mJl6D66d!v*>P*+8>f}>5WDjI&Wa9 za3~f7>5+3%?f`f|`z_{D#TcKmW_vGn|=N|g~?^fNndoLXG$)Rc{ zmPtbEEz(5XH!~c#Yck{Uh_h~a$iVFK$>T{7KyZRE?0j2jXNN&XCR9A4gCn^+md4Ui zqm)4u@Mt!vP#`4Qq@1wNrP2sRToRv$;ZwpHKPFov*H}d;nxguD9W2*VL5dc&F#}<^I=?^*pgF zz9Xp%+QlKe$|3_}j`hgtV08^xjX{^%XQ0tBa7fiAyGv^~vgz1bi%Gh7m-UBkfQ{c@ z{^IE2sqgoE>ftgZ6lg??6JZcG6)O#%EKyf3<&QblWrODz-@bqSlTR&YCjCMii*9z& z@7g}SXW`IyT+yFTl_OyZ&MIXDbc~U>cd|D--shVx*@mKuNvM@B)Q;jL6W7LRIahSUJv8;63J3Tk^n2&g@hap^-QY20u0TLl)W4CK-Jv$^k2 z^NWR2Dxc((10o_DLFU3q{1%D`MHgYn9N0`g@YLDEPc9FCd#?Lv#k+j5@B43o*KhGi ztPi;Aal_L`dhcIMogWCC>2-F-z)w)Xy(T6!+PJu=N|3hbOEG=FQ{WT;c-UtgbW1 zw}o8-qKS*M%dusX=F~E6{6uat8pi-}-1V=s-}@UyNT6Y$e-Y^{hfA)d#cbLk0VN>Z zVmgUb&BC>6m13@%u2j%tPS)~x>5fCA-8qd9$Z!}$JUSAMs718aF&b(NRDzV$>~-

zc`!%EJ9pMH9-I3uy7o* ziO^gRs7VSEG^*xyXTXGXq=GKLS?^X0%yO!pMN%qc+j~9N&-6E+-MRA0+i$#j;gO3+ z$buFOx0%9(;3!ZGy#|CL$1|RnPftI6c<{3m#XsHC|Jp*^llxblo##J$$fcB0NMsy| z60s|H6%1SR*2S#tc(-%hCvpiYkPe9?X%fzSSS?XL^2$%xv{!xgHvT#*ftsR0J#N-S0fSu!w9in#DN-F zSgrsg!V@Z;2?t0c;W5#uC+g8n7kc(@-!(s$-dYTL0BarB2>3D}Bf!AG==c^2zD`Jh zsKw2FQOjr|nD;t@7Gu<^Y6Gee9I4U6I-TkN?yP~>!Y@{SziO?4*3wi9LWx2)mVj=A zK&sdT^xjk1uYSLBb-Dc6k;2Nu-LFlTo_cc{eDdnsFP&UoI(y=-i%0MM`6JK0dT#7+ zO24gW-&xQX%p4uBML>j$TA^YhT+3{>lF-GddGD5ioimkld$ymP8lN4g9NSTx><>$o ztpOEbM^LxD%{`RzhOMk1n*zK`4GYR4z!15RUKbV&hDVc~p?RHI0X10uGb1p$Ss%W#TjWQySf>eS9Q^tx~bv6kZI6P8zOVVzf z=nVBFeaVn99kR7YgIzO*|2|I!}6G;y+1zpKZ-lry`OL?OPG*-|GaGys=@ zONEMY2){#|3h2#B-JS#Gqa)k57xVYeFHFva9Wvo?&iCd+(~s;MetLOowqg^Zp%OwJ z15t}?+5~79Ag?s&;b(O$Ed#;ABP2YOn9!={(miIKUB$5~RHc}t;x}>$jSNh^na53g z90`jy=aSEK7KTgS%QGkcrK{lW&drlVI*b5s#A1LP%ck0Ow;*cP3Ye{K9ag}Fa)|X2 zw`M9I;e$kP3j_mgq~Q^>{oBA7$9MG8kq|naD3CG(ZfV@0k2oadgvF<^XjFI(3c+O% zIcS96u1VP>HW|mN*F;>B;bN>i-B$9Mav^hjE>QGa+g#=pM-`7t)U2xN3WQ=jxY!+-7>TcLTSv zUA?2vI~13zaZm!ZhDB`9P@rH~fLbl0r3O`B+tgUY7g2d!kpzfjhW#B$iJr^vPpQ6o z@Ys6~@BP`so5!;zQgaQpWj(y=mO3D#QFn`w4;!_z+tfr09R(DD)Dn_kDQ{0(N>OV$ zZ1a1yZYjmDma>So_{Lh7M!9ckX4jU%t$oqs^V?=7+fHvDzTVT1-#zoUYbREoyZ4(% z@6UOhcx)9wi_wTJ2uM}!#x+bLOv*%PI5?YuU=z?1VfAoF(5#{vcw{k`C}b02LCZ{! zPs6~fxne-BY;LNaNqVwAeK`{! z%tyPkwsc(EmGtK$0gp{|asS}TZ+`a6_YWLBy!(NN_kfVe)9*j};j5S4etFMp=jKn& zq%$r~Z(9A>!EN`QI(B+_WY(DU7*d_w3F<@5A4`ru{f$;ZV~#%w~X11+sEH0y<%o_hkdcm|{4r?)91mEv!3wypL`!Y%3?lJUG1> z3L@SY5V;E^H?_-lUzH}ziq5PVj@oZlqz-; zvaSkJy}r4A1ForoO+$*vFo#l@4$HzCW}lb)z--TJ5A1pT_D(aG!f1u?@bv*c`aqiR zrwJZ?W8car&;IqPr@y<_e>rzyW+Wt3(Ps-A6mV(%4rxTf@v6x= zug&f<0Ox&k*Qin{G)h4X2W#a~4R(%{PZrRS1~$RY#jAyMk4BTv7aU4%)TMSC`BFAf z$wa2Kn#un1sheB2vXb8K6#74tFcsZB!$HE4)oJGbU&@CIu*k&%1XjZWD zK3=zxch9blyXJgWCQ?p@=XI1rA-&fMTlskPE(6GS@b3?wKD#_$**ch694H+a>+VW< z%|fO_FOJJ}F_UMotUo=y^XI<>fsX6#&wyDeB*Hm_CclZ;mC7Ahn%L4G&D+&kvm#*= zb~)^uM?-@S_MwRSN{{LGPScKrYSJz2GKqR4I){eP6XZWU*OoC!`V;(pJ>mJzHm6R& zYHl!yhzizu*XP@2Lo-lOfV{|$i z+EmS=qtXFj`6%2bc~5V8X25x9JU2BE;(!PnkFB>EgbF$ji{;X>LAS!D15!4`nStVL ze{V<1)7fqTaJa6HP};4Gdwd2x`{k$Kdi<#u0Zbo-y@4@%=h1y{9Z8NExos-^Vlk=` zlF*PfphSW}*xK4qSBkq!Z4Q&1s-<9=$X2%*Jh`+fkE^Ra*y$EU73}UP_mS=0y=K0R z-O{6B%m*cI3H0Nu*YOM=u8!$fuN|8kDCNwBsB5w>(V6y1#Y|v|C~3mVqU+L3>EVOt zSAO%xjp6^z&lYC`Ty&E{h)d)Y_Z_LCluVy5uoJwfsGi_&=vv&h+y}9Db z>F&Rr>HE#?o!4d~b7}QxNYWFtWG#ZFqPE}4EjZ*Iaht=zbF0vM`?~Jf+tu4<)JYft zEhA+V45r+2J}7*JkjP9jtzDt96yy2luRU>kb`Bu46LxXJ?i2xmu=Muqo=Z1QdK_G0_owg2j zgtGz5V28KQn|}1U=dZo;B2&z%gQ2dU@>ee|KeIhPXrfC94Mo4%X8@$b##JCC6^4Uv zh1m;@MZ12NNz3$ReR2GA__sVQRtS6#RNy_;2pqU+3bILKJ zMFV_37>6yEsX?VJkJ0LQape_Z|KtwbENeVF) z$s@t|h$adFR*rkS@`-^$c%joZ>X&!fIQ?dRQcNs6q(z4yDJC6g(VbUme|R!L|*P_aC6d~VxJvO~+M zIF%>II;I8_C+GU7CsHXjSH)pQ?6Ry=7fow#-fq7)9z7c~oAfF$LvbuD0^16~;Tkc>4L&|0=@aGxs`2*Z>~Lf(7w}41 z0FExA;iMcK8;>F&Yee`4GC+hjY;b91gQaM9yS-R&mHm#Me55n(&xc(--NEie+qtuQ zOQm#EP4kT*TKU=Y$GR62YM^xxH3(8pxkpJA;#<|+CKe8>RM8az76Suex3p+woPqY# z*kJp7DSX(=+wBqbIXRo%!kmSXvxsJa45pqjV&{}yg69rw|L)*`XS%wuJpb}J3K_d@ z(@(3`{0FZEmeKMuDnd{|p6zs%3!24o|MB76#rc6}?^)bG7aI%N{YIvc3KLOa9UjTj zFvxKAUfel)Za#W^I30D#P_U-Px(y&Py>a7BJ}qey5MT@-j7ej#D4=pm&Y=S4@#0MM z#7L%?b!20PoJ$Fuh#e)cbVpbqF;lyaPk|;B{BoOF&L!Yln%Q)=fUlMbEWH!G=a)}d zv_g-DtL2k4bfA+16EksaP74VKbxMf_Hfb>J+A^B%&UyPn3a^^(R}!5ff{@a}KvlC5 zjmU;g&2_gr^^Eaia-08kOuLi{ zt*%wzkTM?0p=FI0<1^)SLcx*Jp<+t2lmU0B7(p$C-&zl?TZe&cjG35y#YCk&4v+w) zgtI;FtrWdUyR;Z{`2Zq>fNrW;w|4CYFfyOqm;GR=^t-EDU)Y&CH|#z>JACf+iA!f^ zUOzd1eoI6v!kf$-pIx5}0mp@{<`FA#cXuJs>EZ8pvlrc*?GEk1f$Z&-+?nZ##d(0Z z@QuV=Acp2tvhrbNhe^8fy|?CDy{mF+@Q|7a zKwlyXBJPuQ_~m9UuFXu`5>p-Tbzj*MePW?-bkMiE9Eo^XkVcr^NO|YKU0>GAVReY7xWb`SSW*RHrhDVePF=)oZjAM9IVz5n7;)>xr#ZDHUlX;~hx`fuds}=FY`!z+6eV?1Mew zyvJ{H2*}ONH?6tp!tT<}oODk`cX^@h>D`6rcJ@4e&;6$lj=y+v*ISq8fBwjU>6EeT zc6EjwJ(*A_U#jWq~$}@6&zlPXr=Z=K!3k%7W-#z!;B>(zow|_EjTCL6?qXsOj zh*I+I6Zc)8|2@Ozp132EGPv~&I~y~gWhF#}b`@?CMzwx>q-9UHbf4JUSu@HtExp@xP9 zu<%B^m|)jPe0ILqp>_rdF%|cjU3Yx={-*ZG!LV(hJtpSM zZ&~{vKv_kGhZv|W38`RzzyHY{m1}qJnVpWD+|m8;;rVAS4_>=_-_rKpLR6I2OFRmu zN{Wqmx$~p$rH;sO(Y?Lhc4V;o+GEFog|a%xzrA#!;MpIxOhjycJxj$zgK9^eZT9Qq zKQv$u%e66^(j{fI>1f9$3ga1Nhk>&zZMuJHbTAEWYpj&mNP}&Z(hvb9XR)9>TydNq z4?Ml2?e)EdHxCr<8uIPwcLy{a6m+A3%Wo^@+e@KxHoUKEipGHP35_ZMWD;Y-GL}KW z!0s9yz}ZEhgMeJS$XV@zV-nYp`{~dpBasN^>jOr*q!zH>`ECOJ=)c? z=lIgTy|)K|V~J}7g7829M=cOYSWLTu>okeXGS13(N=`=)j0JsCK=8yjLqTd-DW~%& ztw6cOW05EnED0S;MAWfK^<5t2a8EGP2>_X(E{ru73>HI<&Uh#imv-Ai8XIrp&9_eE zJfm6fNFJmNERz}2{=V47BV*H*$XsXS_Jy9^gO=U#!~{?{QgTF$IyGXG2a8KExx0F^ z=a=W7IKA-Y`%kR=;hEn*cIL_X-Cw`{U$7}YwaV5p~mX0~(dVDD&O4@+1pq_hT=)RckM zu3-l?Op%DOV{~}st=Hds=*{=O{Mq$6aP^tfWat_`3$eSWcY1WL2}mHVoU&bm@oq2W!rpy9pv$hG!tEyp9-J(_fB(hR zVj=&d{q?0s%3l9iuX}l}45a&wbTk)Ew!g({7UnD86qpnq&)uK-WQ~ptX(` zN)KG^Rto738wgJ@Hmtd&q3X7#s+&Q%2UJ590=lVgVm55R)Ih4> z)gQkMCdoe<4Fz{;C=P`{$f1~3)KA|23S^=j2K98sVc>BDWQai`SE#L8wZ$mn_Bb_h zx0Q^pRWS)&ad*V!B!kyK6P3OZt~15@MSQHT569zWiHWW1BDJ_L`F^*K~Ib(*oft#FH0qpl;s2>kx;7l5>dwa*0DhBrNpSa(TSZs3fBY z3to$cLPIvm*nEXTCle}_d`V2J&X!FKd<&nA0#Rs_N=(N>l~l4>L^X@pP8Yu~<+hsG zL~Ij~i4d|eECK=#6fy~L5U6iyu13H&(vXd6A=RRjbn0~pmnz~{r-Rl^K(C@PMKmG} zg}%A^C&i4x0i<)nM{H)=r9Vpxus>z3kMdz|H7v|gVys4U4Q)Q;ZqaY!&{f{yt44% z&dHT0US0Y5Ie;7aNmn&W(S52e<@%&u(j8QGh zhBFlMZz-DxbJ~Pcn+ciAX-9j)XRymXCIg79Qjt((b3LqSqg^hXD~9GWjw`#m?>)F> z+hDX!Eyz0w-;dt^<8`T2DYpAYd*hvHS0!QZP8-sGy^w$gX(qr90jf<{+}f3IR^d7~ z+Gp>)aQlfrtsGf-=$-F(U4Q-6kxKHOEp69s-}}zjo*H zuYd9U%EzB4fkNsY9ok7 zoB5*df@rv%iEP-2ZrB(#se#CbN`n5wxtAY3Gq$@NGf0^{d?U88kr%^s~E`5WJ|rds8`VLb#*4}pqQdxEh>Pk4kBVLn;L3CwXF35 zX7f-`xh<|94oHSO0(Tr9X-}9~bol8##Z}beD);bT+gLUX;v2(-&`^79u#g@sh4V3+ zSb$^TA#5-!Fc4NN|y`^7Um{F7}dkAaj> zU^lmE4)z$gWu#k^ioAu|X;uKzw^~Hqxuav{L%?ucZO?oP@K^t%z4HE(FW;c5fxk`k zM4mkYF4?dC{@U-yhx)ZH*1pB=ayjpa^3+Ur#19JN0Ku+WI+YB{cgY2EHiIRRA`@=0 zmV>2}3Ehdn(sb$keG}jP%qQnF-LZ@e1d#~18X*bC7vMo5NHeo8u9BpZE*XV{YHC0u z)?0;aAlj00MPrFFD0-Ljc+xRr)S(X9)Y~@~6G3-ePiJDmGiwu%+LY8W`n}?~ik@}~ z&+Xf~Yb+vB4O1kWrGQ2z3XGo`G)y-l6+f#O7-j<(=X-dzs* zf>uq?rH!~Pfw0yn7OnjHU#9wM_uu`+gWsLS5C8pn_hj(FJwq=&a_r%gQ=dI^&zo2G z4yM8~W|I)vWEG%uMs6O6b?fPIDXBZ4FL-5ft?cFd9s}gxmB0Ob<@ev7x!C>HgJ*s? z-1Ut=JO$z&FaEIe^zi6kR^Gh+_Ugr1HJ`4e0xX=m&F$1l>3$ZzH0WT!ih{E5R0l}9<`k>>@2u09NKzlTmRl(Z(2hUVjyA?#KXth zxhMq_?p3kpM%$Mc2Mciro&e;~t2i|1Y>(paFJAfb-Fu$8yzA#LJ$CiTmPAlq0*d!G zmD4J)+jymTAny0IyDi;O8>rARK3t6W^v9Rm0V@3azI_9PpjK8HLqAfc`u2fXLbi}MSkwY7 z@V|_F4Ueu7l1(bM90*=iuWbYceKy`&Q*~=Y_052-2HUiucHO#WSiMosciXrJ`ffl$ zzI<};(V32tPZM;?d~(6%3p1;y`T6sI`0B+cp4k`ju)Ca6pNcEu!h4c|LdMY+GUvmF zeAt=@8trDK->NvYd;Ed3^LNaac9iTvIaNk#)G%Q{nFnM)fo_aT!)p(^+k$3=gwz7x zRJ-X`5~^Wn*l~P$WOkzS;_1ahi!**F;J{kIGU3sSt!hy`?2NSm-O^MfnD_=mj;F3oTLazIRAi-3}O3$#zI{^OO?e|_ue zt9NY;D5aTz1f-??N={Ew1%x|dK@e{g2)T4n7YPq*LLqAjNEDroA<>~y0Zqpzh0F@K zQE>bA%5VPf{ik+ywa1;O_Lbj%X7Piww;x;T7|B~Zvi3|w9@NT;ZAOrWU@=fsTq6jr z0}NX;;4a>B)7slWROY5NH{T4zdDgADc>}0mx@P^Rb+>`KKd^>+1{QsAoA=`G^o6}6 zyT&S|yhWn`CKAN25o9wumky_+z@=EYQb3Mrd8MGrB>}Pq-g3^D4mmo2^A)gV1CE^6 z$!B2Oqt+c`vE{MU;ko!^&goaOof3jw!7@l_bV4(x6%YWcNQgR(n8M(-LNOca8aJSu ztAIj+Mn;_KaSQZF$?|a4~KU z#Jx69J3y;r#si9bw%qmo9dZ5jola|Rq^q*DZ|B!54}e`?zkTZG4_?Uvo5~KTWI*2l zqvp4=fea27LT-g)5DoPk*3_=M6;f5zvZ=PYwhmDbN?bM2Fi^k&Ct@LXC2hx8L8-)| zn>L|Q^(0CIQ10q&cLA}csM`oA11cUf5irQ)BnAwQfxsJU*Hy2-ZQaecf|@P2-SXc* z1+{a2aucXieH*Czv*wodYu0T5#VaK~&s1mAWOq*Vm47t)_OV;lm1P>j*jznyh3ndzZicUjQ2(1dS1W4G)=zPDK z`{VKd*LUYml;6Mh)Ug9$pO|1}BUw~18P;Q(tB{C|fELkMeJd6MF>1J8i`1u=`qXl# z0cd@Iu!c0C6)}hq0;W0w%00@+5XeR$j$#xk^*T17gQ&T5C7TR5bY>+Ntdk)ZXFQ)@ z>>G$XEd(5V?Jb~g+D*3r5<0j9xBST5P^yamH*D%zudvn*^=I2CjTTiV;@WVIkh-LlSvaUky=N=MRu z8^|?u#&yY%y_EN+EHvoJKwIf;M?OV$)UOxT80|)O~9`DNs`D_e>&;V*B zgZkH4I1s?uq*5?VR+d7-cDtw#pWHt+UzwQjkXBv3DjE}ohC-38HCz@(tEAa9>P*I= zk+R^>>Sj0$i2#L?VYPMD=;j8!h|Xpb0x|l(j`)w?AMMY6@!`4W?%rA{_zXfUvAGJ~ zv<6&y0=kBZgX5qzJbH6yNY@{;RU+n2pP}HEcLW^8ur}|}83BYES0`pu%qBJsLmus~ z+;w;vxM{NQhTR~Qwk2)TJ?#_YT_auogPZ5C+;?noe;3Hz+;S@@X?4?@wYT8w8#DqC z22z>SoS?*z33UXEF`#NY=1HGg0>ILPr8v09!zqWdsje=l6@2~hsyTP?K>L}g?wLTj zZ!Wm<$p-<8I;3TkqMA;>r<@I?!`9A{C7(7Vz2?5OGwxESqpG}TRiS9O7yz1Rmk-aq z{?vt+pE&)@m7|9idZ%*kKCjXYGPW(XSoj7Y>_Z|nO1OB9h(FR1K0G(Nb*%09mJz*y z*^1kQ!8ZcE)CTApXk!%$Rl{Qtqd^S?g#u-{njqD+_3K-pwas<44I9?b$PHv7oX2T- z{oL*!5BZPYoxOYU55GNoESk!=1$-)+2sG|mYk?M@kOr5~T50GeDXlf_(Ex_8oC%k) z5n>X;sN@7SbWl9PrDx*NRWwX9KqGLd=*5|dy*tLrS-*@)^V%d)lfBd(a`}u-v%EX) zF30u#9(8ZXqXy}fjg2%iWx6u<-Y-7*QEz_y=H7cQO!lTP?H_OR81QIlIwQSvat~m% zz5C#S&HZg&kNDXqhE_gZ@dij<2{CU|#J$E!Hdu*RJEAdwi$ zRwdnSk^l${s<{!&d?daOl)D9rVGJG~OKhSt8lf%K5bPQRss@SLghp1k66+wXYpBcy z3J!w7Z*=IXZ$5tg5&(nthX(AH7gnA*{If0NE3ZE~8ZvprKxPU;LN?Hl^-=}`)SECe z5dxqx%|PeeQaOVN$YU)vH^bN8W|UG&CPQ2=_UqVW#3mxV9u40>Be!-G5R5Hf#0r%`c zEFbk1lkP&q8uq&q5$}A5@6;5ispY@-jx--b){3#4E8&r>rOPLJ;n}ONyt%J0qxJxt zB&k`!!E1y>1s5AIaGfe9;8#Y>;`uHIP#sJt89^RiLTUsex>P^|BsNO$5IMD(0hHXD zHZ|9;M?yBCTk5$qG=qc!VqroK#%5NA{X!b6mBxojWN4tWzLaubIW+slhv)B|A3C)q zbNlS|tUWlP&mK6kaA>BKvhsk$HKLw3QnhmL;y#viVaWLHPnt<|vcEt@u0fwHjh zrY4<`<+SrZdFx*s!+#$58{5;{GjbZcuiLS`*Du5q4I<8?7tVvfz3|%5*^+rNZri^P zMEt&abL47A(%x6L6nctd6J>yv08&Vms5KUK=2QOBw!q?KVY)wFDgcfwRxTpM9D0w1 z_0c^b+JC+MDd{ok83>1fU;zd4C@?b*rDNiN$hL`#FFMo#HA63^SQRv*ilS9wq*N40 z-Dd5)0k z%AY@a{Cvk9CF|MY=>FZ=&)<6H!#6$~+iB_#2J(^i(PF+3c8{iPJrQdzYVIufO27&N zY+FzekVg^V5LP4CEZ_%i?0@|3-~I&fr~w(H)4=I63k!N)+0IMp*cu*ACcx;VBsCWy zBcU}cyh%s~sa`r3(2l4;1!Z@|9rg>N8nIbOGqX@`AyLgmT7?*|lAz@x)pVp)jA!BN zV2x|Yt)Mzy3zLZC0|pkCESJ(GAiUr=cer&L5#49ugft3^o*&S00YpT?LILYjv81C-#*>!gw=j!$>aScZF#YTF4PqA17Tdc8=vo8XRK}yqc{P18 z%W%Z%)sRXiVMV2M(I^%++ATuacxVl!#l*s@nXLj6(k#SRoP4IH73lr|Jt7LKNy;K5 z%@QY{Rxoe}vxa$4&^B)h=@?9Et&9m#(3-fI1_l}`z(bToI1}3_5|QLG3Leu$ryzJF z1OUn!#CQ$?3QE)oNvNQKw>VsiI?O)|$+J)IUnrUP4@5egpjaG^(TY~9ctNXG!Nr05 z#3CSxK@A`_DdH3b6>1Ra*;fI$`pS->(!%Dlm>S)b21~MjQ1SJFj!3E(8pq7b8GzdY_kxjReTB+CllxP#Oxx_{Vs0=713J6dUu^D8R zsg2~{{{BZX>Z6ysYFFW8)&G;AhR6cy? z+G{Ud`SkrSo_X()LkszYM#JJXWkW%!5CPa~S_Xy$uZGmD1+{ave1S}jF9Gmdd#Wqr z?HU>EtMr55b7#cale2b&>=WhqmZ9|jU)Om+w{>0V`mIT1OQIz9Uc_ENkRU;#kmvx> zdk4`g(FqU)NP+~|dle~>l&DIH>dor18B7B5_;~omJ8RcBCaxTdUpN*!IN1)gt4X*1NX*jj5b7zw`a^;$1CM+Z zOiw|V$tW0)hNl&f;T~$>v#I#<6`$lg!}#1M(Azy$SEo$OA)spO*mRPBjHcm=TI#ES zMT0}+Hls^9O_ep3*)U7|&M$8M{Rv7M$6Y^fN%f*r}5#Sz}X^hqpJce zO()kerXHN=efPtAuit;;=f8RXlVAV%%-r6n+ce_mKf8Q`rYK^RgBhfv0^$;KlghG_ z1kF^j1YuFNMuYAMm{xpU{XTERsR~*&y>7j~RX!dv>>248iiF1c)QfRvhg~jcX#i%4 zUPiMUg+@8WE2eg+S$Y*P%=re5!f2H4i@LuZA3QnGYwp(Lm_RHl$N^LZ#MpJKv7|@1 zx-1;A+#p0tu%%7tyeed7aZ&Q_vgArsW(_K{t&MZ`Nc6;{?@&~CVx;Z(jC*|%U})_J z#@r4a9~LGq89_(Guy6=29wBMQ$SJr%mwGN9Tp8;g2iv%vxS;I&?lAuC?YF>`{lg=N zHEfaJP-Axu*hDT8&)y^oFyvia$Dpor((~F+ESyDDSf{hYr5^5c4#mw-Om<3WY6ji~t!FW_*C2_Tp=HlAsQ5OM(&x}$oq{&g ze{F*uxOO2PS(#apn_7~a1m>5RNhAcS5h6uILPd-%Z9E*OsTk@8HRV~jh5`wzl8Px+ z3o!e7ZO8hxdxOgP4#iBny5A}7_b9r80xys;`B0M)jRZtIDheu-3~afKi=K~YE-m*Q zn~xk?>^U_V9&pC~YY^XgMyKH8z1CRR%EL7S0KbKVY((ZkoG3J#CHN|22|UlNFUz7; zXSR{>h8DEGv4RGdyD9{ShThX1|IRCI{_W2ndS;{|2aTd9Fs+R;AJsUdbdPC{&vi~M zv~T|DQ?R7(KQuDiDLoW+{N(AL*Y^6L7ENVVTg2j+Q`Kh>ce|u6vCyt<>GwKg4sDOy z)D<>$h2%XpwOyms>q%E;PaazeciR*usZ1#+Nttw~Q5o{8y&k#SCFKjrA~D);mc)JP zpkD5@Xv{|5!$bcyg@e&94U&5-a(WZ89*pGmB?2+lq2lrQ<-I=df&KCBu$D`%k+f9e z>jA7?2$>F#f`%KiRg{MDfc8kAalXqm)M;Ign!&kgB%wzYl75HYDQoW6aONE%JDZ3{ z7O}`RG8xTmV;!3go!uLm?{V$x3ocHyefL;>`&_E2)u>r~Q9}!=O2uku!B^K+=0X#S zh{$AB?=H+ugoL%~;wl`r!XwsAt+`%!?FIvlKo=E4u2N4MyymxkiqC)d$A?{u8m^7O z(BU`^qO_w$J1S`(X&vj*500vS{KCTruissq4UStxXF6KnS@TT=c%;&DgGRVIH$DWd z8YNZ3!An^PljyLh)ou26S#ydpkC1Hp|x7Z#3ED*rd3C>E9vmD+H{Q55#Q^_$G&Z~ ze>19E^Amj?zDZxacd##D6SXu|La<{KvVc?x$&z{H1t|sjnT1)ojp$Ng9eeZVufT8l zC)+{Q(54gzJnC=8@bydm{13mo(>-cqJ7{bRNfux!dU@J0nLB3eAJC3=>Mx%@b^oQS z^WpHgo&Td5^T$WMSNnWwF1Uz_V6jaZEWCKidH(`u$h9MeOgp@@07Q`ha$gxb`8 zyQUj*%IRb#2@y0&LID*KQ@gt$r@9VMsKrf??L}ztcgXwu%zi_QNy-(l8X^#AY`G53-oo1xW{b#_+RgRJ^pVV{^S0&^LAg z0k3v(WYDVP#ohX;pnfQ9iJAEj*LeHvH%-;euYUdmxbg&bx-P#tpx27n&E+LI>1o^3 z)3*b)cy~4k%!>6ghMI>5T%MAH9&9(g^I&!J7x1Qk`PFNO7oy#RQ|6V?$XwLtQv-CN zh>R{|Vc?@j!Ys7%(VVcjE*mvJ?jF`gA(9A2`SFXlvFH+~Tr@M{|IWx^BKnTMzm4l` zZZS4Vd^C9%TN{;_BJ%NmMZC@O;!8I_`{Lu(*_oi3(5s=1+gZIv4j7Bv4rO21W3dV` zfD$jt!`J18-HN%P_K3eD9=G&(H6gn`Y*GizdM3MBt743I*`2K|=!(+f+zb@599dae zT$GJ2tvoakzHxZIH{ROmQ>z68oqz}$#id^Ra+mR?Q}gfMTzmEW-qS07yR;c}k7QhR zGd6$R@4I)ix=&*Upc3IAiVlkq&2v zPW0iAFT*ds{e#0xy_RmDZX)D~d7N}=LuO`DGSuV~wj?9~wH0Edx7Jk^ds=yIA|~{m z=*^X|Y%xlj6cQAfP(rD$qJXdjSInlBy5)Qdp`se%dMkEABpee{-dI~AVKmsZ#`AYZ z`<(j4UenW09*B2ySQsRz5Z6Zj+dbv)p4nEei^yngWIJi%09zfBm@C?L?m>Jkd?q_?eZ9U+AI_8_K4Hi zVQpWv&sz;Ygh&%+UsH{92U02eAuVgeL75oMbr&$82 zb6kSUz@6x}t_?d@hO85C84OvqQYr_JP|<4z*t)RC^umS34`00T?eTy4#e*{%+H9LH zV%Gut$fjZ#l`T?A9jUrVFKM!hNbLpzm}mkPDUDdacf!8;`A2hcQJaXS=QqY3^0-Ua zgsXt|->w~7wmtU!?T>HW^{6ei1p=pJw7RHEuH+LMA$kc@Dxsk0*y^I}Tx3-~B7ZkF zCx=^G8-UIe3ty0yTAKM78j(hCEW;oQ74*76mo9E(v=S^a@ zHDAcYHdJNdkX0gibDtLUiA|!GMyEy-=+gFf*~3;tr(5cCsH}Fa(MAi~&8PPFXk>V= zN@|m`?PBJbGhos%nyYFtwe{7AGGuv4V?(i4KnI#`d42(uyi8K|bVu9DRM@SS$XlBH zt?ZafI}&r~L?nw0zc3Pmp6?&te_``af89IN7Lu{|b(sJ5-YXBDK0exGgJ3c?1IHj# zHzRi2Rn0OcMJ;d|S* z0UP530n(%maVYW*~nTTnXOgsl$$!^tJ(~*J)<4t z{Z3$4Lh1*(9?8PhjCEP}PxV-IbS;lH9rA7d^#fS-Ja=-*q2Sd6ex$aTiZ6y7Unmv` zx#)mUJZ20oO!wb9a^&Xfxf2JDUffuJ;rzwR2lu-z!UiNFCo?G}_qWgB>g#qTZB5y=Jtb);WM~$Y zXV8hH3%3rvaO3&g_pf~L>Dw!-BOqFv*t=luBi0sVfDwk+SXN&Kl!I(+QIeM{?g){) z)arqVrKe37by&MX0MjsrjICoa%b?DBWn=J}lT%0c&HeD=<@ay@@b!C7Upshcs@)4A zVqurgs^oD9L_rG?gfTXa%qZkCns7uc(xa5k40)EvJw5|RNUinRBywppz}8S@*=Qh{ z*H((iO}+M3msY6cQd@=OV+*iW`t9exda*C0U|}jLH3d>qokPLWD#;x#`TIY5Z}U%| zZ~o;U4wXTp~4Tfe_!>z0JA-`lbMU($CZ)EDRX zc!iht&rQO#thp2xQ26SKfJxHnpnh`8V>ohMYLB6f@-%#BHLbe(OID{+A#8?9Hw**Y4b+Ea>j!=BMT4r2@1%FQc-&Zf$Mm z?zvmHp5J)+xeG(xZVaNxs$?pu^?B*(6_r`2(wwq9=oRkB$x5uP&Jw_xCX@6<%u#FW zTu3wGQpc>|NV3ic9CoGP?Ec8!p@$MMTTMV%xa96O0rZ2$+S@lj{2$7UU;XWqfBWFGLu+FKpf=}b z=I89%wQK96wm5tR9{-meTmL2Tu`Q7GqoARW4z;h0v{4(1(h{~2Q56A~s8jl=qY7Kl zB=Ty7UL)^=+v_lEr8bo^iD*i7Nz6#S7}eb!H|?|Yyn^O^9{Pxe<~PzuWBxyV_~sX% zzdAhPk_#dItu{Lk2qL+=bMniI%hBbv)yUeSf+9YXJT?$N@#KM9>uXWB0bn*PU=Zi$ zA+wIx+XH_&eB^;q{z`A4Gq-VWE)M07ONVYbs(p$9$C zSCpAmnqPz}FF!s#`OU~LuOGd7VC2$rA55GE{5FX4WH#6H@R)eJ|5xw*?5P)TJB-Yz zPBGHyo$Pf^M%%o4fsBWfvZ+Cx$}Z#M&}EQGreKpTtvn5fIubO4R_@M|$L5Ayb|uTs z!&(*iKBpw8Z|!k8H~;;EZ{GBu9{>LNJ3oE<&~mt`8c~pyk(Zy5p1eKr@hu6DZUNul z`Pjeg*z*1ClmxJ8^ypXzraf@=D=$bRpvwJPF1xXDVcq-I56;cRl_9MG2m?!_eg(G~ zhbm+?W3UwkHY)OLm-td#dpM+8X&0>rMP@dNj4BjT8eX`wcT0N z2PTIcyev8+%5JExQA$Y7RW0ZWg|M#k0~jU2O9F*9Iajw=_!(OXYVjsXUzO|v-MzHht> zR8^CK(7D^^CI>?z2BjWZh603dQ#G56KXxRx`Ipb0IyJvC+#R%t^pbiusn{lAdiA1+ zU$*(bzq#f8@t3bW_uY@+_kaCx&*ETJX;FR_{79}s-;l6d@=*GT% zPu)1a`IjGUe)Ko!cIx%A?*7(ZmwB?=b$n^~^2*Z2(M5wstXHs=d={ur5k)!p>T;)? zu{sk!G#{Vu?eb|!Dp75HRUTB8-By8?Px{SEkBaT@UBC0UPk#CcAoP2KGzS1S* znN^(0v2IZ4tYk`#Ur$G;U&|y{D}_{tnFwVopM+?W zkUcuOj)zqe34lYFH6}F3oGmFR0DSuPM&0Db0m@dcfSeu^fB;@%}Xd_qvmv zYtmX3VlIhJtfg0>nyRWWwZ&|Dwcjk-{L4Spf8ev29`nK-w$s9y9S=+aO;1d$Ey?5H zFkmlkKvd|(90tCbL9Pde>R8AcP*D^ubq!^iC}chsUsY3)MQ^AH`;=EV`mY_jT8FEy zEJGI+SMeCE{qYFY_4GV4grLEBtLI_BRs-Uqq@CM$ZQqs%KZ!dMlXr>;wHNn9espH} z>Y<@nk4?kF8l+7DKQHQ0&I~!8I1n4_P)BXTRuEOv>VmDTg-)}WUX3l!g#2`2V==d( z4AP35>hoa-xvV-dy^f2oq}G%Y>x+5xnobWCDv4xFp_YmO(L0Y^&%#wI=|rm#+od23 zxRsE`9W+aId}@DC(Wd3gIaDqS%_G+_2sNm(EC82=+n-&pL`a*6FKZono8dV~=lFKm}Z1GX-YB^ok^dyQ_p(4=E-)IR7ldGrFSR)#FxjYgJ0&NdcP8Pa*q8-ul>eTe4jI@-@%=T-*AYzn%rxW2Zm zwlYh>!A3icZaq!K#fxc8CJ`m>*T;ftfF^;}l!L2j0V5Ks0#gG9(wf@JLNc+0OKL#X z76WZYD*{-t$YbTV8`vFokZ_}9=H^6u)x=?6);gUHbjTlsaQBZ9i`z@jADvW$<#}58eX$WLounj2Cc;D)98fQ z-d4)1$6;#y=fD2u9*1!^9Jz=x(Cs+1C%hPI zn+{t?hP_%n8#KU487ax>Nr@>b$%#qX8A%!mbA5L3?)91bSC&i~5e-#kCenOnzJb#O zr#*nG_~b?ir)a^|TTMchh{Pjfs7>V*EHszP!J&pMD`-X*HJ8GrvXTty#&kdeB(R$r z)L=b9Wj9q9LTAith5e{#sH(uBtFahlLv5LgLD)Cm5%U@4YyyvmW1<1llbEt&OH$%i zs8H?N_E^%kEk(J>Mg?hZ$h$ZiJiRx5W})}QLg>g$c;AHY;8>_Hr1O}0J#MRzjMd4= zGE%*Y*{A|}D!Hy6osKF>Y%a=2S7qaB3&=GEl!i(c$ewZaokrbD7w_J_e0zR$7-q>b zQnikWKqHIlkR?sEpwj1#hCF*?zGJhIh0&-b?z3yVqa(*v zD(aws-6v{w*#r(1Da53U!4+EpmGjKPoRsYJ?M3<9Db1zb?aGVmQ+G}rSQrYMcmyyR zmE^X7U;Msgc0G@xYPM_K}dK)35Ke>sl2o78Q-JDQRjf#n8%gb| zq?K@=8()g?DC6IEHE;L69wherE*S$aYlj31r}E=Sw_5U511pB_>W(8!dOgxrh-M)jkhIvTN#N5u-8$Qs~BfUlse&g@qO zY+zC0EzfuEUmZC*-M+B|9R9Z6E_2AmhvKm{C>rQ6ug`X$p6p+r`_Gl*KL-Ek&uh2$ z)l?T}7bIt;?kLVmE6vF+%-V%TV<@lZVbfxe<;+`dwpHS#D)9BHq~vS>&qh3Mi@OSu@x|*8psmv)%x4owdIde8<4^m=FCUxv%{$jW zd2sW=*^~Q6^-dmUsmFlBRpzAa*jyv{g=6g5ByVjRtLw#+ru*D6kj`~XQ^iB1-mdA(pj&we}HvZis|IPEi{{BgR z6I&}0u<6YLFkn(>945&mlLlMmK_$&ArFCf){Vt_PC*zQsNwt-T(u|U-3`ktF*@OpX zIuDHdT{=3sxtK+(uB|STvMB~54FMDS;ATlWsNJCt~IHKD$efUe@xAO{Ts9Uf7_;4DG_04QlEEU33_d3?wA zooQKV6bc>JP@kE-ZGC9|!GjG7F0V4L*eap7Noair#w+KSfBp8;&)+-s)TN0V{Ug&} z(YfWQLWD2go!+Ju{PM~(=dbs71tiC2JC4tIWU_j2#PkhlKE1Pkex-LRV!U{A{+a9J z8{<94_rxc{zTTkTAm>(9VR2q4*izjb9eDm6yul?&A@7yjeC{Eaxknq^{9Xk?} z(o)NcvwGc<$^MT1Hu;&p@X^I6y9HgEUyzlT7u1UW`E2>=+426nXO`|PufBX|{e`Dj z&#Z;PLPa9fkO_6fdIXP%Ce~G!7G{-!BsD(^`Z8(BTem(2$N#bI2?;xsq0F6_n*|=2 zofRdigoc6)X!|66kBBL7sW5I8dVkb%|K!xY<9iG>t(wr*! zRJ6-#sSwxIKtqHpW|SADmgc5;l&ZrEzLR?c13}N3wb3i5_v{-8ZY;${eC~)xsZ;Q3 z5LrJzcXxet{?wk)>xU=LkF5Oj9)0t--)cCt#fiS-dsj}a%#2UPkF1Tnar4yc&s})) zrE7kzoY0tytI8$R7J|3gtQ3ID-ffXn8TIv;-8L0rPhZ>3!-H$FHZ={6X((f}QT<-t ziSfYUeSw!R-hbuEZw!@>-ss0yZvW=@_lk0=5|iK(1m|P|%y6=DQfOoZk4y{q>Cf%! zyD&Cx(&e!toiOnwNU{NAMqnt^Va2ar8 z0)8DJs%Z74O?4Icy8QgCt&iF|Tej?a{PA50+mn-b}9Zo=6cG50@ z6)UMtdt;i*bMbz!aecafac0kCpL2boFCOys`1MY;h=?l&5yhdI-djfo@9&Kt7-;{_ z9{uOj-@N`})U7)*H@ALretEiQseADB>iiF`9lm{j@9`sj!@gFRPD3LVqe}~k)%h}J zqtnD#nHWE}zBjB>3NTeHY-PVsdH?F%*?IqR$n?s^i7$S3e)E$X|NfKP-+ix-KYV#? z$rfqT?9NTvxobO=7gBcYD%@R)MAg8|f`)EjP?~fmi9R4`!b86)B`-ULQBU}$fxEf+ z!REg`I6XSMJfuIi(Cx7*d}gIpA+?GnzBZ59Em6zp9-FGo05~_8UE}z;Dh8nxr0jXQ zNl-`1%gAh~uYf|OTPHhl{LpXT`Oj10t53j*@M!z}3+K*`d1tybJ&=oPrp#N&du*Je z{l=s7p%x;pAUl0}-1gbS8!w$bc5ue~%&kdCy(HkODhi5LPDE~>J+-&TF&x(J9g2*0 z1*dzBYty}hUF|Wy6{5^2gwpfT@N`)9)RD-o{hqa1$MmFmW7Yf2O8cFu*i+{Qm*?~$ zy|_1^zi?>w)bjMuSkNC)j{4e8&xOzJ={hmjKR4lxw&~gU3g8|PtFqwQV3W4QeX`Ss zV^=RPh9h<+v0lw+I2Lbv zYOC8V5PC`2xqa)l9ozD9vVnis(Pj+#Y)-c#>}mVOwHup%_;PBf58Hr9PfaW?$r#fM zp4@;n$Y0I`#pQ4h;xyBb99%nn;?R>f4*&hLx4u5)|NY_Rr)G~Y4=ncE$0Nq6h<>q0cYG|c zG}y^wG=kF1AJQFPjlKKK)Ej4JfBA6Z-sRbwD^XyY{o=~Riz_EPy_~6T%W}Xq;@3}g z+LlIw6T{9}P-oKdB{JH|k+uGI`PmWM<#F>06otW(OvEtQ)%3)2e^%%I0^foM9Yv zzI8{^<4||ks<{rku*<8OiFvOr%)j~awO3!8J-Rf=BQ>TaCuXH4WECa0$;!VuSkIpV z{i^uLzD^}jdRp)TARFu11O1`7;huoaY*w(M0k6|0l`&x+R8B)9+;+EK%`VxM=P=kF zsXyQPL8Q$v+-sN~vro;L5A>=R-RiS@IyS%h-$3xae|&hk$F|ySSnRYd_BvL_yl3Yl z3quf&(b!m>*ZPajn!86saZoyW_7?fe}MDIv{t4Z3fQT5ou z!=t^Y&K&L^vC~)-Dh|b`kYfS8OTeD*b?ohR461n+4%)y%`?O4xfY{msNN18!NohB8 zPOL@sq=!O!9Ro)bGO2(tcCk_l8oD3hCc7?2b4i)RM=q|0Fxqt7*kDh<(XK%l8 z_wMQXdQ|eR$8xgr)nWpe1Y=?QzM<&LXI6i9?ZmSehQ~X7L|jc~Qlb!pp0W4Nj`qLy z+|uTkR{}14w?%XS658k(1&@k{-Y&L&DC`E1bq_S3q@18r?X@Wtz?W$(Yb+^__IAhT z7mLb^lxoqPhi}H`hbr@vTw3v3EO=vJ_~g>c$%{v>ogUcy>>)^0_ITV@Rh>(N@vHH1 z1Et@_+L(5oSsn>m#I(jLt4{jd{e6G<=)#XL9sl&B^@lGkj)auZ=^Hn79$t%GJ=YP@ zn4;a*z8-tPsMs?e91TF#2)vVGUK8*6H}0KZnm8D>PWd#ERt67O6YR5IIJNZD(&@R; z)tOU03kzKv15?jlyVe&`QyVEQ3>=V)l>8Q}h_`3N8?Zq(QjeI?tK%OS^llu9_?==Qxj{*3sx#uTDpTnk%Pn3dLqEL)An{`US`KYsET|G$ENSf^8oDphQxLs;(?HaMh+gZ<`< zt0PgLQjSF)o$h?&smrf^dgI*9+qrv3r{WHU zkVC*$3uuj0JTxeF%UWtY25!*Ehz68*&yQa^&@Cr5!1t$J-grB{aq^(Q)5h~#L|y}P zJfs?KYt@R$fJ7FH@PG&MnZzCiMMMT{b!DGV`^JOk&mCE9x5@oZ<%6fMKls7T7NAb- z&XqG61~qFuW*-Xco;cpI`F~z|_ou?bXs4WYK$trNU-<>(KeJiNYd6w?4euxxPcskjU<{`-`g1N&kFy?ur@5gkW- z7ofY3Hc&1B7~#Kfzk6+6$8A<}YL$$db}2q2X^I(|&&51X9P2k3*kvWzkPx#v;eGar zJ*QTdUcbMxdLV+R&CM^YY{WPA?;DJ^clJ2cBVq6Wnqof9P`hiW!{gG(90pO)*}8XN z*x$2Rp#TUPR@8NH+yz$1lUp)QPa96uTP!Adpep4-zUW2PG zuPQ4-ROeTh7q#;1wB))WJ^#{lc+koo)Qj$3{#wQK!I62dGNfyZW+U=F> zpB$cvw*TS5#xs}u5fxS0#Z|eL1)=qek;5xT=k;@82ZVt}OopLO$6z$#(DQvJQLxR> z%EyKkEe8hNz;6s&`IySQhPr$fD~np0$*8En)EC#HQkxJ3HW9DisSY~jR6v0gXJu!k zra|uW&h0hG0yP(<5H^|xB&UKf>Jne*bNuqTlRvmTwi?yOeO#}G3gR@2oYHQfy9_OS zP7`FKIJERuKH02i@o6X#8L@Y)>%F&L-#E5(bZOxyuiu6H-LHS~l9)#4@o6D{Ys}*b z`xQ@~2)}+TdMl>AJRh_uoAZ*>I~>+eKe_+bE6;cPr8*W-#ApQg6~J)7EG(cj@)&3q zrp_dwPW78s=KL;;m`!VJAs`72C>92_JRf-bU+BhufIHBpbOC1r5T($Qm3Q)=-sJ}BSQm`0Dz-?1`+6>#s-64tJ%ry4n@5-2Hc4S&T24lKVO~y7L3&nZN=_<-bY>P6=E(#U0HX2%I#rOCm7bEHom5{@ z$s}QcX&k6{Vv$I-=)e$3SDXP0T)p233?}@qi*|nkZPFYaT1u?Kq)&$T$t3j$3 zG>Rz524I^vBh4Z@3tQ$;P+ooR*$;m5uOEK!-nR-gr^iO4A=-`Ffp%xRSF(4;``XFj zV^B!xH$)xkvdq+>zDG`kKfeF&%A!5)H^Bk-ngmg=4ya<_nU%_$WK>p%PdO5BPmTCN z+iwsF#S$i$NwCO3{3AFuJ8|RqsZ$5KI%Lgi0!o0ZAR!8hbCYTd^6JaWYHG?L4N8D( zs4mMX&dq{Vb3tAbxuJ|p00ijnRt{yrA_n`lh(k?I-CCTT-szL5 Float[ndarray, \"{len(self.array)}\"]:\n", - " return np.random.rand(self.length)\n", - "\n", - "abc = ABC(3)\n", - "print(abc.length)\n", - "print(abc.check())" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[array([[0.90476524, 0.10874211, 0.24224782, 0.43911673],\n", - " [0.45333331, 0.56400373, 0.27479027, 0.91849912],\n", - " [0.73771259, 0.66894545, 0.44491045, 0.41501535]]),\n", - " array([[0.52512411, 0.08684256, 0.55091178, 0.61461792, 0.17619967,\n", - " 0.18451527],\n", - " [0.29144337, 0.54634595, 0.95109858, 0.38502176, 0.20480036,\n", - " 0.84316217],\n", - " [0.9438109 , 0.04602024, 0.4003533 , 0.64806409, 0.27990322,\n", - " 0.55300185],\n", - " [0.2204828 , 0.1696394 , 0.39441362, 0.22283286, 0.30059866,\n", - " 0.17681949],\n", - " [0.84324772, 0.13621688, 0.33615577, 0.79379001, 0.53483733,\n", - " 0.18615876]])]" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "@jaxtyped(typechecker=beartype)\n", - "def check(a: List[Float[ndarray, \"n m\"]]):\n", - " return a\n", - "\n", - "check([np.random.rand(3, 4), np.random.rand(5, 6)])" - ] - }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { - "name": "stdout", + "name": "stderr", "output_type": "stream", "text": [ - "boxes: torch.Size([9, 5])\n", - "gbbs: torch.Size([9, 3])\n", - "tensor(0.0357) tensor(0.0295) tensor(0.0031)\n", - "boxes: (9, 5)\n", - "gbbs: (9, 3)\n", - "0.035692394 0.02945243 0.0031444672\n" + "100%|██████████| 10/10 [00:05<00:00, 1.84it/s]\n" ] } ], "source": [ - "import torch\n", - "boxes = torch.rand(9, 5)\n", - "print(f\"boxes: {boxes.shape}\")\n", - "gbbs = torch.cat((boxes[:, 2:4].pow(2) / 12, boxes[:, 4:]), dim=-1)\n", - "print(f\"gbbs: {gbbs.shape}\")\n", - "a, b, c = gbbs.split(1, dim=-1)\n", - "cos = c.cos()\n", - "sin = c.sin()\n", - "cos2 = cos.pow(2)\n", - "sin2 = sin.pow(2)\n", - "ans = a * cos2 + b * sin2, a * sin2 + b * cos2, (a - b) * cos * sin\n", - "print(ans[0].mean(), ans[1].mean(), ans[2].mean())\n", "\n", - "boxes = boxes.numpy()\n", - "print(f\"boxes: {boxes.shape}\")\n", - "gbbs = np.concatenate((boxes[:, 2:4] ** 2 / 12, boxes[:, 4:]), axis=-1)\n", - "print(f\"gbbs: {gbbs.shape}\")\n", - "a, b, c = np.split(gbbs, [1, 2], axis=-1)\n", - "cos = np.cos(c)\n", - "sin = np.sin(c)\n", - "cos2 = cos ** 2\n", - "sin2 = sin ** 2\n", - "ans = a * cos2 + b * sin2, a * sin2 + b * cos2, (a - b) * cos * sin\n", - "print(ans[0].mean(), ans[1].mean(), ans[2].mean())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": 21, - "metadata": {}, - "outputs": [ - { - "data": { - "application/vnd.jupyter.widget-view+json": { - "model_id": "4ae13fc25d9046a381e5e56745a394ed", - "version_major": 2, - "version_minor": 0 - }, - "text/plain": [ - " 0%| | 0/10000 [00:00\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "

<xarray.DataArray (band: 3, y: 128, x: 128)> Size: 49kB\n",
-                            "[49152 values with dtype=uint8]\n",
-                            "Coordinates:\n",
-                            "  * band         (band) int64 24B 1 2 3\n",
-                            "  * x            (x) float64 1kB 7.405e+05 7.405e+05 ... 7.417e+05 7.417e+05\n",
-                            "  * y            (y) float64 1kB 3.741e+06 3.741e+06 ... 3.74e+06 3.74e+06\n",
-                            "    spatial_ref  int64 8B 0\n",
-                            "Attributes: (12/30)\n",
-                            "    AREA_OR_POINT:       Area\n",
-                            "    _FillValue:          0\n",
-                            "    scale_factor:        1.0\n",
-                            "    add_offset:          0.0\n",
-                            "    timestamp:           2023-01-24 16:25:51.024000+00:00\n",
-                            "    visual_href:         https://sentinel2l2a01.blob.core.windows.net/sentine...\n",
-                            "    ...                  ...\n",
-                            "    granule-metadata:    https://sentinel2l2a01.blob.core.windows.net/sentine...\n",
-                            "    inspire-metadata:    https://sentinel2l2a01.blob.core.windows.net/sentine...\n",
-                            "    product-metadata:    https://sentinel2l2a01.blob.core.windows.net/sentine...\n",
-                            "    datastrip-metadata:  https://sentinel2l2a01.blob.core.windows.net/sentine...\n",
-                            "    tilejson:            https://planetarycomputer.microsoft.com/api/data/v1/...\n",
-                            "    rendered_preview:    https://planetarycomputer.microsoft.com/api/data/v1/...
" - ], - "text/plain": [ - " Size: 49kB\n", - "[49152 values with dtype=uint8]\n", - "Coordinates:\n", - " * band (band) int64 24B 1 2 3\n", - " * x (x) float64 1kB 7.405e+05 7.405e+05 ... 7.417e+05 7.417e+05\n", - " * y (y) float64 1kB 3.741e+06 3.741e+06 ... 3.74e+06 3.74e+06\n", - " spatial_ref int64 8B 0\n", - "Attributes: (12/30)\n", - " AREA_OR_POINT: Area\n", - " _FillValue: 0\n", - " scale_factor: 1.0\n", - " add_offset: 0.0\n", - " timestamp: 2023-01-24 16:25:51.024000+00:00\n", - " visual_href: https://sentinel2l2a01.blob.core.windows.net/sentine...\n", - " ... ...\n", - " granule-metadata: https://sentinel2l2a01.blob.core.windows.net/sentine...\n", - " inspire-metadata: https://sentinel2l2a01.blob.core.windows.net/sentine...\n", - " product-metadata: https://sentinel2l2a01.blob.core.windows.net/sentine...\n", - " datastrip-metadata: https://sentinel2l2a01.blob.core.windows.net/sentine...\n", - " tilejson: https://planetarycomputer.microsoft.com/api/data/v1/...\n", - " rendered_preview: https://planetarycomputer.microsoft.com/api/data/v1/..." - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "lat_c = 33.7756\n", - "lon_c = -84.3963\n", - "time_of_interest = \"2023-01-01/2023-01-31\"\n", - "max_cloud_cover = 10.0\n", - "max_items = 10\n", + "def pull_image(x, y):\n", + " params = {\n", + " \"api-version\": \"2.1\",\n", + " \"tilesetId\": \"microsoft.imagery\",\n", + " \"zoom\": 19,\n", + " \"x\": x,\n", + " \"y\": y,\n", + " \"tileSize\": 256,\n", + " }\n", "\n", - "# polygon = box(lon_c - 0.01, lat_c - 0.01, lon_c + 0.01, lat_c + 0.01)\n", - " \n", - "# catalog = Client.open(\"https://planetarycomputer.microsoft.com/api/stac/v1\", modifier=pc.sign_inplace)\n", + " headers = {\n", + " \"Authorization\": \"Bearer eyJ0eXAiOiJKV1QiLCJhbGciOiJSUzI1NiIsIng1dCI6Ikg5bmo1QU9Tc3dNcGhnMVNGeDdqYVYtbEI5dyIsImtpZCI6Ikg5bmo1QU9Tc3dNcGhnMVNGeDdqYVYtbEI5dyJ9.eyJhdWQiOiJodHRwczovL2F0bGFzLm1pY3Jvc29mdC5jb20iLCJpc3MiOiJodHRwczovL3N0cy53aW5kb3dzLm5ldC8xNmIzYzAxMy1kMzAwLTQ2OGQtYWM2NC03ZWRhMDgyMGI2ZDMvIiwiaWF0IjoxNzI2NjAyNzg3LCJuYmYiOjE3MjY2MDI3ODcsImV4cCI6MTcyNjY4OTQ4NywiYWlvIjoiRTJkZ1lKRCtiNmYvOFhQekVXMVdOdGNzMStiUEFBPT0iLCJhcHBpZCI6IjE3NzE1ODE1LTgwNDQtNGZkMi1hMjVkLWViNDVhNzQ1YmEzYSIsImFwcGlkYWNyIjoiMiIsImlkcCI6Imh0dHBzOi8vc3RzLndpbmRvd3MubmV0LzE2YjNjMDEzLWQzMDAtNDY4ZC1hYzY0LTdlZGEwODIwYjZkMy8iLCJpZHR5cCI6ImFwcCIsIm9pZCI6Ijk1MjgwNDczLWU4MzItNGMwMi1iMDlmLTljYjBiNGMyNzZkNyIsInJoIjoiMC5BRVlBRThDekZnRFRqVWFzWkg3YUNDQzIweUtnSHJvSFdOVkJ1LXNwTEg0YzlmYnhBQUEuIiwic3ViIjoiOTUyODA0NzMtZTgzMi00YzAyLWIwOWYtOWNiMGI0YzI3NmQ3IiwidGlkIjoiMTZiM2MwMTMtZDMwMC00NjhkLWFjNjQtN2VkYTA4MjBiNmQzIiwidXRpIjoiNVJfdXRNZGlBVUtOVEpsRDI5QWRBQSIsInZlciI6IjEuMCIsInhtc19pZHJlbCI6IjcgOCJ9.TKZNz289jA_APg_PKEeSvTX77EnPJ4yrvXcgL4m3c2um2UeWeTztj-q5ynPBY9R98qht3JrcQSh6iZIQorWqTX3mdFcXCe7TTAXMk-2j0dpYWL9fn_mNswAJY4I5QhGX-asGPSqcEiOWcogxsGRCuFmQxkhJy5S9IUi0lP9tipdsuCxvV-9M3EKuntyQVur8Ha7pjM9FJYNsBrEFleoP_if5FNthH_TPqpeRod6rQ8u2kak7RhkSAi9_G13V4pGuhqlzYdSQlx8NPqXIOdYHElThXnsZWsB-RpdWoLwvPEYu08kDdDi-IGzQK-N_UsDatXJDobUxFFBo-embjWj8yg\",\n", + " \"X-MS-Client-Id\": \"d069e722-70c3-4dd6-8532-a6f4b18c9bfb\",\n", + " }\n", "\n", - "# search = catalog.search(\n", - "# collections=[\"sentinel-2-l2a\"],\n", - "# intersects=polygon,\n", - "# datetime=time_of_interest,\n", - "# query={\"eo:cloud_cover\": {\"lt\": max_cloud_cover}},\n", - "# max_items=max_items,\n", - "# )\n", + " response = requests.get(url, params=params, headers=headers)\n", + " image_data = response.content\n", + " image = Image.open(io.BytesIO(image_data))\n", + " return np.array(image)\n", "\n", - "# items = search.item_collection()\n", - "# least_cloud_cover_item = min(items, key=lambda x: eo.ext(x).cloud_cover)\n", - "# href = least_cloud_cover_item.assets[\"visual\"].href\n", - "# visual_raster = rioxarray.open_rasterio(pc.sign(href))\n", - "# visual_raster\n", - "ds = get_sentinel2_visual(lat_c, lon_c, 128, 128, time_of_interest, max_cloud_cover)\n", - "ds" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "CRS.from_epsg(32616)" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "ds.rio.crs" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "ds.rio.to_raster(\"test.tif\")" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "EPSG:32616\n", - "{'AOT': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_AOT_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'AREA_OR_POINT': 'Area', 'B01': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R60m/T16SGC_20230124T162551_B01_60m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B02': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_B02_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B03': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_B03_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B04': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_B04_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B05': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_B05_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B06': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_B06_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B07': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_B07_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B08': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_B08_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B09': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R60m/T16SGC_20230124T162551_B09_60m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B11': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_B11_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B12': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_B12_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'B8A': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_B8A_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'datastrip-metadata': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/DATASTRIP/DS_ESRI_20230126T150710_S20230124T162913/MTD_DS.xml?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'granule-metadata': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/MTD_TL.xml?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'inspire-metadata': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/INSPIRE.xml?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'preview': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/QI_DATA/T16SGC_20230124T162551_PVI.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'product-metadata': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/MTD_MSIL2A.xml?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'rendered_preview': 'https://planetarycomputer.microsoft.com/api/data/v1/item/preview.png?collection=sentinel-2-l2a&item=S2A_MSIL2A_20230124T162551_R040_T16SGC_20230126T150708&assets=visual&asset_bidx=visual%7C1%2C2%2C3&nodata=0&format=png', 'safe-manifest': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/manifest.safe?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'SCL': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R20m/T16SGC_20230124T162551_SCL_20m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'tilejson': 'https://planetarycomputer.microsoft.com/api/data/v1/item/tilejson.json?collection=sentinel-2-l2a&item=S2A_MSIL2A_20230124T162551_R040_T16SGC_20230126T150708&assets=visual&asset_bidx=visual%7C1%2C2%2C3&nodata=0&format=png', 'timestamp': '2023-01-24 16:25:51.024000+00:00', 'visual': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_TCI_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'visual_href': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_TCI_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D', 'WVP': 'https://sentinel2l2a01.blob.core.windows.net/sentinel2-l2/16/S/GC/2023/01/24/S2A_MSIL2A_20230124T162551_N0400_R040_T16SGC_20230126T150708.SAFE/GRANULE/L2A_T16SGC_A039648_20230124T162913/IMG_DATA/R10m/T16SGC_20230124T162551_WVP_10m.tif?st=2024-08-24T10%3A36%3A57Z&se=2024-08-25T11%3A21%3A57Z&sp=rl&sv=2024-05-04&sr=c&skoid=9c8ff44a-6a2c-4dfb-b298-1c9212f64d9a&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2024-08-25T02%3A23%3A12Z&ske=2024-09-01T02%3A23%3A12Z&sks=b&skv=2024-05-04&sig=auUiWOkdO/NOaTlBNCWm/5KMLyYFd2phAcVDusZCHRI%3D'}\n" - ] - } - ], - "source": [ - "loaded_tif = rio.open(\"test.tif\")\n", - "# get crs\n", - "print(loaded_tif.crs)\n", - "# get timestamp\n", - "print(loaded_tif.tags())" + "np_list = []\n", + "for x in tqdm(range(380000, 380010)):\n", + " sleep_time = np.random.uniform(0.1, 0.5)\n", + " # sleep(sleep_time)\n", + " buffer_list = []\n", + " for y in range(221700, 221710):\n", + " buffer_list.append(delayed(pull_image)(x, y))\n", + " results = Parallel(n_jobs=16)(buffer_list)\n", + " np_list.append(results)\n", + " \n", + "transformed_list = tree_map(lambda x: rearrange(x, 'h w c -> c w h'), np_list)\n", + "big_img = np.block(transformed_list)\n", + "big_img = rearrange(big_img, 'c h w -> h w c')\n", + "print(big_img.shape)" ] } ], diff --git a/requirements.txt b/requirements.txt index c406cda..b172bd8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,4 +8,4 @@ planetary-computer shapely rioxarray opencv-python -supervision \ No newline at end of file +einops \ No newline at end of file diff --git a/tests/test_ops.py b/tests/test_utils.py similarity index 98% rename from tests/test_ops.py rename to tests/test_utils.py index 877327c..866b778 100644 --- a/tests/test_ops.py +++ b/tests/test_utils.py @@ -4,7 +4,7 @@ from ultralytics.utils import ops from ultralytics.utils.metrics import probiou -from garuda.ops import webm_pixel_to_geo, geo_to_webm_pixel, obb_to_aa, label_studio_csv_to_obb, obb_iou, xyxyxyxy2xywhr, xywhr2xyxyxyxy +from garuda.utils import webm_pixel_to_geo, geo_to_webm_pixel, obb_to_aa, label_studio_csv_to_obb, obb_iou, xyxyxyxy2xywhr, xywhr2xyxyxyxy # test zoom levels test_zoom_levels = [0, 8, 17, 20]