# -*- coding: utf-8 -*- """satellite.ipynb Automatically generated by Colaboratory. Original file is located at https://colab.research.google.com/github/opengeos/segment-geospatial/blob/main/docs/examples/satellite.ipynb # Segment Anything Model for Geospatial Data [![image](https://studiolab.sagemaker.aws/studiolab.svg)](https://studiolab.sagemaker.aws/import/github/opengeos/segment-geospatial/blob/main/docs/examples/satellite.ipynb) [![image](https://img.shields.io/badge/Open-Planetary%20Computer-black?style=flat&logo=microsoft)](https://pccompute.westeurope.cloudapp.azure.com/compute/hub/user-redirect/git-pull?repo=https://github.com/opengeos/segment-geospatial&urlpath=lab/tree/segment-geospatial/docs/examples/satellite.ipynb&branch=main) [![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/segment-geospatial/blob/main/docs/examples/satellite.ipynb) This notebook shows how to use segment satellite imagery using the Segment Anything Model (SAM) with a few lines of code. Make sure you use GPU runtime for this notebook. For Google Colab, go to `Runtime` -> `Change runtime type` and select `GPU` as the hardware accelerator. ## Install dependencies Uncomment and run the following cell to install the required dependencies. """ # %pip install segment-geospatial """## Import libraries""" from osgeo import gdal from PIL import Image import numpy as np from socket import * from oct2py import octave from io import BytesIO import traceback import os import leafmap from samgeo import SamGeo, tms_to_geotiff, get_basemaps #from samgeo.hq_sam import SamGeo, tms_to_geotiff, get_basemaps import time from shapely.geometry import Point, Polygon import json import pyproj import requests import cv2 as cv """## Create an interactive map""" start_time = time.time() def make_dir(uid): if not os.path.exists(uid): os.makedirs(uid) def convert_coordinates_to_epsg4326(coordinates, source_epsg): # Define the source and target CRS source_crs = pyproj.CRS.from_epsg(source_epsg) target_crs = pyproj.CRS.from_epsg(4326) # Create a PyProj Transformer transformer = pyproj.Transformer.from_crs(source_crs, target_crs, always_xy=True) # Convert the coordinates converted_coordinates = [] for coord_pair in coordinates: x, y = coord_pair lon, lat = transformer.transform(x, y) converted_coordinates.append([lon, lat]) return converted_coordinates def submit_farm(uid, coords, field_name): obj = { "UID": uid, "Points":coords, "PaymentType":3, "CropCode":999, "FieldName":field_name } url = 'https://us-central1-farmbase-b2f7e.cloudfunctions.net/submitField' json_data = json.dumps(obj) # Set the headers to specify that you are sending JSON data headers = {'Content-Type': 'application/json'} # Send the POST request response = requests.post(url, data=json_data, headers=headers) # Check the response if response.status_code == 200: print("POST request was successful.") print("Response content:") print(response.text) else: print(f"POST request failed with status code {response.status_code}") def create_geojson(coordinates, geometry_type="Point", properties=None): """ Create a GeoJSON object from a coordinates array. Args: coordinates (list): List of coordinates in the form [longitude, latitude]. geometry_type (str): Type of geometry (e.g., "Point", "LineString", "Polygon"). properties (dict): Optional properties to include in the GeoJSON properties field. Returns: str: GeoJSON string. """ if geometry_type not in ["Point", "LineString", "Polygon"]: raise ValueError("Invalid geometry type. Supported types are 'Point', 'LineString', and 'Polygon'.") # Create a GeoJSON feature feature = { "type": "Feature", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::4326" } }, "geometry": { "type": geometry_type, "coordinates": [coordinates] } } # Add properties if provided if properties is not None: feature["properties"] = properties # Create a GeoJSON object geojson = { "type": "FeatureCollection", "features": [feature] } # Serialize the GeoJSON object to a JSON string geojson_string = json.dumps(geojson) return geojson_string def find_polygon(point, uid, field_name): # Replace 'your_geojson_file.geojson' with the path to your GeoJSON file file_name = uid + "/" + field_name + "segment.geojson" with open(file_name, 'r') as geojson_file: geojson_data = json.load(geojson_file) # Access the type of the GeoJSON object (e.g., "FeatureCollection", "Feature", etc.) geojson_type = geojson_data['type'] # print(f"GeoJSON Type: {geojson_type}") # Access features in a FeatureCollection (if applicable) if geojson_type == 'FeatureCollection': features = geojson_data['features'] for feature in features: coordinates = feature['geometry']['coordinates'] source_epsg = 3857 # Example source EPSG code (UTM Zone 33N) converted_coordinates = convert_coordinates_to_epsg4326(coordinates[0], source_epsg) #print(converted_coordinates) if is_point_inside_polygon(point, converted_coordinates) == True: new_geojson = create_geojson(converted_coordinates, "Polygon",None) # print(new_geojson) return new_geojson, converted_coordinates def geotiff_to_png(geotiff_file, png_file): # Open the GeoTIFF file dataset = gdal.Open(geotiff_file) if dataset is None: raise ValueError("Unable to open GeoTIFF file.") # Get the number of bands in the GeoTIFF num_bands = dataset.RasterCount # Create an empty list to store the band data band_data = [] # Read each band and append it to the list for i in range(1, num_bands + 1): band = dataset.GetRasterBand(i) data = band.ReadAsArray() band_data.append(data) # Create a PNG image from the band data image = Image.fromarray(np.array(band_data).transpose(1, 2, 0)) # Save the PNG image image.save(png_file, 'PNG') octave.makeETCI2(png_file, png_file) dataset = None def png_to_geotiff(png_file, geotiff_file, original_geotiff_file): # Open the original GeoTIFF file to copy the geospatial metadata original_dataset = gdal.Open(original_geotiff_file) if original_dataset is None: raise ValueError("Unable to open original GeoTIFF file.") # Get the number of bands in the original GeoTIFF num_bands = original_dataset.RasterCount # Create a new GeoTIFF file with the same geospatial metadata driver = gdal.GetDriverByName('GTiff') dataset = driver.Create( geotiff_file, original_dataset.RasterXSize, original_dataset.RasterYSize, num_bands, # Number of bands gdal.GDT_Byte # Data type ) dataset.SetProjection(original_dataset.GetProjection()) dataset.SetGeoTransform(original_dataset.GetGeoTransform()) # Open the PNG image and write each band to the new GeoTIFF file png_image = Image.open(png_file) png_data = np.array(png_image) for i in range(num_bands): band = dataset.GetRasterBand(i + 1) band.WriteArray(png_data[:, :, i]) # Close the new GeoTIFF and original GeoTIFF datasets dataset = None original_dataset = None # Example usage: # Convert GeoTIFF to PNG # Convert PNG back to GeoTIFF def gray_img(uid, field_name): input_file = uid + "/" + field_name + "satellite.tif" output_file = uid + "/" + field_name + "satellite.png" geotiff_to_png(input_file, output_file) png_to_geotiff(output_file, input_file, input_file) def is_point_inside_polygon(point, polygon_coordinates): # Create a Shapely Point object from the point coordinates point = Point(point) # Create a Shapely Polygon object from the polygon coordinates polygon = Polygon(polygon_coordinates) # Check if the point is inside the polygon return point.within(polygon) def create_rectangle_coords(bbox): min_lng, min_lat, max_lng, max_lat = bbox rectangle_coords = [ [min_lng, min_lat], [min_lng, max_lat], [max_lng, max_lat], [max_lng, min_lat], [min_lng, min_lat] # Closing the rectangle by returning to the starting point ] return rectangle_coords def simplify_coordinates(coordinates, tolerance): coordinates = np.array(coordinates) # Convert to NumPy array if len(coordinates) <= 2: return coordinates.tolist() # Convert back to list # Find the point with the maximum distance max_distance = 0 index = 0 end = len(coordinates) - 1 for i in range(1, end): distance = perpendicular_distance(coordinates[i], coordinates[0], coordinates[end]) if distance > max_distance: max_distance = distance index = i # If the maximum distance is greater than the tolerance, recursively simplify if max_distance > tolerance: simplified1 = simplify_coordinates(coordinates[:index+1], tolerance) simplified2 = simplify_coordinates(coordinates[index:], tolerance) return simplified1[:-1] + simplified2 # No need to convert back to list else: return [coordinates[0], coordinates[end]] # No need to convert back to list def perpendicular_distance(point, start, end): if np.all(start == end): return np.linalg.norm(np.array(point) - np.array(start)) line_length = np.linalg.norm(np.array(end) - np.array(start)) t = np.dot(np.array(point) - np.array(start), np.array(end) - np.array(start)) / (line_length * line_length) t = np.clip(t, 0, 1) closest_point = np.array(start) + t * (np.array(end) - np.array(start)) return np.linalg.norm(np.array(point) - closest_point) def get_main_segment(lat,lng,uid, field_name): make_dir(uid) m = leafmap.Map(center=[lat,lng], zoom=19) m.add_basemap("SATELLITE") m coord_span = 0.001 min_lng, min_lat, max_lng, max_lat = (lng-coord_span),(lat-coord_span),(lng+coord_span),(lat+coord_span) bbox = [min_lng,min_lat,max_lng,max_lat] boundary_field = create_rectangle_coords(bbox) image = uid + "/" + field_name +"satellite.tif" tms_to_geotiff(output=image, bbox=bbox, zoom=18, source="Satellite", overwrite=True) #gray_img(uid, field_name) sam_kwargs = { "points_per_side": 32, "pred_iou_thresh": 0.9, "stability_score_thresh": 0.9, "crop_n_layers":1, "crop_n_points_downscale_factor": 1, "min_mask_region_area": 200, } sam = SamGeo( model_type="vit_h", checkpoint="sam_vit_h_4b8939.pth", sam_kwargs=sam_kwargs ) mask = uid + "/" + field_name + "segment.tif" sam.generate( image, mask, batch=True, foreground=True, erosion_kernel=(3,3), mask_multiplier=255 ) kmlfile = uid + "/" + field_name + "segment.geojson" sam.tiff_to_geojson(mask, kmlfile) try: polygon_geojson, found_coordinates = find_polygon([lng, lat], uid, field_name) # Adjust this value to control the level of simplification tolerance = 0.00005 found_coordinates = simplify_coordinates(found_coordinates, tolerance) print(len(found_coordinates)) print(found_coordinates) print(boundary_field) if len(found_coordinates) < 4: found_coordinates = boundary_field new_geojson_name = uid + "/" + field_name + ".geojson" #print(polygon_geojson) with open(new_geojson_name, "w") as outfile: outfile.write(polygon_geojson) #submit_farm(uid, [found_coordinates], field_name) return found_coordinates #print(found_coordinates) except: found_coordinates = boundary_field print(traceback.format_exc()) return found_coordinates # fields_data = [[26.382823,91.16466,'SA_PAK_ERA_01'], # [26.384481,91.14941,'SA_PAK_ERG_020'], # [26.38304,91.148953,'SA_PAK_ERG_033'], # [26.383156,91.151065,'SA_PAK_ERG_025'], # [26.38179,91.150906,'SA_PAK_ERG_019'], # [26.381722,91.165892,'SA_PAK_ERA_024'], # [26.379582,91.165988,'SA_PAK_ERA_006'], # [26.339145,91.482569,'SA_BAR_RAI_026'], # [26.340619,91.484231,'SA_BAR_RAI_001'], # [26.525938,92.228413,'MK_DAL_BAD_015'], # [26.526865,92.226614,'MK_DAL_BAD_005'], # [26.527791,92.226911,'MK_DAL_BAD_006'], # [26.528246,92.229085,'MK_DAL_BAD_010'], # [26.529123,92.227946,'MK_DAL_BAD_020'] # ] # fields_data = [[22.066756000000002,87.514644000000004, 'AS-PAT-I-GOK-016'], # [22.067325000000000,87.512485999999996,'AS-PAT-I-GOK-095'], # [22.066718999999999,87.511988000000002,'AS-PAT--GOK-038'], # [22.067743000000000,87.511157999999995,'AS-PAT-I-GOK-086'], # [25.343056000000001,87.859827999999993,'SB-HAR2-TAL-001'], # [26.176372000000001,89.176372000000001,'SS-DIN-1-BHA-091'], # [26.317613999999999,89.290958000000003,'SS-MTB2-SB-031'], # [26.167496000000000,89.431129999999996,'SS-DIN-1-BHA-075'], # [26.169864000000000,89.435573000000005,'SS-DIN-1-BHA-033'], [26.171410999999999,89.436895000000007,'SS-DIN-1-BHA-007'], # [26.170113000000001,89.437944000000002,'SS-DIN-1-BHA-041'], [26.175407000000000,89.435907999999998,'SS-DIN-1BHA-013'], # [26.172504000000000,89.437334000000007,'SS-DIN-1-BHA-049'], # [26.171040000000001,89.438339999999997,'SS-DIN-1-BHA-079'], [26.169896000000001,89.439606999999995,'SS-DIN-1-BHA-040']] # #fields_data = [[16.1336841,75.0937387,'MELLIKERI_1'],[16.1335475,75.0877041, 'MELLIKERI_2'],[16.1431067,75.0760271,'MELLIKERI_3'],[16.1275247,75.063244,'MELLIKERI_4'],[16.0928061,75.3923695,'BATAKURKI_1']] # uid = 'F34CLcmyzzczz7gt764WRQKXda72' # #uid = 'R4pbfLc1FRaEDRPMHA8uc1OnexH2' # #octave.makeETCI2((uid + "/SA_PAK_ERG_020satellite.tif"),(uid + "/SA_PAK_ERG_020satellite2.tif")) # for field in fields_data: # print(field) # start_time = time.time() # get_main_segment(field[0], field[1],uid, field[2]) # total_time = time.time() - start_time # print(total_time) #main(12.968766, 77.791663,'BpkwnSJdwHTjKhdm8ZWKJBO1HUn2', 'sample11')