Merge pull request #60 from NOAA-GSL/bug/idsse-638/test-syracuse

Bug fix to vectaster

python/idsse_common/idsse/common/sci/geo_image.py

@@ -41,10 +41,18 @@ class ColorPalette(NamedTuple):
     under_idx: int
     over_idx: int
     fill_idx: int
+    min_value: float = None
+    max_value: float = None
 
     # pylint: disable=invalid-name
     @classmethod
-    def linear(cls, colors: Sequence[Color], anchors: Sequence[int] = None) -> Self:
+    def linear(
+        cls,
+        colors: Sequence[Color],
+        anchors: Sequence[int] = None,
+        min_value: float = None,
+        max_value: float = None
+    ) -> Self:
         """Create a color palette by linearly interpolating between colors
 
         Args:
@@ -53,6 +61,8 @@ def linear(cls, colors: Sequence[Color], anchors: Sequence[int] = None) -> Self:
             anchors (Sequence[int], optional): A list to define what index each color
                                       should be mapped to. Required: anchors[0] must equal 0,
                                       and anchor[-1] must 255. Defaults to None.
+            min_value (float, optional): If provided, used to shift the anchors down.
+            max_value (float, optional): If provided, used to scale the anchors.
 
         Raises:
             ValueError: Raised when the length of the colors and anchors do not match
@@ -64,13 +74,17 @@ def linear(cls, colors: Sequence[Color], anchors: Sequence[int] = None) -> Self:
         if anchors is not None:
             if len(anchors) != num:
                 raise ValueError('Colors and Anchors must be of the same length')
+            if min_value:
+                anchors = [x-min_value for x in anchors]
+            if max_value:
+                anchors = [x/max_value*255 for x in anchors]
             xp = anchors
         else:
             xp = [round_(pos, rounding='floor') for pos in np.linspace(0, 255, num=num)]
         lut = list(
             (round_(r, 0, 'floor'), round_(g, 0, 'floor'), round_(b, 0, 'floor')) for (r, g, b) in
             zip(*list(np.interp(range(256), xp, fp) for fp in np.array(colors).T)))
-        return ColorPalette(lut, 256, 0, 255, 0)
+        return ColorPalette(lut, 256, 0, 255, 0, min_value, max_value)
 
     @classmethod
     def grey(cls, with_excludes: bool = True) -> Self:
@@ -100,6 +114,28 @@ def __init__(self, proj: GridProj, rgb_array: np.ndarray, scale: int = 1):
         self.rgb_array = rgb_array
         self.scale = scale
 
+    @classmethod
+    def from_proj(
+        cls,
+        proj: GridProj,
+        fill_color: Color = (255, 255, 255),
+        scale: int = 1
+    ) -> Self:
+        """Method for building a geographical image without background data.
+
+        Args:
+            proj (GridProj): Grid projection to be used for this geo image, must match data_grid.
+            fill_color (Color): The color to fill the image with. Default is white.
+            scale (int, optional): The height and width that a grid cell will be scaled to in the
+                                   image. Defaults to 1.
+
+        Returns:
+            Self: GeoImage
+        """
+        rgb_array = np.zeros((proj.width*scale, proj.height*scale, 3), np.uint8)
+        rgb_array[...] = fill_color
+        return GeoImage(proj, rgb_array, scale)
+
     @classmethod
     def from_index_grid(
         cls,
@@ -167,6 +203,10 @@ def from_data_grid(  # pylint: disable=too-many-arguments
         """
         if colors is None:
             colors = ColorPalette.grey()
+        if not min_value:
+            min_value = colors.min_value
+        if not max_value:
+            max_value = colors.max_value
         norm_array = normalize(data_array, min_value, max_value, fill_value)
         index_array = scale_to_color_palette(norm_array, colors.num_colors)
         return GeoImage.from_index_grid(proj, index_array, colors, scale)
