Merge pull request #217 from Radio-Spectrum/feat/random_grid_transform

feat: service grid random transformation as parameter option

Author: brunohcfaria

sharc/parameters/imt/parameters_imt_mss_dc.py

@@ -136,6 +136,8 @@ class ParametersServiceGrid(ParametersBase):
 
         country_names: list[str] = field(default_factory=lambda: list([""]))
 
+        transform_grid_randomly: bool = False
+
         # margin from inside of border [km]
         # if positive, makes border smaller by x km
         # if negative, makes border bigger by x km
@@ -197,7 +199,7 @@ def validate(self, ctx: str):
                 raise ValueError(
                     f"{ctx}.eligible_sats_margin_from_border needs to be a number")
 
-            self.reset_grid(ctx)
+            self._load_geom_from_file_if_needed(ctx)
 
             super().validate(ctx)
 
@@ -218,13 +220,30 @@ def load_from_active_sat_conditions(
             if self.eligible_sats_margin_from_border is None:
                 self.eligible_sats_margin_from_border = sat_is_active_if.lat_long_inside_country.margin_from_border
 
-        def reset_grid(self, ctx: str, force_update=False):
+        def reset_grid(
+            self,
+            ctx: str,
+            rng: np.random.RandomState,
+            force_update=False,
+        ):
             """
             After creating grid, there are some features that can only be implemented
-            with knowledge of other parts of the simulator. This method's purpose is
-            to run only once at the start of the simulation
+            with knowledge of other parts of the simulator.
             """
-            if self.lon_lat_grid is not None and not force_update:
+            self._load_geom_from_file_if_needed(ctx, force_update)
+
+            self.lon_lat_grid = generate_grid_in_multipolygon(
+                self.grid_borders_polygon,
+                self.beam_radius,
+                self.transform_grid_randomly,
+                rng
+            )
+
+            self.ecef_grid = lla2ecef(
+                self.lon_lat_grid[1], self.lon_lat_grid[0], 0)
+
+        def _load_geom_from_file_if_needed(self, ctx: str, force_update=False):
+            if self.eligibility_polygon is not None and not force_update:
                 return
             filtered_gdf = load_gdf(
                 self.country_shapes_filename,
@@ -239,17 +258,13 @@ def reset_grid(self, ctx: str, force_update=False):
             shrinked = shrink_countries_by_km(
                 filtered_gdf.geometry.values, self.grid_margin_from_border
             )
-            polygon = shp.ops.unary_union(shrinked)
-            assert polygon.is_valid, shp.validation.explain_validity(polygon)
-            assert not polygon.is_empty, "Can't have a empty polygon as filter"
+            self.grid_borders_polygon = shp.ops.unary_union(shrinked)
 
-            self.lon_lat_grid = generate_grid_in_multipolygon(
-                polygon,
-                self.beam_radius
-            )
+            assert self.grid_borders_polygon.is_valid, \
+                shp.validation.explain_validity(self.grid_borders_polygon)
 
-            self.ecef_grid = lla2ecef(
-                self.lon_lat_grid[1], self.lon_lat_grid[0], 0)
+            assert not self.grid_borders_polygon.is_empty, \
+                "Can't have a empty grid_borders_polygon as filter"
 
             self.eligibility_polygon = shp.ops.unary_union(shrink_countries_by_km(
                 filtered_gdf.geometry.values, self.eligible_sats_margin_from_border

sharc/satellite/utils/sat_utils.py

@@ -115,6 +115,25 @@ def calc_elevation(Le: np.ndarray,
     return np.degrees(elev_angle)
 
 
+def haversine(
+    lon1: np.ndarray,
+    lat1: np.ndarray,
+    lon2: np.ndarray,
+    lat2: np.ndarray,
+    R=EARTH_RADIUS_M
+):
+    """Calculates great-circle distance between 2 points on the surface of Earth.
+    Considers spherical earth.
+
+    Returns np.ndarray with N distances in meters
+    """
+    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
+    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2
+    return 2 * R * np.arcsin(np.sqrt(a))
+
+
 if __name__ == "__main__":
     r1 = ecef2lla(7792.1450, 0, 0)
     print(r1)

sharc/support/sharc_geom.py

@@ -608,6 +608,8 @@ def shrink_countries_by_km(
 def generate_grid_in_polygon(
     polygon: shp.geometry.Polygon,
     hexagon_radius: float,
+    rotation_deg: typing.Optional[float] = None,
+    translation: typing.Optional[typing.Tuple[float, float]] = None,
 ):
     """Generate a hexagonal grid inside a polygon and return points in EARTH_DEFAULT_CRS.
 
@@ -621,6 +623,11 @@ def generate_grid_in_polygon(
         Polygon to fill with a grid.
     hexagon_radius : float
         Radius of the hexagons in the grid (in meters).
+    rotation_deg : float or None, optional
+        Rotation of the grid in degrees (default: None, no rotation).
+    translation : tuple of float or None, optional
+        Translation of the grid in meters (dx, dy) (default: None, no translation).
+        The translation must not exceed the grid spacing in magnitude.
 
     Returns
     -------
@@ -644,14 +651,18 @@ def generate_grid_in_polygon(
     # Transform to projection where unit is meters
     polygon_proj = shp.ops.transform(to_proj, polygon)
 
-    # create a bounding box, afterwards we filter to polygon site
-    minx, miny, maxx, maxy = polygon_proj.bounds
+    # Determine x/y spacing
     x_spacing = 3 * hexagon_radius
     y_spacing = hexagon_radius * np.sqrt(3) / 2
 
-    x_vals = np.arange(minx, maxx + x_spacing, x_spacing)
+    # Create a bounding circle, afterwards we filter to polygon site
+    cx, cy = polygon_proj.centroid.coords[0]
+    minx, miny, maxx, maxy = polygon_proj.bounds
+    bbox_diag = np.hypot(maxx - minx, maxy - miny)
+    bound_radius = bbox_diag / 2 + 3 * hexagon_radius
 
-    y_vals = np.arange(miny, maxy + y_spacing, y_spacing)
+    x_vals = np.arange(cx - bound_radius, cx + bound_radius + x_spacing, x_spacing)
+    y_vals = np.arange(cy - bound_radius, cy + bound_radius + y_spacing, y_spacing)
 
     x_vals, y_vals = np.meshgrid(x_vals, y_vals)
 
@@ -661,6 +672,20 @@ def generate_grid_in_polygon(
     x_vals = x_vals.ravel()
     y_vals = y_vals.ravel()
 
+    if rotation_deg is not None:
+        theta = np.deg2rad(rotation_deg)
+        cos_t, sin_t = np.cos(theta), np.sin(theta)
+        x_rot = cos_t * (x_vals - cx) - sin_t * (y_vals - cy) + cx
+        y_rot = sin_t * (x_vals - cx) + cos_t * (y_vals - cy) + cy
+        x_vals, y_vals = x_rot, y_rot
+
+    if translation is not None:
+        dx, dy = translation
+        if abs(dx) > x_spacing or abs(dy) > y_spacing:
+            raise ValueError("generate_grid_in_polygon.translation (dx, dy) must not exceed grid spacing (x_spacing, y_spacing) in magnitude")
+        x_vals += dx
+        y_vals += dy
+
     # Return to EARTH_DEFAULT_CRS
     xt, yt = from_proj(x_vals, y_vals)
 
@@ -676,34 +701,66 @@ def generate_grid_in_polygon(
 
 def generate_grid_in_multipolygon(
     poly: typing.Union[shp.geometry.MultiPolygon, shp.geometry.Polygon],
-    km: float
+    km: float,
+    random_transform_on_grid: bool = False,
+    rng: np.random.RandomState = None,
 ) -> list[shp.geometry.MultiPolygon]:
-    """Generate a grid in a MultiPolygon or Polygon, shrinking each by a given number of kilometers.
+    """Generate a hexagonal grid in a MultiPolygon or Polygon,
+    considering a hexagon radius in km.
 
-    For each polygon, create a grid and return a single 2xN array of longitudes and latitudes.
+    For each polygon, create a grid and return a single 2xN array of longitudes and latitudes
+    containing all grids.
 
     Parameters
     ----------
     poly : typing.Union[shp.geometry.MultiPolygon, shp.geometry.Polygon]
         The MultiPolygon or Polygon to process.
     km : float
-        Number of kilometers to shrink each polygon by.
+        Hexagon radius in km
+    random_transform_on_grid : bool, optional
+        Whether to apply a random rotation and translation to the grid (default: False).
+    rng : np.random.RandomState, optional
+        Random number generator to use if random_transform_on_grid is True (default: None).
 
     Returns
     -------
     np.ndarray
         2xN array: first row is longitudes, second row is latitudes.
     """
+    if random_transform_on_grid:
+        assert rng is not None
+
     lons = []
     lats = []
 
     if poly.geom_type == 'Polygon':
-        x, y = generate_grid_in_polygon(poly, km)
+        if random_transform_on_grid:
+            x, y = generate_grid_in_polygon(
+                poly, km,
+                rng.uniform(-180., 180.),
+                (
+                    3 * rng.uniform(-km, km),
+                    rng.uniform(-km, km) * np.sqrt(3) / 2
+                )
+            )
+        else:
+            x, y = generate_grid_in_polygon(poly, km)
         lons.extend(x)
         lats.extend(y)
     elif poly.geom_type == 'MultiPolygon':
         for p in poly.geoms:
-            x, y = generate_grid_in_polygon(p, km)
+            if random_transform_on_grid:
+                x, y = generate_grid_in_polygon(
+                    p, km,
+                    rng.uniform(-180., 180.),
+                    (
+                        3 * rng.uniform(-km, km),
+                        rng.uniform(-km, km) * np.sqrt(3) / 2
+                    )
+                )
+            else:
+                x, y = generate_grid_in_polygon(p, km)
+
             lons.extend(x)
             lats.extend(y)
 

sharc/topology/topology_imt_mss_dc.py

@@ -373,7 +373,10 @@ def get_satellite_pointing(
                 (before sharc's coordinate transformation, with distances in km)
         """
         if orbit_params.beam_positioning.type == "SERVICE_GRID":
-            orbit_params.beam_positioning.service_grid.validate("service_grid")
+            # TODO: remove this reset_grid call from here.
+            # Since it should be called once per drop, it shouldn't
+            # be this deep in the program
+            orbit_params.beam_positioning.service_grid.reset_grid("service_grid", random_number_gen)
 
             # 2xN, (lon, lat)
             grid = orbit_params.beam_positioning.service_grid.lon_lat_grid

tests/parameters/parameters_for_testing.yaml

@@ -232,6 +232,7 @@ imt:
                         min: 0.0
                         max: 66.1
                 service_grid:
+                    transform_grid_randomly: true # add per drop rand rotation + transl. to grid
                     # by default this already gets shapefile from natural earth
                     country_shapes_filename: sharc/data/countries/ne_110m_admin_0_countries.shp
                     # By default, these are the same as taken from `sat_is_active_if.lat_long_inside_country`
@@ -1008,6 +1009,7 @@ mss_d2d:
                 min: 0.0
                 max: 66.1
         service_grid:
+            transform_grid_randomly: true # add per drop rand rotation + transl. to grid
             # by default this already gets shapefile from natural earth
             country_shapes_filename: sharc/data/countries/ne_110m_admin_0_countries.shp
             # By default, these are the same as taken from `sat_is_active_if.lat_long_inside_country`

tests/parameters/test_parameters.py

@@ -280,6 +280,9 @@ def test_parameters_imt(self):
         self.assertEqual(
             self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.beam_radius,
             19000)
+        self.assertEqual(
+            self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.transform_grid_randomly,
+            True)
         self.assertEqual(
             self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.grid_margin_from_border,
             0.11)
@@ -630,6 +633,9 @@ def test_parametes_mss_d2d(self):
         self.assertEqual(
             self.parameters.mss_d2d.beam_positioning.service_grid.beam_radius,
             19001)
+        self.assertEqual(
+            self.parameters.mss_d2d.beam_positioning.service_grid.transform_grid_randomly,
+            True)
         self.assertEqual(
             self.parameters.mss_d2d.beam_positioning.service_grid.grid_margin_from_border,
             0.11)