voxcity.geoprocessor.utils

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

floor_height
logger

Functions

tile_from_lat_lon(lat, lon, level_of_detail)	Convert latitude/longitude coordinates to tile coordinates at a given zoom level.
quadkey_to_tile(quadkey)	Convert a quadkey string to tile coordinates.
initialize_geod()	Return the module-level WGS84 Geod singleton for geodetic calculations.
calculate_distance(geod, lon1, lat1, lon2, lat2)	Calculate geodetic distance between two points on the Earth's surface.
normalize_to_one_meter(vector, distance_in_meters)	Normalize a vector to represent one meter in geographic space.
setup_transformer(from_crs, to_crs)	Set up a coordinate transformer between two Coordinate Reference Systems (CRS).
transform_coords(transformer, lon, lat)	Transform coordinates using provided transformer with error handling.
create_polygon(vertices)	Create a Shapely polygon from a list of vertices.
create_geodataframe(polygon[, crs])	Create a GeoDataFrame from a Shapely polygon.
haversine_distance(lon1, lat1, lon2, lat2)	Calculate great-circle distance between two points using Haversine formula.
get_raster_bbox(raster_path)	Get the bounding box of a raster file in its native coordinate system.
raster_intersects_polygon(raster_path, polygon)	Check if a raster file's extent intersects with a given polygon.
save_raster(input_path, output_path)	Create a copy of a raster file at a new location.
merge_geotiffs(geotiff_files, output_dir)	Merge multiple GeoTIFF files into a single mosaic.
convert_format_lat_lon(input_coords)	Convert coordinate format and close polygon.
get_coordinates_from_cityname(place_name)	Geocode a city name to get its coordinates using OpenStreetMap's Nominatim service.
get_city_country_name_from_rectangle(coordinates)	Get the city and country name for a location defined by a rectangle.
get_timezone_info(rectangle_coords)	Get timezone and central meridian information for a location.
validate_polygon_coordinates(geometry)	Validate and ensure proper closure of polygon coordinate rings.
create_building_polygons(filtered_buildings)	Create building polygons with properties from filtered GeoJSON features.
get_country_name(lon, lat)	Get country name from coordinates using reverse geocoding.

Module Contents

voxcity.geoprocessor.utils.floor_height = 2.5

voxcity.geoprocessor.utils.logger

voxcity.geoprocessor.utils.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.

Parameters:

• lat (float) – Latitude in degrees (-90 to 90)

• lon (float) – Longitude in degrees (-180 to 180)

• level_of_detail (int) – Zoom level (0-23, where 0 is the entire world)

Returns:

(tile_x, tile_y) tile coordinates in the global tile grid

Return type:

tuple

Example

>>> tile_x, tile_y = tile_from_lat_lon(35.6762, 139.6503, 12)  # Tokyo at zoom 12

voxcity.geoprocessor.utils.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

Parameters:

quadkey (str) – Quadkey string (e.g., “120” for zoom level 3)

Returns:

(tile_x, tile_y, level_of_detail) tile coordinates and zoom level

Return type:

tuple

Example

>>> x, y, zoom = quadkey_to_tile("120")  # Returns coordinates at zoom level 3

voxcity.geoprocessor.utils.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

Return type:

Geod

Example

>>> geod = initialize_geod()
>>> fwd_az, back_az, dist = geod.inv(lon1, lat1, lon2, lat2)

voxcity.geoprocessor.utils.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.

Parameters:

• geod (Geod) – Geod object for calculations, initialized with WGS84

• lon1 (float) – Coordinates of first point in decimal degrees

• lat1 (float) – Coordinates of first point in decimal degrees

• lon2 (float) – Coordinates of second point in decimal degrees

• lat2 (float) – Coordinates of second point in decimal degrees

Returns:

Distance in meters between the two points along the ellipsoid

Return type:

float

Example

>>> geod = initialize_geod()
>>> distance = calculate_distance(geod, 139.6503, 35.6762,
...                             -74.0060, 40.7128)  # Tokyo to NYC

voxcity.geoprocessor.utils.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.

Parameters:

• vector (numpy.ndarray) – Vector to normalize, typically a direction vector

• distance_in_meters (float) – Current distance in meters that the vector represents

Returns:

Normalized vector where magnitude represents 1 meter

Return type:

numpy.ndarray

Example

>>> direction = np.array([3.0, 4.0])  # Vector of length 5
>>> unit_meter = normalize_to_one_meter(direction, 5.0)

voxcity.geoprocessor.utils.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

Parameters:

• from_crs – Source coordinate reference system (EPSG code, proj4 string, or CRS dict)

• to_crs – Target coordinate reference system (EPSG code, proj4 string, or CRS dict)

Returns:

Initialized transformer object for coordinate conversion

Return type:

Transformer

Example

>>> transformer = setup_transformer("EPSG:4326", "EPSG:3857")
>>> x, y = transformer.transform(longitude, latitude)

voxcity.geoprocessor.utils.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.

Parameters:

