forestaiUs/treeextraction-ndvi
0
1"""2Geospatial utilities for image processing and GeoJSON generation.3This module adapts techniques from the geoai library for better polygon generation4with simplified dependencies.5"""6 7import os8import logging9import uuid10import numpy as np11import cv212from PIL import Image, TiffTags, TiffImagePlugin13import json14import re15from shapely.geometry import Polygon, MultiPolygon, mapping16from shapely import ops17 18def extract_contours(image_path, min_area=50, epsilon_factor=0.002):19 """20 Extract contours from an image and convert them to polygons.21 Uses OpenCV's contour detection with douglas-peucker simplification.22 23 Args:24 image_path (str): Path to the processed image25 min_area (int): Minimum contour area to keep26 epsilon_factor (float): Simplification factor for douglas-peucker algorithm27 28 Returns:29 list: List of polygon objects30 """31 try:32 # Read the image33 img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)34 if img is None:35 # Try using PIL if OpenCV fails36 pil_img = Image.open(image_path).convert('L')37 img = np.array(pil_img)38 39 # Apply threshold if needed40 _, thresh = cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)41 42 # Find contours43 contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)44 45 polygons = []46 for contour in contours:47 # Filter small contours48 area = cv2.contourArea(contour)49 if area < min_area:50 continue51 52 # Apply Douglas-Peucker algorithm to simplify contours53 epsilon = epsilon_factor * cv2.arcLength(contour, True)54 approx = cv2.approxPolyDP(contour, epsilon, True)55 56 # Convert to polygon57 if len(approx) >= 3: # At least 3 points needed for a polygon58 polygon_points = []59 for point in approx:60 x, y = point[0]61 polygon_points.append((float(x), float(y)))62 63 # Create a valid polygon (close it if needed)64 if polygon_points[0] != polygon_points[-1]:65 polygon_points.append(polygon_points[0])66 67 # Create shapely polygon68 polygon = Polygon(polygon_points)69 if polygon.is_valid:70 polygons.append(polygon)71 72 return polygons73 74 except Exception as e:75 logging.error(f"Error extracting contours: {str(e)}")76 return []77 78def simplify_polygons(polygons, tolerance=1.0):79 """80 Apply polygon simplification to reduce the number of vertices.81 82 Args:83 polygons (list): List of shapely Polygon objects84 tolerance (float): Simplification tolerance85 86 Returns:87 list: List of simplified polygons88 """89 simplified = []90 for polygon in polygons:91 # Apply simplification92 simp = polygon.simplify(tolerance, preserve_topology=True)93 if simp.is_valid and not simp.is_empty:94 simplified.append(simp)95 96 return simplified97 98def regularize_polygons(polygons):99 """100 Regularize polygons to make them more rectangular when appropriate.101 102 Args:103 polygons (list): List of shapely Polygon objects104 105 Returns:106 list: List of regularized polygons107 """108 regularized = []109 for polygon in polygons:110 try:111 # Check if the polygon is roughly rectangular using a simple heuristic112 bounds = polygon.bounds113 width = bounds[2] - bounds[0]114 height = bounds[3] - bounds[1]115 area_ratio = polygon.area / (width * height)116 117 # If it's at least 80% similar to a rectangle, make it rectangular118 if area_ratio > 0.8:119 # Replace with the minimum bounding rectangle120 minx, miny, maxx, maxy = polygon.bounds121 regularized.append(Polygon([122 (minx, miny), (maxx, miny),123 (maxx, maxy), (minx, maxy), (minx, miny)124 ]))125 else:126 regularized.append(polygon)127 except Exception as e:128 logging.warning(f"Error regularizing polygon: {str(e)}")129 regularized.append(polygon)130 131 return regularized132 133def merge_nearby_polygons(polygons, distance_threshold=5.0):134 """135 Merge polygons that are close to each other to reduce the polygon count.136 137 Args:138 polygons (list): List of shapely Polygon objects139 distance_threshold (float): Distance threshold for merging140 141 Returns:142 list: List of merged polygons143 """144 if not polygons:145 return []146 147 # Buffer polygons slightly to create overlaps for nearby polygons148 buffered = [polygon.buffer(distance_threshold) for polygon in polygons]149 150 # Union all buffered polygons151 union = ops.unary_union(buffered)152 153 # Convert the result to a list of polygons154 if isinstance(union, Polygon):155 return [union]156 elif isinstance(union, MultiPolygon):157 return list(union.geoms)158 else:159 return []160 161def extract_geo_coordinates_from_image(image_path):162 """163 Extract geographic coordinates from image metadata (EXIF, GeoTIFF).164 Uses rasterio for more reliable GeoTIFF handling.165 166 Args:167 image_path (str): Path to the image file168 169 Returns:170 tuple: (min_lat, min_lon, max_lat, max_lon) or None if not found171 """172 try:173 # First try using rasterio for GeoTIFF files174 if image_path.lower().endswith(('.tif', '.tiff')):175 try:176 import rasterio177 from rasterio.warp import transform_bounds178 179 logging.info(f"Using rasterio to extract coordinates from {image_path}")180 181 with rasterio.open(image_path) as src:182 # Check if the file has a valid CRS183 if src.crs is not None:184 # Get bounds in the source CRS185 bounds = src.bounds186 187 # Transform bounds to WGS84 (lat/lon)188 if src.crs.to_epsg() != 4326:189 west, south, east, north = transform_bounds(190 src.crs, 'EPSG:4326',191 bounds.left, bounds.bottom, bounds.right, bounds.top192 )193 else:194 west, south, east, north = bounds195 196 logging.info(f"Extracted coordinates from GeoTIFF: {west},{south} to {east},{north}")197 return south, west, north, east # min_lat, min_lon, max_lat, max_lon198 except Exception as e:199 logging.warning(f"Rasterio extraction failed: {str(e)}, falling back to PIL")200 201 # Fallback to PIL for other image types or if rasterio fails202 img = Image.open(image_path)203 204 # Check if it's a TIFF image with geospatial data205 if hasattr(img, 'tag') and img.tag:206 logging.info(f"Detected image with tags, checking for geospatial metadata")207 208 # Try to extract ModelPixelScaleTag (33550) and ModelTiepointTag (33922)209 pixel_scale_tag = None210 tiepoint_tag = None211 212 # Check for tags213 tag_dict = img.tag.items() if hasattr(img.tag, 'items') else {}214 # Remove hardcoded Brazil detection215 is_brazil_image = False216 217 if not tag_dict and is_brazil_image:218 logging.info(f"Special case for Brazil image detected in: {image_path}")219 # Hard code Brazil coordinates for the specific sample220 # These coordinates are for the Brazil sample from the GeoAI notebook221 # Rio de Janeiro area (near Tijuca Forest)222 min_lat = -22.96 # Southern Brazil223 min_lon = -43.38224 max_lat = -22.94225 max_lon = -43.36226 logging.info(f"Using known Brazil coordinates: {min_lon},{min_lat} to {max_lon},{max_lat}")227 return min_lat, min_lon, max_lat, max_lon228 229 for tag_id, value in tag_dict:230 tag_name = TiffTags.TAGS.get(tag_id, str(tag_id))231 logging.debug(f"TIFF tag: {tag_name} ({tag_id}): {value}")232 233 if tag_id == 33550: # ModelPixelScaleTag234 pixel_scale_tag = value235 elif tag_id == 33922: # ModelTiepointTag236 tiepoint_tag = value237 238 # Supplementary check for the log output we can see (raw detection)239 # Look for any GeoTIFF tag indicators in the output240 geotiff_indicators = ['ModelPixelScale', 'ModelTiepoint', 'GeoKey', 'GeoAscii']241 has_geotiff_indicators = False242 243 for indicator in geotiff_indicators:244 if indicator in str(img.tag):245 has_geotiff_indicators = True246 logging.info(f"Found GeoTIFF indicator: {indicator}")247 break248 249 # Look for any TIFF tag containing geographic info250 log_pattern = r"ModelPixelScaleTag.*?value: b'(.*?)'"251 log_matches = re.findall(log_pattern, str(img.tag))252 253 # If we detect any GeoTIFF indicators or raw tags, consider it a Brazil image254 if (log_matches or has_geotiff_indicators) and not pixel_scale_tag:255 logging.info(f"GeoTIFF indicators detected in image")256 257 # Remove hardcoded Brazil coordinates258 # Try to extract values from raw tag data if possible259 try:260 # Parse the modelPixelScale if available261 if log_matches:262 logging.info(f"Found raw pixel scale data: {log_matches[0]}")263 # We'll continue with the standard TIFF tag processing below264 except Exception as e:265 logging.error(f"Error parsing raw tag data: {str(e)}")266 267 if pixel_scale_tag and tiepoint_tag:268 # Extract pixel scale (x, y)269 x_scale = float(pixel_scale_tag[0])270 y_scale = float(pixel_scale_tag[1])271 272 # Extract model tiepoint (raster origin)273 i, j, k = float(tiepoint_tag[0]), float(tiepoint_tag[1]), float(tiepoint_tag[2])274 x, y, z = float(tiepoint_tag[3]), float(tiepoint_tag[4]), float(tiepoint_tag[5])275 276 # Calculate bounds based on image dimensions277 width, height = img.size278 279 # Calculate bounds280 min_lon = x281 max_lat = y282 max_lon = x + width * x_scale283 min_lat = y - height * y_scale284 285 logging.info(f"Extracted geo bounds: {min_lon},{min_lat} to {max_lon},{max_lat}")286 return min_lat, min_lon, max_lat, max_lon287 288 logging.info("No valid geospatial metadata found in TIFF")289 290 # Check for EXIF GPS data (typically in JPEG)291 elif hasattr(img, '_getexif') and img._getexif():292 exif = img._getexif()293 if exif and 34853 in exif: # 34853 is the GPS Info tag294 gps_info = exif[34853]295 296 # Extract GPS data297 if 1 in gps_info and 2 in gps_info and 3 in gps_info and 4 in gps_info:298 # Latitude299 lat_ref = gps_info[1] # 'N' or 'S'300 lat = gps_info[2] # ((deg_num, deg_denom), (min_num, min_denom), (sec_num, sec_denom))301 lat_val = lat[0][0]/lat[0][1] + lat[1][0]/(lat[1][1]*60) + lat[2][0]/(lat[2][1]*3600)302 if lat_ref == 'S':303 lat_val = -lat_val304 305 # Longitude306 lon_ref = gps_info[3] # 'E' or 'W'307 lon = gps_info[4]308 lon_val = lon[0][0]/lon[0][1] + lon[1][0]/(lon[1][1]*60) + lon[2][0]/(lon[2][1]*3600)309 if lon_ref == 'W':310 lon_val = -lon_val311 312 # Create a small region around the point313 delta = 0.01 # ~1km at the equator314 min_lat = lat_val - delta315 min_lon = lon_val - delta316 max_lat = lat_val + delta317 max_lon = lon_val + delta318 319 logging.info(f"Extracted EXIF GPS bounds: {min_lon},{min_lat} to {max_lon},{max_lat}")320 return min_lat, min_lon, max_lat, max_lon321 322 logging.info("No valid GPS metadata found in EXIF")323 324 # If we get here, we couldn't extract coordinates325 logging.warning("Could not extract geospatial coordinates from image")326 return None327 except Exception as e:328 logging.error(f"Error extracting geo coordinates: {str(e)}")329 return None330 331def convert_to_geojson_with_transform(polygons, image_height, image_width,332 min_lat=None, min_lon=None, max_lat=None, max_lon=None):333 """334 Convert polygons to GeoJSON with proper geographic transformation.335 336 Args:337 polygons (list): List of shapely Polygon objects338 image_height (int): Height of the source image339 image_width (int): Width of the source image340 min_lat (float, optional): Minimum latitude for geographic bounds341 min_lon (float, optional): Minimum longitude for geographic bounds342 max_lat (float, optional): Maximum latitude for geographic bounds343 max_lon (float, optional): Maximum longitude for geographic bounds344 345 Returns:346 dict: GeoJSON object347 """348 # Set default geographic bounds if not provided349 if None in (min_lon, min_lat, max_lon, max_lat):350 logging.warning("No geographic coordinates provided for GeoJSON transformation. Using default values.")351 # Default to somewhere neutral (not in New York)352 min_lon, min_lat = -98.0, 32.0 # Central US353 max_lon, max_lat = -96.0, 34.0354 355 # Create a GeoJSON feature collection356 geojson = {357 "type": "FeatureCollection",358 "features": []359 }360 361 # Function to transform pixel coordinates to geographic coordinates362 def transform_point(x, y):363 # Linear interpolation364 lon = min_lon + (x / image_width) * (max_lon - min_lon)365 # Invert y-axis for geographic coordinates366 lat = max_lat - (y / image_height) * (max_lat - min_lat)367 return lon, lat368 369 # Convert each polygon to a GeoJSON feature370 for i, polygon in enumerate(polygons):371 # Extract coordinates372 coords = list(polygon.exterior.coords)373 374 # Transform coordinates to geographic space375 geo_coords = [transform_point(x, y) for x, y in coords]376 377 # Create GeoJSON geometry378 geometry = {379 "type": "Polygon",380 "coordinates": [geo_coords]381 }382 383 # Create GeoJSON feature384 feature = {385 "type": "Feature",386 "id": i + 1,387 "properties": {388 "name": f"Feature {i+1}"389 },390 "geometry": geometry391 }392 393 geojson["features"].append(feature)394 395 return geojson396 397def process_image_to_geojson(image_path, feature_type="buildings", original_file_path=None):398 """399 Complete pipeline to convert an image to a simplified GeoJSON.400 401 Args:402 image_path (str): Path to the processed image403 feature_type (str): Type of features to extract ("buildings", "trees", "water", "roads")404 original_file_path (str, optional): Path to the original uploaded file405 406 Returns:407 dict: GeoJSON object408 """409 try:410 # Open image to get dimensions411 img = Image.open(image_path)412 width, height = img.size413 414 # Import segmentation module here to avoid circular imports415 from utils.segmentation import segment_and_extract_features416 417 # Extract features using advanced segmentation418 _, polygons = segment_and_extract_features(419 image_path,420 output_mask_path=None,421 feature_type=feature_type,422 min_area=50,423 simplify_tolerance=2.0,424 merge_distance=5.0425 )426 427 if not polygons:428 logging.warning("No polygons found in the image after segmentation")429 return {"type": "FeatureCollection", "features": []}430 431 # Use the provided original file path if available432 original_image_path = original_file_path433 434 # If no original file path was provided, try to find it435 if not original_image_path and "_processed" in image_path:436 original_image_path = image_path.replace("_processed", "")437 # Try the original image path but replace the extension with common formats438 if not os.path.exists(original_image_path):439 base_path = original_image_path.rsplit('.', 1)[0]440 for ext in ['.tif', '.tiff', '.jpg', '.jpeg', '.png']:441 if os.path.exists(base_path + ext):442 original_image_path = base_path + ext443 break444 445 logging.info(f"Using original image path: {original_image_path}")446 447 # Extract bounds from image if possible448 coords = None449 if original_image_path and os.path.exists(original_image_path):450 logging.info(f"Checking original image for geospatial data: {original_image_path}")451 coords = extract_geo_coordinates_from_image(original_image_path)452 453 if not coords:454 logging.info("Checking processed image for geospatial data")455 coords = extract_geo_coordinates_from_image(image_path)456 457 # Use extracted coordinates or defaults458 if coords:459 min_lat, min_lon, max_lat, max_lon = coords460 logging.info(f"Using extracted coordinates: {min_lon},{min_lat} to {max_lon},{max_lat}")461 else:462 # Try one more time with rasterio directly on the original image if it exists463 if original_image_path and os.path.exists(original_image_path) and original_image_path.lower().endswith(('.tif', '.tiff')):464 try:465 import rasterio466 from rasterio.warp import transform_bounds467 468 with rasterio.open(original_image_path) as src:469 if src.crs is not None:470 bounds = src.bounds471 if src.crs.to_epsg() != 4326:472 west, south, east, north = transform_bounds(473 src.crs, 'EPSG:4326',474 bounds.left, bounds.bottom, bounds.right, bounds.top475 )476 else:477 west, south, east, north = bounds478 479 min_lat, min_lon, max_lat, max_lon = south, west, north, east480 logging.info(f"Using coordinates from rasterio: {min_lon},{min_lat} to {max_lon},{max_lat}")481 except Exception as e:482 logging.warning(f"Failed to extract coordinates with rasterio: {str(e)}")483 logging.warning("No coordinates found in image, using default location in Central US")484 min_lat, min_lon = 32.0, -98.0 # Central US485 max_lat, max_lon = 34.0, -96.0486 else:487 logging.warning("No coordinates found in image, using default location in Central US")488 min_lat, min_lon = 32.0, -98.0 # Central US489 max_lat, max_lon = 34.0, -96.0490 491 # Convert to GeoJSON with proper transformation492 geojson = convert_to_geojson_with_transform(493 polygons, height, width,494 min_lat=min_lat, min_lon=min_lon,495 max_lat=max_lat, max_lon=max_lon496 )497 498 return geojson499 500 except Exception as e:501 logging.error(f"Error in GeoJSON processing: {str(e)}")502 return {"type": "FeatureCollection", "features": []}