voxcity.geoprocessor.utils ========================== .. py:module:: voxcity.geoprocessor.utils .. autoapi-nested-parse:: Utility functions for geographic operations and coordinate transformations. This module provides various utility functions for working with geographic data, including coordinate transformations, distance calculations, geocoding, and building polygon processing. It supports operations such as: - Tile coordinate calculations and quadkey conversions - Geographic distance calculations (Haversine and geodetic) - Coordinate system transformations - Polygon and GeoDataFrame operations - Raster file processing and merging - Geocoding and reverse geocoding - Timezone and location information retrieval - Building polygon validation and processing The module uses several external libraries for geographic operations: - pyproj: For coordinate transformations and geodetic calculations - geopandas: For handling geographic data frames - rasterio: For raster file operations - shapely: For geometric operations - geopy: For geocoding services - timezonefinder: For timezone lookups Attributes ---------- .. autoapisummary:: voxcity.geoprocessor.utils.floor_height voxcity.geoprocessor.utils.logger Functions --------- .. autoapisummary:: voxcity.geoprocessor.utils.tile_from_lat_lon voxcity.geoprocessor.utils.quadkey_to_tile voxcity.geoprocessor.utils.initialize_geod voxcity.geoprocessor.utils.calculate_distance voxcity.geoprocessor.utils.normalize_to_one_meter voxcity.geoprocessor.utils.setup_transformer voxcity.geoprocessor.utils.transform_coords voxcity.geoprocessor.utils.create_polygon voxcity.geoprocessor.utils.create_geodataframe voxcity.geoprocessor.utils.haversine_distance voxcity.geoprocessor.utils.get_raster_bbox voxcity.geoprocessor.utils.raster_intersects_polygon voxcity.geoprocessor.utils.save_raster voxcity.geoprocessor.utils.merge_geotiffs voxcity.geoprocessor.utils.convert_format_lat_lon voxcity.geoprocessor.utils.get_coordinates_from_cityname voxcity.geoprocessor.utils.get_city_country_name_from_rectangle voxcity.geoprocessor.utils.get_timezone_info voxcity.geoprocessor.utils.validate_polygon_coordinates voxcity.geoprocessor.utils.create_building_polygons voxcity.geoprocessor.utils.get_country_name Module Contents --------------- .. py:data:: floor_height :value: 2.5 .. py:data:: logger .. py:function:: tile_from_lat_lon(lat, lon, level_of_detail) Convert latitude/longitude coordinates to tile coordinates at a given zoom level. Uses the Web Mercator projection (EPSG:3857) commonly used in web mapping. :param lat: Latitude in degrees (-90 to 90) :type lat: float :param lon: Longitude in degrees (-180 to 180) :type lon: float :param level_of_detail: Zoom level (0-23, where 0 is the entire world) :type level_of_detail: int :returns: (tile_x, tile_y) tile coordinates in the global tile grid :rtype: tuple .. rubric:: Example >>> tile_x, tile_y = tile_from_lat_lon(35.6762, 139.6503, 12) # Tokyo at zoom 12 .. py:function:: quadkey_to_tile(quadkey) Convert a quadkey string to tile coordinates. A quadkey is a string of digits (0-3) that identifies a tile at a certain zoom level. Each digit in the quadkey represents a tile at a zoom level, with each subsequent digit representing a more detailed zoom level. The quadkey numbering scheme: - 0: Top-left quadrant - 1: Top-right quadrant - 2: Bottom-left quadrant - 3: Bottom-right quadrant :param quadkey: Quadkey string (e.g., "120" for zoom level 3) :type quadkey: str :returns: (tile_x, tile_y, level_of_detail) tile coordinates and zoom level :rtype: tuple .. rubric:: Example >>> x, y, zoom = quadkey_to_tile("120") # Returns coordinates at zoom level 3 .. py:function:: initialize_geod() Return the module-level WGS84 Geod singleton for geodetic calculations. The Geod object provides methods for: - Forward geodetic calculations (direct) - Inverse geodetic calculations (inverse) - Area calculations - Line length calculations :returns: Initialized Geod object for WGS84 calculations :rtype: Geod .. rubric:: Example >>> geod = initialize_geod() >>> fwd_az, back_az, dist = geod.inv(lon1, lat1, lon2, lat2) .. py:function:: calculate_distance(geod, lon1, lat1, lon2, lat2) Calculate geodetic distance between two points on the Earth's surface. Uses inverse geodetic computation to find the shortest distance along the ellipsoid, which is more accurate than great circle (spherical) calculations. :param geod: Geod object for calculations, initialized with WGS84 :type geod: Geod :param lon1: Coordinates of first point in decimal degrees :type lon1: float :param lat1: Coordinates of first point in decimal degrees :type lat1: float :param lon2: Coordinates of second point in decimal degrees :type lon2: float :param lat2: Coordinates of second point in decimal degrees :type lat2: float :returns: Distance in meters between the two points along the ellipsoid :rtype: float .. rubric:: Example >>> geod = initialize_geod() >>> distance = calculate_distance(geod, 139.6503, 35.6762, ... -74.0060, 40.7128) # Tokyo to NYC .. py:function:: normalize_to_one_meter(vector, distance_in_meters) Normalize a vector to represent one meter in geographic space. Useful for creating unit vectors in geographic calculations, particularly when working with distance-based operations or scaling geographic features. :param vector: Vector to normalize, typically a direction vector :type vector: numpy.ndarray :param distance_in_meters: Current distance in meters that the vector represents :type distance_in_meters: float :returns: Normalized vector where magnitude represents 1 meter :rtype: numpy.ndarray .. rubric:: Example >>> direction = np.array([3.0, 4.0]) # Vector of length 5 >>> unit_meter = normalize_to_one_meter(direction, 5.0) .. py:function:: setup_transformer(from_crs, to_crs) Set up a coordinate transformer between two Coordinate Reference Systems (CRS). Results are cached to avoid repeated expensive Transformer construction. The always_xy=True parameter ensures consistent handling of coordinate order by always using (x,y) or (longitude,latitude) order regardless of CRS definition. Common CRS codes: - EPSG:4326 - WGS84 (latitude/longitude) - EPSG:3857 - Web Mercator - EPSG:2263 - NY State Plane :param from_crs: Source coordinate reference system (EPSG code, proj4 string, or CRS dict) :param to_crs: Target coordinate reference system (EPSG code, proj4 string, or CRS dict) :returns: Initialized transformer object for coordinate conversion :rtype: Transformer .. rubric:: Example >>> transformer = setup_transformer("EPSG:4326", "EPSG:3857") >>> x, y = transformer.transform(longitude, latitude) .. py:function:: transform_coords(transformer, lon, lat) Transform coordinates using provided transformer with error handling. Includes validation for infinite values that may result from invalid transformations or coordinates outside the valid range for the target CRS. :param transformer: Coordinate transformer from setup_transformer() :type transformer: Transformer :param lon: Input coordinates in the source CRS :type lon: float :param lat: Input coordinates in the source CRS :type lat: float :returns: (x, y) transformed coordinates in the target CRS, or (None, None) if transformation fails :rtype: tuple .. rubric:: Example >>> transformer = setup_transformer("EPSG:4326", "EPSG:3857") >>> x, y = transform_coords(transformer, -74.0060, 40.7128) # NYC coordinates >>> if x is not None: ... print(f"Transformed coordinates: ({x}, {y})") .. py:function:: create_polygon(vertices) Create a Shapely polygon from a list of vertices. Input vertices must be in (longitude, latitude) format as required by Shapely. The polygon will be automatically closed if the first and last vertices don't match. :param vertices: List of (longitude, latitude) coordinate pairs forming the polygon. The coordinates should be in counter-clockwise order for exterior rings and clockwise order for interior rings (holes). :type vertices: list :returns: Shapely polygon object that can be used for spatial operations :rtype: Polygon .. rubric:: Example >>> vertices = [(0, 0), (1, 0), (1, 1), (0, 1)] # Square >>> polygon = create_polygon(vertices) >>> print(f"Polygon area: {polygon.area}") .. py:function:: create_geodataframe(polygon, crs=4326) Create a GeoDataFrame from a Shapely polygon. Default CRS is WGS84 (EPSG:4326) for geographic coordinates. The GeoDataFrame provides additional functionality for spatial operations, data analysis, and export to various geographic formats. :param polygon: Shapely polygon object to convert :type polygon: Polygon :param crs: Coordinate reference system EPSG code (default: 4326 for WGS84) :type crs: int :returns: GeoDataFrame containing the polygon with specified CRS :rtype: GeoDataFrame .. rubric:: Example >>> vertices = [(0, 0), (1, 0), (1, 1), (0, 1)] >>> polygon = create_polygon(vertices) >>> gdf = create_geodataframe(polygon) >>> gdf.to_file("polygon.geojson", driver="GeoJSON") .. py:function:: haversine_distance(lon1, lat1, lon2, lat2) Calculate great-circle distance between two points using Haversine formula. This is an approximation that treats the Earth as a perfect sphere. :param lon1: Coordinates of first point :type lon1: float :param lat1: Coordinates of first point :type lat1: float :param lon2: Coordinates of second point :type lon2: float :param lat2: Coordinates of second point :type lat2: float :returns: Distance in kilometers :rtype: float .. py:function:: get_raster_bbox(raster_path) Get the bounding box of a raster file in its native coordinate system. Returns a rectangular polygon representing the spatial extent of the raster, which can be used for spatial queries and intersection tests. :param raster_path: Path to the raster file (GeoTIFF, IMG, etc.) :type raster_path: str :returns: Shapely box representing the raster bounds in the raster's CRS :rtype: box .. rubric:: Example >>> bbox = get_raster_bbox("elevation.tif") >>> print(f"Raster extent: {bbox.bounds}") # (minx, miny, maxx, maxy) .. py:function:: raster_intersects_polygon(raster_path, polygon) Check if a raster file's extent intersects with a given polygon. Automatically handles coordinate system transformations by converting the raster bounds to WGS84 (EPSG:4326) if needed before the intersection test. :param raster_path: Path to the raster file to check :type raster_path: str :param polygon: Shapely polygon to test intersection with (in WGS84) :type polygon: Polygon :returns: True if raster intersects or contains the polygon, False otherwise :rtype: bool .. rubric:: Example >>> aoi = create_polygon([(lon1, lat1), (lon2, lat2), ...]) # Area of interest >>> if raster_intersects_polygon("dem.tif", aoi): ... print("Raster covers the area of interest") .. py:function:: save_raster(input_path, output_path) Create a copy of a raster file at a new location. Performs a direct file copy without any transformation or modification, preserving all metadata, georeferencing, and pixel values. :param input_path: Source raster file path :type input_path: str :param output_path: Destination path for the copied raster :type output_path: str .. rubric:: Example >>> save_raster("original.tif", "backup/copy.tif") >>> print("Copied original file to: backup/copy.tif") .. py:function:: merge_geotiffs(geotiff_files, output_dir) Merge multiple GeoTIFF files into a single mosaic. Handles edge matching and overlapping areas between adjacent rasters. The output will have the same coordinate system and data type as the input files. Important considerations: - All input files should have the same coordinate system - All input files should have the same data type - Overlapping areas are handled by taking the first value encountered :param geotiff_files: List of paths to GeoTIFF files to merge :type geotiff_files: list :param output_dir: Directory where the merged output will be saved :type output_dir: str .. rubric:: Example >>> files = ["tile1.tif", "tile2.tif", "tile3.tif"] >>> merge_geotiffs(files, "output_directory") >>> print("Merged output saved to: output_directory/lulc.tif") .. py:function:: convert_format_lat_lon(input_coords) Convert coordinate format and close polygon. Input coordinates are already in [lon, lat] format. :param input_coords: List of [lon, lat] coordinates :type input_coords: list :returns: List of [lon, lat] coordinates with first point repeated at end :rtype: list .. py:function:: get_coordinates_from_cityname(place_name) Geocode a city name to get its coordinates using OpenStreetMap's Nominatim service. Includes rate limiting and error handling to comply with Nominatim's usage policy. Note: - Results may vary based on the specificity of the place name - For better results, include country or state information - Service has usage limits and may timeout :param place_name: Name of the city to geocode (e.g., "Tokyo, Japan") :type place_name: str :returns: (latitude, longitude) coordinates or None if geocoding fails :rtype: tuple .. rubric:: Example >>> coords = get_coordinates_from_cityname("Paris, France") >>> if coords: ... lat, lon = coords ... print(f"Paris coordinates: {lat}, {lon}") .. py:function:: get_city_country_name_from_rectangle(coordinates) Get the city and country name for a location defined by a rectangle. Uses reverse geocoding to find the nearest named place to the rectangle's center. The function: 1. Calculates the center point of the rectangle 2. Performs reverse geocoding with rate limiting 3. Extracts city and country information from the result :param coordinates: List of (longitude, latitude) coordinates defining the rectangle :type coordinates: list :returns: String in format "city/ country" or fallback value if lookup fails :rtype: str .. rubric:: Example >>> coords = [(139.65, 35.67), (139.66, 35.67), ... (139.66, 35.68), (139.65, 35.68)] >>> location = get_city_country_name_from_rectangle(coords) >>> print(f"Location: {location}") # e.g., "Shibuya/ Japan" .. py:function:: get_timezone_info(rectangle_coords) Get timezone and central meridian information for a location. Uses the rectangle's center point to determine the local timezone and calculates the central meridian based on the UTC offset. The function provides: 1. Local timezone identifier (e.g., "America/New_York") 2. UTC offset (e.g., "UTC-04:00") 3. Central meridian longitude for the timezone :param rectangle_coords: List of (longitude, latitude) coordinates defining the area :type rectangle_coords: list :returns: (timezone string with UTC offset, central meridian longitude string) :rtype: tuple .. rubric:: Example >>> coords = [(139.65, 35.67), (139.66, 35.67), ... (139.66, 35.68), (139.65, 35.68)] >>> tz, meridian = get_timezone_info(coords) >>> print(f"Timezone: {tz}, Meridian: {meridian}") # e.g., "UTC+09:00, 135.00000" .. py:function:: validate_polygon_coordinates(geometry) Validate and ensure proper closure of polygon coordinate rings. Performs validation and correction of GeoJSON polygon geometries according to the GeoJSON specification requirements. Validation checks: 1. Geometry type (Polygon or MultiPolygon) 2. Ring closure (first point equals last point) 3. Minimum number of points (4, including closure) :param geometry: GeoJSON geometry object with 'type' and 'coordinates' properties :type geometry: dict :returns: True if polygon coordinates are valid or were successfully corrected, False if validation failed :rtype: bool .. rubric:: Example >>> geom = { ... "type": "Polygon", ... "coordinates": [[[0,0], [1,0], [1,1], [0,1]]] # Not closed ... } >>> if validate_polygon_coordinates(geom): ... print("Polygon is valid") # Will close the ring automatically .. py:function:: create_building_polygons(filtered_buildings) Create building polygons with properties from filtered GeoJSON features. Processes a list of GeoJSON building features to create Shapely polygons with associated height and other properties, while also building a spatial index. Processing steps: 1. Extract and validate coordinates 2. Create Shapely polygons 3. Process building properties (height, levels, etc.) 4. Build spatial index for efficient querying Height calculation rules: - Use explicit height if available - Calculate from levels * floor_height if height not available - Calculate from floors * floor_height if levels not available - Use NaN if no height information available :param filtered_buildings: List of GeoJSON building features with properties :type filtered_buildings: list :returns: ( list of tuples (polygon, height, min_height, is_inner, feature_id), rtree spatial index for the polygons ) :rtype: tuple .. rubric:: Example >>> buildings = [ ... { ... "type": "Feature", ... "geometry": {"type": "Polygon", "coordinates": [...]}, ... "properties": {"height": 30, "levels": 10} ... }, ... # ... more buildings ... ... ] >>> polygons, spatial_idx = create_building_polygons(buildings) .. py:function:: get_country_name(lon, lat) Get country name from coordinates using reverse geocoding. Uses a local database for fast reverse geocoding to country level, then converts the country code to full name using pycountry. Results are cached to avoid repeated lookups for nearby coordinates. :param lon: Longitude in decimal degrees :type lon: float :param lat: Latitude in decimal degrees :type lat: float :returns: Full country name or None if lookup fails :rtype: str .. rubric:: Example >>> country = get_country_name(139.6503, 35.6762) >>> print(f"Country: {country}") # "Japan"