Merge pull request #368 from ungarj/cli_cp_point

add --point and --point-crs option to cp command

mapchete/cli/default/cp.py

@@ -13,6 +13,8 @@
 @options.opt_area_crs
 @options.opt_bounds
 @options.opt_bounds_crs
+@options.opt_point
+@options.opt_point_crs
 @options.opt_overwrite
 @options.opt_verbose
 @options.opt_no_pbar

mapchete/cli/options.py

@@ -329,7 +329,7 @@ def _cb_none_concurrency(ctx, param, value):
     "--username", "-u", type=click.STRING, help="Username for HTTP Auth."
 )
 opt_http_password = click.option(
-    "--password", "-p", type=click.STRING, help="Password for HTTP Auth."
+    "--password", type=click.STRING, help="Password for HTTP Auth."
 )
 opt_force = click.option("-f", "--force", is_flag=True, help="Don't ask, just do.")
 opt_src_fs_opts = click.option(

mapchete/commands/_cp.py

@@ -4,7 +4,7 @@
 from multiprocessing import cpu_count
 import os
 from rasterio.crs import CRS
-from shapely.geometry import box, shape
+from shapely.geometry import box, Point, shape
 from shapely.geometry.base import BaseGeometry
 from shapely.ops import unary_union
 from tilematrix import TilePyramid
@@ -28,6 +28,8 @@ def cp(
     area_crs: Union[CRS, str] = None,
     bounds: Tuple[float] = None,
     bounds_crs: Union[CRS, str] = None,
+    point: Tuple[float, float] = None,
+    point_crs: Tuple[float, float] = None,
     overwrite: bool = False,
     workers: int = None,
     multi: int = None,
@@ -59,6 +61,10 @@ def cp(
         Override bounds or area provided in process configuration.
     bounds_crs : CRS or str
         CRS of area (default: process CRS).
+    point : iterable
+        X and y coordinates of point whose corresponding output tile bounds will be used.
+    point_crs : str or CRS
+        CRS of point (defaults to process pyramid CRS).
     overwrite : bool
         Overwrite existing output.
     workers : int
@@ -155,6 +161,8 @@ def _empty_callback(*args):
                     workers,
                     src_fs,
                     dst_fs,
+                    point,
+                    point_crs,
                     overwrite,
                 ),
                 executor_concurrency=concurrency,
@@ -164,7 +172,7 @@ def _empty_callback(*args):
                     dask_client=dask_client,
                 ),
                 as_iterator=as_iterator,
-                total=src_mp.count_tiles(),
+                total=1 if point else src_mp.count_tiles(),
             )
 
 
@@ -176,19 +184,28 @@ def _copy_tiles(
     workers,
     src_fs,
     dst_fs,
+    point,
+    point_crs,
     overwrite,
     executor=None,
 ):
     for z in src_mp.config.init_zoom_levels:
         msg_callback(f"copy tiles for zoom {z}...")
+
         # materialize all tiles
-        aoi_geom = src_mp.config.area_at_zoom(z)
-        tiles = [
-            t
-            for t in tp.tiles_from_geom(aoi_geom, z)
-            # this is required to omit tiles touching the config area
-            if aoi_geom.intersection(t.bbox).area
-        ]
+        if point:
+            point_geom = reproject_geometry(
+                Point(point), src_crs=point_crs or tp.crs, dst_crs=tp.crs
+            )
+            tiles = [tp.tile_from_xy(point_geom.x, point_geom.y, z)]
+        else:
+            aoi_geom = src_mp.config.area_at_zoom(z)
+            tiles = [
+                t
+                for t in tp.tiles_from_geom(aoi_geom, z)
+                # this is required to omit tiles touching the config area
+                if aoi_geom.intersection(t.bbox).area
+            ]
 
         # check which source tiles exist
         logger.debug("looking for existing source tiles...")

test/test_cli.py

@@ -786,12 +786,6 @@ def test_convert_geobuf(landpoly, mp_tmpdir):
                 assert shape(f["geometry"]).area
 
     # convert from geobuf
-    # NOTE: if shapely was built using GEOS 3.8.0 or smaller, there is one more feature
-    if version_is_greater_equal(shapely.geos.geos_version, (3, 9, 0)):
-        zoom9_control = 31
-    else:
-        zoom9_control = 32
-
     geojson_outdir = os.path.join(mp_tmpdir, "geojson")
     run_cli(
         [
@@ -806,12 +800,16 @@ def test_convert_geobuf(landpoly, mp_tmpdir):
             "none",
         ]
     )
-    for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, zoom9_control]):
+    for (zoom, row, col), control in zip([(4, 0, 7), (4, 1, 7)], [9, [31, 32]]):
         out_file = os.path.join(
             *[geojson_outdir, str(zoom), str(row), str(col) + ".geojson"]
         )
         with fiona.open(out_file, "r") as src:
-            assert len(src) == control
+            if isinstance(control, list):
+                assert len(src) in control
+            else:
+                assert len(src) == control
+
             for f in src:
                 assert shape(f["geometry"]).is_valid
 
@@ -1300,6 +1298,21 @@ def test_cp(mp_tmpdir, cleantopo_br, wkt_geom):
     )
     out_path = os.path.join(TESTDATA_DIR, cleantopo_br.dict["output"]["path"])
 
+    # copy tiles and subset by point
+    run_cli(
+        [
+            "cp",
+            out_path,
+            os.path.join(mp_tmpdir, "all"),
+            "-z",
+            "5",
+            "-p",
+            "170",
+            "-85",
+            "--concurrency",
+            "none",
+        ]
+    )
     # copy tiles and subset by bounds
     run_cli(
         [