import rasterio
import numpy as np
import cv2
from shapely.geometry import Polygon, mapping
import geojson

def read_geotiff(file_path):
    with rasterio.open(file_path) as src:
        img = src.read([1, 2, 3])  # Read RGB channels
        img = np.stack(img, axis=2)  # Stack channels to create an RGB image
        transform = src.transform  # Get the transformation matrix
        crs = src.crs  # Get the coordinate reference system
    return img, transform, crs

def detect_oil_palm_plantations(img):
    # Convert the image to HSV color space
    hsv_img = cv2.cvtColor(img, cv2.COLOR_RGB2HSV)
    
    # Define the range for green color in HSV space
    lower_green = np.array([30, 40, 40])
    upper_green = np.array([80, 255, 255])
    
    # Create a binary mask where green colors are white and the rest are black
    mask = cv2.inRange(hsv_img, lower_green, upper_green)
    
    # Find contours in the mask
    contours, _ = cv2.findContours(mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    
    # Approximate contours to polygons
    polygons = [cv2.approxPolyDP(contour, 0.01 * cv2.arcLength(contour, True), True) for contour in contours]
    
    # Convert contours to polygons
    polygons = [Polygon(polygon[:, 0, :]) for polygon in polygons if len(polygon) > 2]
    
    return polygons

def save_geojson(polygons, transform, crs, output_path):
    features = []
    for polygon in polygons:
        # Transform polygon coordinates from image to geographical coordinates
        transformed_polygon = [rasterio.transform.xy(transform, point[1], point[0]) for point in polygon.exterior.coords]
        geojson_polygon = Polygon(transformed_polygon)
        features.append(geojson.Feature(geometry=geojson_polygon))
    
    feature_collection = geojson.FeatureCollection(features)
    
    with open(output_path, 'w') as f:
        geojson.dump(feature_collection, f)

def detect_and_save_oil_palm_plantations(file_path, output_path):
    img, transform, crs = read_geotiff(file_path)
    polygons = detect_oil_palm_plantations(img)
    save_geojson(polygons, transform, crs, output_path)

# Example usage:
file_path = 'Downloads/satellite_new8.tif'
output_path = 'plantations.geojson'
detect_and_save_oil_palm_plantations(file_path, output_path)