• transformer (Transformer) – Coordinate transformer from setup_transformer()

• lon (float) – Input coordinates in the source CRS

• lat (float) – Input coordinates in the source CRS

Returns:

(x, y) transformed coordinates in the target CRS, or (None, None) if transformation fails

Return type:

tuple

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})")

voxcity.geoprocessor.utils.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.

Parameters:

vertices (list) – 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).

Returns:

Shapely polygon object that can be used for spatial operations

Return type:

Polygon

Example

>>> vertices = [(0, 0), (1, 0), (1, 1), (0, 1)]  # Square
>>> polygon = create_polygon(vertices)
>>> print(f"Polygon area: {polygon.area}")

voxcity.geoprocessor.utils.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.

Parameters:

• polygon (Polygon) – Shapely polygon object to convert

• crs (int) – Coordinate reference system EPSG code (default: 4326 for WGS84)

Returns:

GeoDataFrame containing the polygon with specified CRS

Return type:

GeoDataFrame

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")

voxcity.geoprocessor.utils.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.

Parameters:

• lon1 (float) – Coordinates of first point

• lat1 (float) – Coordinates of first point

• lon2 (float) – Coordinates of second point

• lat2 (float) – Coordinates of second point

Returns:

Distance in kilometers

Return type:

float

voxcity.geoprocessor.utils.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.

Parameters:

raster_path (str) – Path to the raster file (GeoTIFF, IMG, etc.)

Returns:

Shapely box representing the raster bounds in the raster’s CRS

Return type:

box

Example

>>> bbox = get_raster_bbox("elevation.tif")
>>> print(f"Raster extent: {bbox.bounds}")  # (minx, miny, maxx, maxy)

voxcity.geoprocessor.utils.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.

Parameters:

• raster_path (str) – Path to the raster file to check

• polygon (Polygon) – Shapely polygon to test intersection with (in WGS84)

Returns:

True if raster intersects or contains the polygon, False otherwise

Return type:

bool

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")

voxcity.geoprocessor.utils.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.

Parameters:

• input_path (str) – Source raster file path

• output_path (str) – Destination path for the copied raster

Example

>>> save_raster("original.tif", "backup/copy.tif")
>>> print("Copied original file to: backup/copy.tif")

voxcity.geoprocessor.utils.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

Parameters:

• geotiff_files (list) – List of paths to GeoTIFF files to merge

• output_dir (str) – Directory where the merged output will be saved

Example

>>> files = ["tile1.tif", "tile2.tif", "tile3.tif"]
>>> merge_geotiffs(files, "output_directory")
>>> print("Merged output saved to: output_directory/lulc.tif")

voxcity.geoprocessor.utils.convert_format_lat_lon(input_coords)

Convert coordinate format and close polygon. Input coordinates are already in [lon, lat] format.

Parameters:

input_coords (list) – List of [lon, lat] coordinates

Returns:

List of [lon, lat] coordinates with first point repeated at end

Return type:

list

voxcity.geoprocessor.utils.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

Parameters:

place_name (str) – Name of the city to geocode (e.g., “Tokyo, Japan”)

Returns:

(latitude, longitude) coordinates or None if geocoding fails

Return type:

tuple

Example

>>> coords = get_coordinates_from_cityname("Paris, France")
>>> if coords:
...     lat, lon = coords
...     print(f"Paris coordinates: {lat}, {lon}")

voxcity.geoprocessor.utils.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

Parameters:

coordinates (list) – List of (longitude, latitude) coordinates defining the rectangle

Returns:

String in format “city/ country” or fallback value if lookup fails

Return type:

str

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"

voxcity.geoprocessor.utils.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

Parameters:

rectangle_coords (list) – List of (longitude, latitude) coordinates defining the area

Returns:

(timezone string with UTC offset, central meridian longitude string)

Return type:

tuple

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"

voxcity.geoprocessor.utils.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)

Parameters:

geometry (dict) – GeoJSON geometry object with ‘type’ and ‘coordinates’ properties

Returns:

True if polygon coordinates are valid or were successfully corrected,

False if validation failed

Return type:

bool

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

voxcity.geoprocessor.utils.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

Parameters:

filtered_buildings (list) – List of GeoJSON building features with properties

Returns:

(

list of tuples (polygon, height, min_height, is_inner, feature_id), rtree spatial index for the polygons

)

Return type:

tuple

Example

>>> buildings = [
...     {
...         "type": "Feature",
...         "geometry": {"type": "Polygon", "coordinates": [...]},
...         "properties": {"height": 30, "levels": 10}
...     },
...     # ... more buildings ...
... ]
>>> polygons, spatial_idx = create_building_polygons(buildings)

voxcity.geoprocessor.utils.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.

Parameters:

• lon (float) – Longitude in decimal degrees

• lat (float) – Latitude in decimal degrees

Returns:

Full country name or None if lookup fails

Return type:

str

Example

>>> country = get_country_name(139.6503, 35.6762)
>>> print(f"Country: {country}")  # "Japan"