From ed22d1da5640eccfce23a10105994af40d794b04 Mon Sep 17 00:00:00 2001 From: Henry Lydecker Date: Fri, 23 Feb 2024 15:35:25 +1100 Subject: [PATCH 1/3] feature: bbox to geojson tile grid --- aigis/convert/tiles.py | 108 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) diff --git a/aigis/convert/tiles.py b/aigis/convert/tiles.py index c8141f2..80eecea 100644 --- a/aigis/convert/tiles.py +++ b/aigis/convert/tiles.py @@ -1,4 +1,5 @@ # -*- coding: utf-8 -*- +import json import glob import itertools import logging @@ -6,6 +7,8 @@ import rasterio as rio import rasterio.windows as riow +from shapely.geometry import Polygon, box +from pyproj import Proj, transform """A collection of functions for manipulating raster tiles.""" @@ -13,6 +16,111 @@ log = logging.getLogger(__name__) +def create_grid_geojson(bbox, tile_size=1, crs='EPSG:4326', feature_collection_name="GSU_LCZ_subset"): + """ + Create a GeoJSON feature collection representing a grid of tiles within a given bounding box. + + Parameters: + - bbox (tuple): The bounding box coordinates in the format (min_lon, min_lat, max_lon, max_lat). + - tile_size (float): The size of each tile in degrees of latitude and longitude. Default is 1. + - crs (str): The coordinate reference system of the bounding box. Default is 'EPSG:4326'. + - feature_collection_name (str): The name of the feature collection. Default is "GSU_LCZ_subset". + + Returns: + - str: The GeoJSON feature collection as a JSON string. + """ + min_lon, min_lat, max_lon, max_lat = bbox + + # Check if the bounding box is large enough to create a grid + # Convert tile_size to degrees of latitude and longitude + tile_size_deg = tile_size / (111.32 * 1000) # Approximate length of one degree of latitude in meters + + if max_lon - min_lon < tile_size_deg or max_lat - min_lat < tile_size_deg: + # Create a polygon encompassing the extent of the bounding box + bbox_polygon = box(min_lon, min_lat, max_lon, max_lat) + + # Create a GeoJSON feature for the bounding box + feature = { + "type": "Feature", + "properties": { + "id": 1, + "left": min_lon, + "top": max_lat, + "right": max_lon, + "bottom": min_lat + }, + "geometry": bbox_polygon.__geo_interface__ + } + + # Create a GeoJSON feature collection with the bounding box feature + feature_collection = { + "type": "FeatureCollection", + "name": feature_collection_name, + "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}}, + "features": [feature] + } + + return json.dumps(feature_collection) + + # Transform coordinates if needed + if crs != 'EPSG:4326': + in_proj = Proj(init='epsg:4326') + out_proj = Proj(init=crs) + min_lon, min_lat = transform(in_proj, out_proj, min_lon, min_lat) + max_lon, max_lat = transform(in_proj, out_proj, max_lon, max_lat) + + # Calculate the number of tiles in x and y directions + num_tiles_x = int((max_lon - min_lon) / tile_size) + num_tiles_y = int((max_lat - min_lat) / tile_size) + + features = [] + + for i in range(num_tiles_x): + for j in range(num_tiles_y): + # Calculate the coordinates of the current tile + tile_min_lon = min_lon + i * tile_size + tile_max_lon = min_lon + (i + 1) * tile_size + tile_min_lat = min_lat + j * tile_size + tile_max_lat = min_lat + (j + 1) * tile_size + + # Create a polygon for the current tile + tile_polygon = box(tile_min_lon, tile_min_lat, tile_max_lon, tile_max_lat) + + # Transform back to EPSG:4326 if needed + if crs != 'EPSG:4326': + in_proj = Proj(init=crs) + out_proj = Proj(init='epsg:4326') + tile_polygon = transform(in_proj, out_proj, tile_polygon) + + # Create a GeoJSON feature for the current tile with the desired properties + feature = { + "type": "Feature", + "properties": { + "LCZ": 6.0, + "id": f"{i}{j}", + "left": tile_min_lon, + "top": tile_max_lat, + "right": tile_max_lon, + "bottom": tile_min_lat, + "row_index": j, + "col_index": i + }, + "geometry": tile_polygon.__geo_interface__ + } + + features.append(feature) + + # Create a GeoJSON feature collection + feature_collection = { + "type": "FeatureCollection", + "name": feature_collection_name, + "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}}, + "features": features + } + + return json.dumps(feature_collection) + + def get_tiles( geotiff: rio.DatasetReader, tile_width: int = 2000, From bb39c8924ccacc265eaf1364ee6f02df194725e8 Mon Sep 17 00:00:00 2001 From: Henry Lydecker Date: Fri, 23 Feb 2024 16:43:57 +1100 Subject: [PATCH 2/3] Update tiles.py --- aigis/convert/tiles.py | 91 ++++++++++++++---------------------------- 1 file changed, 31 insertions(+), 60 deletions(-) diff --git a/aigis/convert/tiles.py b/aigis/convert/tiles.py index 80eecea..0d3c225 100644 --- a/aigis/convert/tiles.py +++ b/aigis/convert/tiles.py @@ -16,58 +16,24 @@ log = logging.getLogger(__name__) -def create_grid_geojson(bbox, tile_size=1, crs='EPSG:4326', feature_collection_name="GSU_LCZ_subset"): +def create_grid_geojson(bbox, tile_size): """ - Create a GeoJSON feature collection representing a grid of tiles within a given bounding box. - + Create a GeoJSON representation of a grid of tiles within the given bounding box. + Parameters: - - bbox (tuple): The bounding box coordinates in the format (min_lon, min_lat, max_lon, max_lat). - - tile_size (float): The size of each tile in degrees of latitude and longitude. Default is 1. - - crs (str): The coordinate reference system of the bounding box. Default is 'EPSG:4326'. - - feature_collection_name (str): The name of the feature collection. Default is "GSU_LCZ_subset". - + bbox (tuple): A tuple containing the minimum and maximum longitude and latitude values of the bounding box. + tile_size (float): The size of each tile in degrees. + Returns: - - str: The GeoJSON feature collection as a JSON string. + str: A JSON string representing the GeoJSON feature collection. """ min_lon, min_lat, max_lon, max_lat = bbox - # Check if the bounding box is large enough to create a grid - # Convert tile_size to degrees of latitude and longitude - tile_size_deg = tile_size / (111.32 * 1000) # Approximate length of one degree of latitude in meters - - if max_lon - min_lon < tile_size_deg or max_lat - min_lat < tile_size_deg: - # Create a polygon encompassing the extent of the bounding box - bbox_polygon = box(min_lon, min_lat, max_lon, max_lat) - - # Create a GeoJSON feature for the bounding box - feature = { - "type": "Feature", - "properties": { - "id": 1, - "left": min_lon, - "top": max_lat, - "right": max_lon, - "bottom": min_lat - }, - "geometry": bbox_polygon.__geo_interface__ - } - - # Create a GeoJSON feature collection with the bounding box feature - feature_collection = { - "type": "FeatureCollection", - "name": feature_collection_name, - "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}}, - "features": [feature] - } - - return json.dumps(feature_collection) - - # Transform coordinates if needed - if crs != 'EPSG:4326': - in_proj = Proj(init='epsg:4326') - out_proj = Proj(init=crs) - min_lon, min_lat = transform(in_proj, out_proj, min_lon, min_lat) - max_lon, max_lat = transform(in_proj, out_proj, max_lon, max_lat) + # Transform the coordinates to EPSG:3857 (Web Mercator) for easier calculations + in_proj = Proj(init='epsg:4326') + out_proj = Proj(init='epsg:3857') + min_lon, min_lat = transform(in_proj, out_proj, min_lon, min_lat) + max_lon, max_lat = transform(in_proj, out_proj, max_lon, max_lat) # Calculate the number of tiles in x and y directions num_tiles_x = int((max_lon - min_lon) / tile_size) @@ -83,29 +49,29 @@ def create_grid_geojson(bbox, tile_size=1, crs='EPSG:4326', feature_collection_n tile_min_lat = min_lat + j * tile_size tile_max_lat = min_lat + (j + 1) * tile_size + # Convert the coordinates back to EPSG:4326 + tile_min_lon, tile_min_lat = transform(out_proj, in_proj, tile_min_lon, tile_min_lat) + tile_max_lon, tile_max_lat = transform(out_proj, in_proj, tile_max_lon, tile_max_lat) + # Create a polygon for the current tile tile_polygon = box(tile_min_lon, tile_min_lat, tile_max_lon, tile_max_lat) - # Transform back to EPSG:4326 if needed - if crs != 'EPSG:4326': - in_proj = Proj(init=crs) - out_proj = Proj(init='epsg:4326') - tile_polygon = transform(in_proj, out_proj, tile_polygon) - - # Create a GeoJSON feature for the current tile with the desired properties + # Create a GeoJSON feature for the current tile feature = { "type": "Feature", "properties": { - "LCZ": 6.0, - "id": f"{i}{j}", + "id": i * num_tiles_y + j, "left": tile_min_lon, "top": tile_max_lat, "right": tile_max_lon, "bottom": tile_min_lat, - "row_index": j, - "col_index": i + "row_index": i, + "col_index": j }, - "geometry": tile_polygon.__geo_interface__ + "geometry": { + "type": "Polygon", + "coordinates": [list(tile_polygon.exterior.coords)] + } } features.append(feature) @@ -113,8 +79,13 @@ def create_grid_geojson(bbox, tile_size=1, crs='EPSG:4326', feature_collection_n # Create a GeoJSON feature collection feature_collection = { "type": "FeatureCollection", - "name": feature_collection_name, - "crs": {"type": "name", "properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}}, + "name": "GSU_grid_1", + "crs": { + "type": "name", + "properties": { + "name": "urn:ogc:def:crs:EPSG::3857" + } + }, "features": features } From 2b2405124a089226e68f58ae1a738c497ad78a7d Mon Sep 17 00:00:00 2001 From: Henry Lydecker Date: Fri, 23 Feb 2024 16:44:40 +1100 Subject: [PATCH 3/3] Update coco2geojson.py --- scripts/coco2geojson.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/scripts/coco2geojson.py b/scripts/coco2geojson.py index 419e36c..773e5a1 100644 --- a/scripts/coco2geojson.py +++ b/scripts/coco2geojson.py @@ -486,5 +486,28 @@ def main(args=None): polygons_df.to_parquet(geopardquet_path) -if __name__ == "__main__": - main() + if len(tiles_list) == 1: + # Handle single input tile + single_tile_df = tiles_df.loc[tiles_df["tile_name"] == os.path.basename(tiles_list[0]).split(".")[0]] + if single_tile_df.empty: + log.warning("No annotations found for the single input tile.") + else: + single_tile_polygons_df = merge_class_polygons_shapely([single_tile_df], crs) + single_tile_polygons_df = multipolygon_to_polygons(single_tile_polygons_df) + single_tile_polygons_df["geometry"] = single_tile_polygons_df["geometry"].progress_apply( + lambda x: shape_regulariser( + x, simplify_tolerance, minimum_rotated_rectangle, orthogonalisation + ) + ) + single_tile_polygons_df = single_tile_polygons_df.to_crs(original_crs) + single_tile_polygons_df.Name = meta_name + + if geojson_path is not None: + single_tile_polygons_df.to_file(geojson_path, driver="GeoJSON") + if geopardquet_path is not None: + single_tile_polygons_df.to_parquet(geopardquet_path) + else: + log.warning("Multiple input tiles detected. The code will handle overlapping and adjacent tiles by default.") + + if __name__ == "__main__": + main()