Add size analysis that is appropriate for geospatial data

[HaleySchuhl]

The traits getting stored with pcv.analyze.size are numerous and most make very little sense in the context of a geospatial workflow. We propose new geospatial analysis tools to answer questions like canopy coverage, stand count, and doing so in a vectorized manner over all regions in a GeoJSON. Example below:

# Imports and setup debug visualization 
from rasterstats import zonal_stats
from plantcv import plantcv as pcv 
import plantcv.geospatial as geo
from plantcv import annotate 
import fiona 
import numpy as np
import os
%matplotlib widget
pcv.params.debug = "plot"

# Read in ortho
tif_filename = "./ortho.tif"  # multispectral data
spectral = geo.read_geotif(filename=tif_filename, bands=[650,560,480,842,717])

# Use RGB for simple segmentation
rgb_img = spectral.pseudo_rgb
gray = pcv.rgb2gray_lab(rgb_img=rgb_img, channel="a")
plant_mask = pcv.threshold.binary(gray_img=gray, threshold=126, object_type="dark") # good for left side of the field

# Geo analysis with Rasterstats package
affine = spectral.metadata["transform"]
# sum gives the sum of pixel values, so change from [0,255] to [0,1] 
b = plant_mask.astype(float) / 255

stats = zonal_stats('./shapefiles/2024-06-06-Baxter_Altum_Left.shp',
                    b,
                    affine=affine,
                   stats="sum")

# Gather IDs from geoJSON & connect to pcv.Outputs
# Grab shapefile ID labels 
img = spectral
geojson = './geojson.shp'
geo_transform = img.metadata["transform"]
ids = [] 
# Gather list of IDs 
with fiona.open(geojson, 'r') as shapefile:
    for row in shapefile:
        temp_list = []
        # Polygon
        if len(row.geometry["coordinates"][0]) > 1:
            square = row.geometry["coordinates"][0][:-1]
            ids.append((row['properties']["ID"]))

for i, id_lbl in enumerate(ids):
    pcv.outputs.add_observation(sample=id_lbl, variable="pixel_count", trait="count", method="rasterstats.zonal_stats", scale="pixels", datatype=int,
                                value=stats[i]["sum"], label="pixels")