Create plot_extraction.md

README.md

@@ -11,5 +11,10 @@ Developers should set up a PlantCV conda environment from source code as normal,
 
 ```bash
 cd plantcv-geospatial
-pip install -e .
+<<<<<<< Updated upstream
+pip install -e . --config-settings editable_mode=strict
+
+=======
+pip install -e . 
+>>>>>>> Stashed changes
 ```

docs/plot_extraction.md

@@ -0,0 +1,34 @@
+## Create a rectangular shapefile to extract area of interest from agricultural research plots
+
+
+**plantcv.geospatial.extract_plot**(*input_filename, output_filename, ranges, columns, row_length, vertical_alley, horizontal_alley)
+
+**returns** [Multi-Polygon Shapefile]
+
+- **Parameters:**
+    - input_filename - Filepath to input shapefile consiting four corners of the field
+    - output_filename - Filepath to output shapefile containing individual rectangular ploygon corresponding to specific research plot
+    - ranges - Number of ranges
+    - columns - Number of columns
+    - row_length - Row longth along column in meters
+    - vertical_alley - Length of alley between ranges in meters
+    - horziontal_alley - Length of alley between columns in meters
+
+- **Example use:**
+    - below
+
+
+```python
+import plantcv.plantcv as pcv 
+import plantcv.geospatial.extract_plot as gridcell
+
+# Path of input and output shapefile
+input_shapefile_path = '712_corner_points.shp'
+output_shapefile_path = '712_polygon.shp'  
+output_shapefile = gridcell.write_shapefile(input_shapefile_path, output_shapefile_path, number_of_ranges=45, number_of_columns=12, 
+row_length_along_column=2.4384, vertical_alley=0.9144, 
+horizontal_alley=0.0)
+
+```
+
+**Source Code:** [Here](https://github.com/danforthcenter/plantcv-geospatial/blob/main/plantcv/geospatial/plot_extract.py)

mkdocs.yml

@@ -18,6 +18,7 @@ nav:
   - 'PlantCV Namespace':
     - 'Geopspatial Tools':
         - Read Geo-tif Data: read_geotif.md
+        - Extract Plot : plot_extraction.md
 markdown_extensions:
   - toc:
       permalink: True

plantcv/geospatial/plot_extraction.py

@@ -0,0 +1,82 @@
+import math
+import fiona
+from shapely.geometry import mapping, Polygon, LineString
+
+
+def create_rectangle_from_points(shapefile_path):
+    with fiona.open(shapefile_path, 'r') as shapefile:
+        points = [shape['geometry']['coordinates'] for shape in shapefile]
+    rectangle = Polygon(points)
+    return rectangle
+
+
+def calculate_line_length_in_meters(line):
+    total_length = 0
+    for i in range(len(line.coords)-1):
+        x1, y1 = line.coords[i]
+        x2, y2 = line.coords[i+1]
+        segment_length = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
+        total_length += segment_length
+    return total_length
+
+
+def get_field_edges(rectangle):
+    coords = list(rectangle.exterior.coords)
+    edge1 = LineString([coords[0], coords[1]])
+    edge2 = LineString([coords[0], coords[3]])
+    return edge1, edge2
+
+
+def create_grid(rectangle, edge1, edge2, columns, row_length_along_column, vertical_gap):
+    edge1_length = calculate_line_length_in_meters(edge1)
+    edge2_length = calculate_line_length_in_meters(edge2)
+    horizontal_threshold = edge1_length/columns
+    vertical_threshold = row_length_along_column
+    divisions_along_edge1 = round(edge1_length / horizontal_threshold)
+    divisions_along_edge2 = round((edge2_length - vertical_gap) / (vertical_threshold + vertical_gap)) + 1
+    edge1_start = edge1.coords[0]
+    edge2_start = edge2.coords[0]
+
+    edge1_dir = ((edge1.coords[1][0] - edge1_start[0]) / edge1.length,
+                 (edge1.coords[1][1] - edge1_start[1]) / edge1.length)
+    edge2_dir = ((edge2.coords[1][0] - edge2_start[0]) / edge2.length,
+                 (edge2.coords[1][1] - edge2_start[1]) / edge2.length)
+
+    grid_cells = []
+
+    for column_number in range(divisions_along_edge1):
+        for row_number in range(divisions_along_edge2):
+            start_x = edge1_start[0] + column_number * horizontal_threshold * edge1_dir[0] + row_number * (vertical_threshold + vertical_gap) * edge2_dir[0]
+            start_y = edge1_start[1] + column_number * horizontal_threshold * edge1_dir[1] + row_number * (vertical_threshold + vertical_gap) * edge2_dir[1]
+            bottom_left = (start_x, start_y)
+            bottom_right = (bottom_left[0] + horizontal_threshold * edge1_dir[0],
+                            bottom_left[1] + horizontal_threshold * edge1_dir[1])
+            top_left = (bottom_left[0] + vertical_threshold * edge2_dir[0],
+                        bottom_left[1] + vertical_threshold * edge2_dir[1])
+            top_right = (bottom_right[0] + vertical_threshold * edge2_dir[0],
+                         bottom_right[1] + vertical_threshold * edge2_dir[1])
+
+            cell = Polygon([bottom_left, bottom_right, top_right, top_left, bottom_left])
+            if cell.intersects(rectangle):
+                grid_cells.append({"polygon": cell, "row": row_number + 1, "column": column_number + 1})
+
+    return grid_cells
+
+
+def write_shapefile(number_of_columns, row_length_along_column, output_shapefile_path, input_shapefile_path, vertical_gap):
+    rectangle = create_rectangle_from_points(input_shapefile_path)
+    edge1, edge2 = get_field_edges(rectangle)
+    grid_cells = create_grid(rectangle, edge1, edge2, number_of_columns, row_length_along_column, vertical_gap)
+
+    with fiona.open(input_shapefile_path, 'r') as input_shapefile:
+        driver = input_shapefile.driver
+        crs = input_shapefile.crs
+        schema = {'geometry': 'Polygon', 'properties': {'id': 'int', 'row': 'int', 'column': 'int', 'label': 'str'}}
+
+    with fiona.open(output_shapefile_path, 'w', driver=driver, crs=crs, schema=schema) as output_shapefile:
+        for idx, cell in enumerate(grid_cells):
+            label = f"({cell['row']},{cell['column']})"
+            output_shapefile.write({
+                'geometry': mapping(cell["polygon"]),
+                'properties': {'id': idx, 'row': cell["row"], 'column': cell["column"], 'label': label},
+            })

pyproject.toml

@@ -12,6 +12,9 @@ dependencies = [
     "plantcv",
     "matplotlib",
     "rasterio",
+    "pyshp",
+    "shapely",
+    "fiona",
 
 ]
 requires-python = ">=3.6"

tests/testdata/0615shapefile.cpg

@@ -0,0 +1 @@
+UTF-8

tests/testdata/0615shapefile.prj

@@ -0,0 +1 @@
+PROJCS["WGS_1984_UTM_Zone_15N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-93.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]