@@ -276,15 +316,15 @@ def draw_linestring(
 
     def draw_shape(
         self,
-        shape: Geometry,
+        shape: Geometry | str | dict,
         color: Color,
         geo: bool = True
     ):
         """Draw a shape onto the image, points and lines will only be on image cell wide
         and polygon will be filled, doesn't make use of scale.
 
         Args:
-            shape (Geometry): A Shapely geometry
+            shape (Geometry): A Shapely geometry. If string, must be wkt; if dict, must be geojson.
             color (Color): RGB value as a tuple of three values between 0 and 255
             geo (bool): Indication that the shape is specified by geographic coordinates (lon/lat).
                         Defaults to True.
@@ -295,6 +335,9 @@ def draw_shape(
         if isinstance(shape, str):
             shape = from_wkt(shape)
 
+        if isinstance(shape, dict):
+            shape = from_geojson(json.dumps(shape))
+
         if geo:
             shape = geographic_to_pixel(shape, self.proj)
 

python/idsse_common/idsse/common/sci/vectaster.py

@@ -9,13 +9,16 @@
 #
 # ----------------------------------------------------------------------------------
 
+import json
 import logging
 from collections.abc import Iterable, Sequence
 from math import floor
 from numbers import Number
 
 import numpy
-from shapely import Geometry, LinearRing, LineString, MultiPolygon, Point, Polygon, from_wkt
+from shapely import (Geometry, LinearRing, LineString,
+                     MultiPolygon, Point, Polygon,
+                     from_geojson, from_wkt)
 
 from idsse.common.sci.grid_proj import GridProj
 from idsse.common.sci.utils import Pixel, Coord, Coords, coordinate_pairs_to_axes
@@ -100,10 +103,10 @@ def rasterize_point(
 
     if grid_proj is not None:
         return coordinate_pairs_to_axes([grid_proj.map_geo_to_pixel(*coord, rounding)],
-                                         dtype=numpy.int64)
+                                        dtype=numpy.int64)
 
     return coordinate_pairs_to_axes([tuple(round_values(*coord, rounding=rounding))],
-                                     dtype=numpy.int64)
+                                    dtype=numpy.int64)
 
 
 def rasterize_linestring(
@@ -188,6 +191,9 @@ def rasterize_polygon(
     if isinstance(polygon, str):
         polygon = from_wkt(polygon)
 
+    if isinstance(polygon, dict):
+        polygon = from_geojson(json.dumps(polygon))
+
     if isinstance(polygon, Polygon):
         coords = [list(interior.coords) for interior in polygon.interiors]
         coords.insert(0, list(coord for coord in polygon.exterior.coords))
@@ -368,9 +374,13 @@ def _pixels_for_linestring(
         list[Pixel]: List of x,y tuples in pixel space
     """
     coords = linestring.coords
-    pixels = _pixels_for_line_seg(coords[0], coords[1])
+    pixels = set(_pixels_for_line_seg(coords[0], coords[1]))
     for pnt1, pnt2 in zip(coords[1:-1], coords[2:]):
-        pixels.extend(_pixels_for_line_seg(pnt1, pnt2, exclude_first=True))
+        pixels.update(_pixels_for_line_seg(pnt1, pnt2, exclude_first=True))
+
+    pixels = list(pixels)
+    pixels.sort()
+
     return pixels
 
 
@@ -406,14 +416,17 @@ def _pixels_for_line_seg(
         step = 1 if x1 < x2 else -1
         if exclude_first:
             x1 += step
-        pixels = [(x, round_(slope * x + intercept)) for x in range(x1, x2+step, step)]
+        pixels = set((x, round_(slope * x + intercept)) for x in range(x1, x2+step, step))
     else:
         slope = dx / dy
         intercept = x1 - slope * y1
         step = 1 if y1 < y2 else -1
         if exclude_first:
             y1 += step
-        pixels = [(round_(slope * y + intercept), y) for y in range(y1, y2+step, step)]
+        pixels = set((round_(slope * y + intercept), y) for y in range(y1, y2+step, step))
+
+    pixels = list(pixels)
+    pixels.sort()
 
     return pixels
 
@@ -435,7 +448,6 @@ def _pixels_for_polygon(
     x_offset = 0 if xmin >= 0 else int(1-xmin)
     ymin = floor(ymin)
     ymax = floor(ymax)
-
     edge_info = []
     for (x1, y1), (x2, y2) in zip(polygon_boundary.coords[:-1], polygon_boundary.coords[1:]):
         int_y1, int_y2 = int(y1), int(y2)
@@ -447,7 +459,7 @@ def _pixels_for_polygon(
 
             edge_info.append((int_y1, int_y2, x1+x_offset, (x2 - x1)/(y2 - y1)))
 
-    pixels = []
+    pixels = set()
     for y in range(ymin, ymax+1):
         x_list = []
         for y1, y2, x, slope in edge_info:
@@ -457,7 +469,12 @@ def _pixels_for_polygon(
         x_list.sort()
 
         for x1, x2 in zip(x_list[:-1], x_list[1:]):
-            pixels.extend([(x, y) for x in range(x1-x_offset, x2-x_offset+1)])
+            pixels.update([(x, y) for x in range(x1-x_offset, x2-x_offset+1)])
+
+    # make sure the edges are included in list of pixels
+    pixels.update(_pixels_for_linestring(polygon_boundary))
+    pixels = list(pixels)
+    pixels.sort()
 
     return pixels
 

python/idsse_common/test/sci/test_geo_image.py

@@ -26,6 +26,19 @@ def proj():
     return GridProj.from_proj_grid_spec(proj_spec, grid_spec)
 
 
+def test_geo_image_without_data_grid(proj):
+    fill_color = (255, 0, 0)
+    geo_image = GeoImage.from_proj(proj, fill_color)
+    values, _, counts = numpy.unique(geo_image.rgb_array,
+                                     return_inverse=True,
+                                     return_counts=True)
+    # the only value in the grid should be 0 and 255, where there should be twice as many
+    # 0 as 255 and the total count should be 11234895 (proj width*height)
+    numpy.testing.assert_array_equal(values, [0, 255])
+    numpy.testing.assert_array_equal(counts, [7489930, 3744965])
+    assert sum(counts) == 11234895
+
+
 def test_geo_image_from_data_grid(proj):
     data = numpy.array([[0, 1, 2, 3],
                         [4, 5, 6, 7],
@@ -182,8 +195,8 @@ def test_draw_polygon(proj):
 
     # values will be 0 or 100 (for polygon) and 0 everywhere else
     numpy.testing.assert_array_equal(values, [0, 100])
-    numpy.testing.assert_array_equal(counts, [237, 6])
-    expected_indices = [91, 94, 118, 121, 124, 145]
+    numpy.testing.assert_array_equal(counts, [236, 7])
+    expected_indices = [91, 94, 118, 121, 124, 145, 148]
     numpy.testing.assert_array_equal(numpy.where(indices == 1)[0], expected_indices)
 
 
@@ -203,9 +216,9 @@ def test_draw_multi_polygon(proj):
 
     # values will be 0 or 100 (for polygon) and 0 everywhere else
     numpy.testing.assert_array_equal(values, [0, 100])
-    numpy.testing.assert_array_equal(counts, [417, 15])
+    numpy.testing.assert_array_equal(counts, [416, 16])
     expected_indices = [118, 121, 124, 154, 157, 160, 190, 193,
-                        196, 235, 238, 271, 274, 277, 307]
+                        196, 235, 238, 271, 274, 277, 307, 310]
     numpy.testing.assert_array_equal(numpy.where(indices == 1)[0], expected_indices)
 
 
@@ -228,11 +241,11 @@ def test_draw_geo_polygon(proj: GridProj):
 
     # values will be 0 or 100 (for polygon) and 0 everywhere else
     numpy.testing.assert_array_equal(values, [0, 100])
-    # numpy.testing.assert_array_equal(counts, [7392, 108])
+    # numpy.testing.assert_array_equal(counts, [[7378, 122])
 
     # the "equality" workarounds below are needed due to counts and indices arrays having
     # different results when run with pytest locally vs. in GitHub Actions runner
-    assert counts.tolist() == approx([7392, 108], rel=0.10)  # counts can be up to 10% off
+    assert counts.tolist() == approx([7378, 122], rel=0.10)  # counts can be up to 10% off
     expected_indices = [2006, 2156, 2306, 2309, 2456, 2459, 2606, 2609, 2753, 2756, 2759,
                         2762, 2903, 2906, 2909, 2912, 3053, 3056, 3059, 3062, 3203, 3206,
                         3209, 3212, 3215, 3353, 3356, 3359, 3362, 3365, 3500, 3503, 3506,
@@ -330,7 +343,7 @@ def test_draw_state(proj):
 
     values, counts = numpy.unique(geo_image.rgb_array, return_counts=True)
     numpy.testing.assert_array_equal(values, [0, 255])
-    numpy.testing.assert_array_equal(counts, [11234304, 591])
+    numpy.testing.assert_array_equal(counts, [11234296, 599])
 
 
 def test_add_one_state(proj):