# -*- coding: utf-8 -*- """Mapping_palm_trees_in eluru_clean.ipynb Automatically generated by Colab. For reference contact - Honey Jain. """ # Commented out IPython magic to ensure Python compatibility. # %pip install leafmap # %pip install ultralytics # %pip install alphashape # %pip install hdbscan # %pip install rasterio # %pip install concave_hull ###################################GENERATING BOUNDING BOXES COVERING ELURU################################################ import matplotlib.pyplot as plt from shapely.geometry import Polygon, box from matplotlib.patches import Polygon as mpl_polygon import json # file_path = '/content/drive/MyDrive/Eluru/dict_grid_bboxes.json' def extract_bounding_boxes_from_city_coordinates(eluru_coordinates_ori, file_path_json): eluru_coordinates = [] sorted_keys = sorted(eluru_coordinates_ori.keys(), key=lambda x: (x == 'a', int(x.split('_')[1]) if '_' in x else 0)) eluru_coordinates_sorted = {key: eluru_coordinates_ori[key] for key in sorted_keys} for point, coord_data in eluru_coordinates_sorted.items(): longitude = coord_data.get('Longitude', None) latitude = coord_data.get('Latitude', None) if longitude is not None and latitude is not None: eluru_coordinates.append((longitude, latitude)) # Ensure the polygon is closed by adding the first point again at the end if len(eluru_coordinates) > 0: eluru_coordinates.append(eluru_coordinates[0]) eluru_polygon_coords = eluru_coordinates # Define the size of each grid cell in pixels cell_size_pixels = 24 # Pixel size # Convert pixel size to degrees (latitude and longitude) # This conversion depends on the region and latitude where the polygon is located # Example: In India, 1 degree of latitude ~= 111 km, 1 degree of longitude ~= 111 km * cos(latitude) # Adjust this conversion factor based on the specific location # For simplicity, assuming a conversion factor based on latitude (approximation) conversion_factor = 0.0001 # Adjust as per your region and specific needs # Convert pixel size to degrees (latitude and longitude) cell_size_latlon = cell_size_pixels * conversion_factor # Create a Polygon object for the Eluru district eluru_polygon = Polygon(eluru_polygon_coords) # Determine bounding box of the entire Eluru polygon min_x, min_y, max_x, max_y = eluru_polygon.bounds # Calculate the range of latitude and longitude range_lat = max_y - min_y range_lon = max_x - min_x # Calculate the number of divisions needed along latitude and longitude axes num_divisions_lat = int(range_lat / cell_size_latlon) + 1 num_divisions_lon = int(range_lon / cell_size_latlon) + 1 # Initialize an empty dictionary to store bounding boxes bounding_boxes = {} # Initialize figure and axis for plotting #fig, ax = plt.subplots(figsize=(10, 10)) # Plot the Eluru polygon x, y = eluru_polygon.exterior.xy # ax.plot(x, y, color='blue', alpha=0.7, linewidth=2, solid_capstyle='round', zorder=2) # Iterate over each grid cell for i in range(num_divisions_lon): for j in range(num_divisions_lat): # Calculate bounding box for current cell in latitude and longitude cell_min_lon = min_x + i * cell_size_latlon cell_max_lon = min_x + (i + 1) * cell_size_latlon cell_min_lat = min_y + j * cell_size_latlon cell_max_lat = min_y + (j + 1) * cell_size_latlon # Create a box Polygon for the current cell in lat/lon cell_bbox = box(cell_min_lon, cell_min_lat, cell_max_lon, cell_max_lat) # Intersect the Eluru polygon with the current cell bounding box intersected_polygon = eluru_polygon.intersection(cell_bbox) # Check if intersection resulted in a valid polygon if not intersected_polygon.is_empty: # Store bounding box coordinates in dictionary bounding_boxes[f"{i}_{j}"] = intersected_polygon.bounds # Plot grid cell as a polygon grid_polygon = mpl_polygon([(cell_min_lon, cell_min_lat), (cell_max_lon, cell_min_lat), (cell_max_lon, cell_max_lat), (cell_min_lon, cell_max_lat)], edgecolor='gray', linewidth=0.5, facecolor='none', zorder=1) #ax.add_patch(grid_polygon) # Set plot title and labels # ax.set_title('Eluru District Polygon with Grid Cells') # ax.set_xlabel('Longitude') # ax.set_ylabel('Latitude') # Show plot # plt.tight_layout() # plt.show() with open(file_path, 'w') as file: json.dump(bounding_boxes, file) file.close() print(f"{len(bounding_boxes.keys())} bounding boxes were generated to extract images of Eluru District.") print(f"These bounding boxes are saved at {file_path}.") return bounding_boxes ######################################################################################################################## import leafmap from samgeo import SamGeo, tms_to_geotiff from shapely.geometry import Polygon def pull_image_leafmap(polygonNo, min_lat, max_lat, min_lng, max_lng, output_folder): # Expand the bounding box slightly to ensure coverage buffer = 0.00000 min_lng -= buffer min_lat -= buffer max_lng += buffer max_lat += buffer bbox=[min_lng, min_lat, max_lng, max_lat] image = f"{output_folder}{polygonNo}.png" with open(os.devnull, 'w') as devnull: tms_to_geotiff(output=image, bbox=[min_lng, min_lat, max_lng, max_lat], zoom=20, source="Satellite", overwrite=True) import rasterio from pyproj import Transformer ########################################################################################################################### def save_pixel_as_lat_lon(image_path, row, col): # Open the .ong file using rasterio with rasterio.open(image_path) as dataset: # Read the affine transformation and CRS (coordinate reference system) transform = dataset.transform crs = dataset.crs # Transformer to convert coordinates to WGS84 (lat/lon) transformer = Transformer.from_crs(crs, 'EPSG:4326', always_xy=True) x, y = transform * (col, row) lon, lat = transformer.transform(x, y) return lon, lat def pixel_to_geo_coordinates(pixel_coords, ref_lat, ref_lon, lat_deg_per_pixel, lon_deg_per_pixel): # Convert pixel coordinates to degrees delta_x = pixel_coords[:, 0] * lon_deg_per_pixel delta_y = pixel_coords[:, 1] * lat_deg_per_pixel # Calculate latitude and longitude latitudes = ref_lat + delta_y longitudes = ref_lon + delta_x return np.column_stack((longitudes, latitudes)) import geojson import json def save_as_geojson(polygon_coordinates_geo, output_file): features = [] for polygon_coords in polygon_coordinates_geo: # Create a GeoJSON Polygon feature polygon_feature = geojson.Feature(geometry=geojson.Polygon([polygon_coords])) features.append(polygon_feature) # Create a Feature Collection feature_collection = geojson.FeatureCollection(features) # Write to a GeoJSON file with open(output_file, 'w') as f: json.dump(feature_collection, f) from PIL import Image def get_image_dimensions(image_path): with Image.open(image_path) as img: width, height = img.size return width, height ########################################################################################################################### ############################DETECTING PALM TREES USING YOLOv8 pretrained model############################################# #NOTE: Change the paths for files #path_model_file- contains the model file with weights yolov8_palm_detection_best.pt import os import cv2 from torchvision.models.detection import fasterrcnn_resnet50_fpn from ultralytics import YOLO #import hdbscan #import alphashape import numpy as np from shapely.geometry import Polygon, MultiPolygon, Point, mapping def detect_palm_trees(path_model_file, image_path, confidence_thres=0.15, iou_thres=0.20): # Load the trained model model = YOLO(path_model_file) tree_count = {} matched_locations = [] frame = cv2.imread(image_path) height, width, channel = frame.shape model_imgsz = max(height, width) results = model(frame, conf=confidence_thres, imgsz=model_imgsz, iou=iou_thres, max_det=10000, save=True, save_txt=True, save_json=True)[0] detect_count = len(results.boxes) for result in results.boxes.data.tolist(): x1, y1, x2, y2, score, class_id = result center_x = int((x1 + x2) / 2) center_y = int((y1 + y2) / 2) matched_locations.append((center_x, center_y)) return matched_locations, detect_count ########################################################################################################################## #################################### CLUSTERING AND MAKING POLYGONS AROUND CLUSTERS####################################### import matplotlib.pyplot as plt import numpy as np from scipy.spatial import ConvexHull from concave_hull import concave_hull, concave_hull_indexes import cv2 from sklearn.cluster import DBSCAN from shapely.geometry import MultiPoint, Polygon def make_concave_polygon(points): points = np.array(points) convex_hull = ConvexHull(points[:, :2]) # it's already N-by-2, I'm just emphasizing idxes = concave_hull_indexes( points[:, :2], length_threshold=50, ) # Ensure the polygon is closed by connecting the last point back to the first point idxes = np.append(idxes, idxes[0]) # Append the first index to close the polygon # Plot the concave hull polygon_points = points[idxes] # Extract the points forming the polygon #plt.plot(polygon_points[:, 0], polygon_points[:, 1], "r-", alpha=1) #plt.show() return polygon_points def palm_tree_polygon_detection(image_path, matched_locations, eps_def, min_samples_def, alpha): matched_locations = np.array(matched_locations) image = cv2.imread(image_path) # Apply DBSCAN clustering db = DBSCAN(eps=eps_def, min_samples=min_samples_def).fit(matched_locations) # Adjust eps and min_samples as needed labels = db.labels_ # Number of clusters (excluding noise) num_clusters = len(set(labels)) - (1 if -1 in labels else 0) # Create polygons around each cluster clusters = [matched_locations[labels == i] for i in range(num_clusters)] polygons = [] for cluster in clusters: cluster_array = np.array(cluster) cluster_tuples = [tuple(point) for point in cluster_array] #print("cluster", cluster_tuples) if len(cluster_tuples) > 3: polygon = make_concave_polygon(cluster_tuples) polygon = Polygon(polygon) if polygon is not None and not polygon.is_empty: if polygon.geom_type == 'Polygon': polygons.append(polygon) elif polygon.geom_type == 'MultiPolygon': polygons.extend(polygon.geoms) for matched_location in matched_locations: center_x, center_y = matched_location # Draw circle at the center of the bounding box # cv2.circle(image, (center_x, center_y), radius=5, color=(0, 255, 255), thickness=-1) # -1 thickness fills the circle # #Draw polygons on the image # for polygon in polygons: # if polygon is not None and not polygon.is_empty: # exterior = np.array(polygon.exterior.coords).astype(int) # cv2.polylines(image, [exterior], isClosed=True, color=(0, 255, 0), thickness=3) # Save and display the result #cv2.imwrite(f"{output_path}", image) #plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB)) #plt.show() # List to store polygon coordinates polygon_coordinates = [] for polygon in polygons: if polygon is not None and not polygon.is_empty: exterior_coords = np.array(polygon.exterior.coords) exterior_coords = np.round(exterior_coords).astype(np.int32) # Round and convert to int32 polygon_coordinates.append(exterior_coords) return polygon_coordinates ############################################################################################################################################### #######################################################CHANGE VARIABLES AS REQUIRED############################################################ # # # eluru_coordinates_ori = {'P_25': {'Latitude': 16.840212019313366, 'Longitude': 81.32598024542183}, 'P_23': {'Latitude': 16.790507310842816, 'Longitude': 81.3778960895139}, 'P_4': {'Longitude': 81.08118721355505, 'Latitude': 16.572988033089736}, 'P_58': {'Longitude': 80.97052550336397, 'Latitude': 16.71005962580638}, 'P_63': {'Longitude': 80.96365904828585, 'Latitude': 16.66072951834851}, 'P_64': {'Latitude': 16.66641839668976, 'Longitude': 80.94248482097693}, 'P_7': {'Latitude': 16.59141451932586, 'Longitude': 81.10109993328162}, 'P_60': {'Longitude': 81.01172423383272, 'Latitude': 16.660071697620122}, 'P_17': {'Longitude': 81.29961850162327, 'Latitude': 16.68661417971704}, 'P_34': {'Longitude': 81.34863954717964, 'Latitude': 17.021736585095542}, 'P_53': {'Latitude': 16.809338874377612, 'Longitude': 81.02820372602022}, 'P_38': {'Longitude': 81.29027467901558, 'Latitude': 16.991532106237308}, 'P_30': {'Latitude': 16.966576909463615, 'Longitude': 81.30057436163277}, 'P_31': {'Latitude': 16.982338472506378, 'Longitude': 81.32872682745308}, 'P_52': {'Longitude': 80.99867796918429, 'Latitude': 16.825113580740112}, 'P_33': {'Longitude': 81.3678656213984, 'Latitude': 16.99415877594965}, 'P_16': {'Longitude': 81.26597287174046, 'Latitude': 16.650435575989544}, 'P_13': {'Longitude': 81.22076526064241, 'Latitude': 16.5623053172597}, 'P_14': {'Latitude': 16.57941663401863, 'Longitude': 81.23243823427522}, 'P_42': {'Latitude': 16.92059813141939, 'Longitude': 81.2250433557734}, 'P_1': {'Longitude': 80.9356183658988, 'Latitude': 16.64405188875733}, 'P_3': {'Longitude': 80.9960431705863, 'Latitude': 16.54666142219228}, 'P_36': {'Latitude': 17.034210926210108, 'Longitude': 81.29782777960152}, 'P_29': {'Latitude': 16.948843569751762, 'Longitude': 81.3239203088984}, 'P_18': {'Latitude': 16.669512419515833, 'Longitude': 81.31815793033421}, 'P_46': {'Latitude': 16.877236428245077, 'Longitude': 81.20169740850777}, 'P_59': {'Longitude': 80.98151183148897, 'Latitude': 16.669280982108884}, 'P_54': {'Longitude': 81.03163695355929, 'Latitude': 16.790276020296464}, 'P_61': {'Latitude': 16.63836234611751, 'Longitude': 81.00897765180147}, 'P_49': {'Longitude': 81.08702760870308, 'Latitude': 16.902203475016783}, 'P_22': {'Longitude': 81.36759640689671, 'Latitude': 16.756978380714404}, 'P_57': {'Latitude': 16.710717274533256, 'Longitude': 80.99455809613741}, 'P_51': {'Latitude': 16.85620898610416, 'Longitude': 80.9943304651484}, 'P_11': {'Latitude': 16.613634702231362, 'Longitude': 81.1816264666971}, 'P_56': {'Latitude': 16.735048683614252, 'Longitude': 81.00623106977022}, 'P_28': {'Latitude': 16.948843569751762, 'Longitude': 81.30950075323433}, 'P_39': {'Longitude': 81.27310854132027, 'Latitude': 16.995472097005653}, 'P_24': {'Latitude': 16.84284079979227, 'Longitude': 81.36717897589058}, 'P_55': {'Longitude': 81.00279784223116, 'Latitude': 16.76003446479851}, 'P_44': {'Longitude': 81.21749025518746, 'Latitude': 16.873293961084812}, 'P_6': {'Latitude': 16.610497233144432, 'Longitude': 81.07912727703162}, 'P_62': {'Longitude': 80.98906493207491, 'Latitude': 16.62980950218376}, 'P_37': {'Longitude': 81.29988771612496, 'Latitude': 17.011231231825114}, 'P_10': {'Latitude': 16.60918124475655, 'Longitude': 81.15603157390662}, 'P_26': {'Longitude': 81.32323366339058, 'Latitude': 16.918402603910245}, 'P_48': {'Longitude': 81.13303285772652, 'Latitude': 16.895633519851334}, 'a': {'Longitude': 80.94248482097693, 'Latitude': 16.66641839668976}, 'P_8': {'Latitude': 16.60194314747709, 'Longitude': 81.1223859440238}, 'P_32': {'Latitude': 16.970517424293444, 'Longitude': 81.34451967413277}, 'P_50': {'Latitude': 16.862123190672566, 'Longitude': 81.05887514288277}, 'P_45': {'Latitude': 16.873293961084812, 'Longitude': 81.21749025518746}, 'P_9': {'Latitude': 16.580885314689894, 'Longitude': 81.14092537273474}, 'P_20': {'Latitude': 16.73396491100082, 'Longitude': 81.32983090396702}, 'P_41': {'Latitude': 17.01779714670282, 'Longitude': 81.24976259405464}, 'P_40': {'Latitude': 17.013857625405922, 'Longitude': 81.27242189581246}, 'P_43': {'Latitude': 16.889063335852544, 'Longitude': 81.23877626592964}, 'P_19': {'Longitude': 81.34699704166233, 'Latitude': 16.685956447991973}, 'P_21': {'Latitude': 16.75895083439045, 'Longitude': 81.34081723209202}, 'P_47': {'Latitude': 16.871322696647127, 'Longitude': 81.14539247686714}, 'P_12': {'Latitude': 16.574151775445337, 'Longitude': 81.1981059588846}, 'P_2': {'Longitude': 80.98368355144568, 'Latitude': 16.594046730421766}, 'P_5': {'Latitude': 16.594704777565404, 'Longitude': 81.06608101238318}, 'P_15': {'Latitude': 16.59455229999004, 'Longitude': 81.26196399111116}, 'P_35': {'Longitude': 81.33902651007027, 'Latitude': 17.03880652638177}, 'P_27': {'Longitude': 81.28272157842964, 'Latitude': 16.919059530368393}} # ####Extract bounding boxes from Eluru for pulling satellite images and store these bounding boxes in a json file.########## # bounding_boxes = extract_bounding_boxes_from_city_coordinates(eluru_coordinates_ori, file_path) #########If .json file for bounding boxes is already available, run the following code################# file_path = "dict_grid_bboxes.json" # Opening JSON file with open(file_path) as json_file: bounding_boxes= json.load(json_file) ########################################################################################################################################### # Path to the trained model file path_model_file = "yolov8_palm_detection_best.pt" num_elements = len(bounding_boxes.keys()) file_ran_successfully = 0 satellite_image_folder = "satellite_images_test/" output_file = 'Eluru_palm_tree_polygon_coordinates.geojson' import time import traceback from itertools import islice # Assuming polygon_coordinates is a list of numpy arrays of pixel coordinates polygon_coordinates_geo = [] total_detect_count = 0 for key, bbox in islice(bounding_boxes.items(), num_elements): try: file_path = 'eluru_meta.dat' # Specify the file name file_name = "output.dat"# Open the file in write mode ('w') with open(file_path, 'w') as file: # Write the data to the file file.write(f"{total_detect_count}, {file_ran_successfully}") # if len(polygon_coordinates_geo) > 0: # save_as_geojson(polygon_coordinates_geo[(len(polygon_coordinates_geo) -1)], (f"intermediate_files/{file_ran_successfully}_palm.geojson")) #print(total_detect_count) image_path = f"{satellite_image_folder}{key}.png" # output_path = f"/content/drive/MyDrive/Eluru/polygon_detection/{key}.png" fieldminlong,fieldminlat, fieldmaxlong, fieldmaxlat = bounding_boxes[key] pull_image_leafmap(f"{key}", fieldminlat, fieldmaxlat, fieldminlong, fieldmaxlong, output_folder=satellite_image_folder) image_width, image_height = get_image_dimensions(image_path) top_left_lat = fieldmaxlat # Latitude of top-left corner of the image top_left_lon = fieldminlong # Longitude of top-left corner of the image bottom_right_lat = fieldminlat # Latitude of bottom-right corner of the image bottom_right_lon = fieldmaxlong # Longitude of bottom-right corner of the image # Calculate differences in latitude and longitude lat_diff = abs(top_left_lat - bottom_right_lat) lon_diff = abs(top_left_lon - bottom_right_lon) # Calculate degrees per pixel lat_deg_per_pixel = lat_diff / image_height lon_deg_per_pixel = lon_diff / image_width matched_locations, detect_count = detect_palm_trees(path_model_file, image_path) matched_locations = np.array(matched_locations) total_detect_count = total_detect_count +detect_count if matched_locations.size == 0: #print("No palm tree detected") os.remove(image_path) else: alpha = 0.5 polygon_coordinates = palm_tree_polygon_detection(image_path, matched_locations, eps_def=80, min_samples_def=2, alpha=alpha) for polygon_coords in polygon_coordinates: geo_coords = [] for coordinate in polygon_coords: lon, lat = save_pixel_as_lat_lon(image_path, coordinate[1], coordinate[0]) geo_coords.append((lon, lat)) polygon_coordinates_geo.append(geo_coords) os.remove(image_path) file_ran_successfully = file_ran_successfully + 1 except: print(traceback.format_exc()) time.sleep(5) time.sleep(0.25) save_as_geojson(polygon_coordinates_geo, output_file) # Open and read the file