๐๏ธ VoxCity Complete Demoยถ
VoxCity is a Python package for generating 3D voxel-based urban models from open geospatial data.
What This Notebook Coversยถ
Environment Setup - Install VoxCity and authenticate Google Earth Engine
Area Selection - Define your study area using coordinates or interactive map
Configuration - Select data sources and parameters
Model Generation - Create 3D voxel city models
Visualization - View and analyze your urban model
Prerequisitesยถ
Python 3.8+
Google Earth Engine account (for some data sources)
Internet connection for data download
โ ๏ธ Important Notesยถ
Execute cells in order from Step 1 to Step 5 by clicking the โถ๏ธ icons
In Step 2, choose only ONE option (Option 1, 2, or 3) for defining your study area
Some data sources require Google Earth Engine authentication
Model generation time depends on area size and voxel resolution
๐ฆ Step 1: Prepare Environmentยถ
Install VoxCity and set up authentication for data sources.
#@title 1.1. Install voxcity
!pip install voxcity
!apt-get update
!apt-get install -y xvfb libgl1-mesa-glx
#@title 1.2. Authenticate Google Earth Engine (Click link, generate token, copy and paste the token)
!earthengine authenticate --auth_mode=notebook
#@title 1.3. Authenticate Google Earth Engine on local environment
!earthengine authenticate
import ee
ee.Authenticate()
ee.Initialize(project='your_project')
๐ Step 2: Define Target Areaยถ
Choose ONE of the following three options to define your study area.
Option |
Method |
Best For |
|---|---|---|
Option 1 |
Manual coordinates |
When you know exact coordinates |
Option 2 |
Draw rectangle |
Interactive selection by corners |
Option 3 |
Center + dimensions |
Define by center point and size |
Option 1: Set Coordinates Directlyยถ
Define the four corners of your study area as (longitude, latitude) pairs. Coordinates must be in WGS84 (EPSG:4326) format.
rectangle_vertices = [
(-74.02034270713835, 40.69992881162822), # Southwest corner (longitude, latitude)
(-74.02034270713835, 40.7111851828668), # Northwest corner (longitude, latitude)
(-74.00555129286164, 40.7111851828668), # Northeast corner (longitude, latitude)
(-74.00555129286164, 40.69992881162822) # Southeast corner (longitude, latitude)
]
Option 2: Draw Rectangle on Mapยถ
Interactive selection - click the rectangle tool and draw your area on the map.
#@title 2.2.1. Set target city (you need to fill in required values before executing!)
cityname = "new york" #@param {type:"string"}
#@title 2.2.2. Draw a rectangle on a map (Click โ on the left side, then click the north west corner and the south east corner of the rectangle.)
from voxcity.geoprocessor.draw import draw_rectangle_map_cityname
m, rectangle_vertices = draw_rectangle_map_cityname(cityname, zoom=15)
m
Option 3: Center Point + Dimensionsยถ
Define by clicking a center point and specifying width (east-west) and height (north-south) in meters.
#@title 2.3.1. Set width (m) and height (m) (you need to fill in required values before executing!)
width = 1250 #@param {type:"number"}
height = 1250 #@param {type:"number"}
#@title 2.3.2. Draw a center point (circlemarker) on a map (Click โ on the left side and then click the center location on the map.)
from voxcity.geoprocessor.draw import center_location_map_cityname
m, rectangle_vertices = center_location_map_cityname(cityname, width, height, zoom=15)
m
โ๏ธ Step 3: Configure Data Sourcesยถ
Select which data sources to use for buildings, terrain, land cover, and vegetation.
Available Data Sourcesยถ
Category |
Options |
|---|---|
Buildings |
OpenStreetMap, Overture, EUBUCCO v0.1, Open Building 2.5D, Microsoft Building Footprints |
DEM/Terrain |
USGS 3DEP 1m, DeltaDTM, FABDEM, NASA, Copernicus |
Land Cover |
OpenStreetMap, ESA WorldCover, ESRI 10m, Dynamic World |
Canopy Height |
High Resolution 1m, ETH Global Sentinel-2, Static |
#@title 3.1. Set data sources and meshsize (m)
building_source = 'OpenStreetMap' #@param ['OpenStreetMap', 'Overture', 'EUBUCCO v0.1', 'Open Building 2.5D Temporal', 'Microsoft Building Footprints', 'Local file']
building_complementary_source = "None" #@param ['None', 'Open Building 2.5D Temporal', 'Microsoft Building Footprints', 'England 1m DSM - DTM', 'Netherlands 0.5m DSM - DTM', 'OpenMapTiles', 'Local file', 'OpenStreetMap', 'Overture', 'EUBUCCO v0.1']
land_cover_source = 'OpenStreetMap' #@param ['OpenStreetMap', 'Urbanwatch', 'OpenEarthMapJapan', 'ESA WorldCover', 'ESRI 10m Annual Land Cover', 'Dynamic World V1']
canopy_height_source = 'High Resolution 1m Global Canopy Height Maps' #@param ['High Resolution 1m Global Canopy Height Maps', 'ETH Global Sentinel-2 10m Canopy Height (2020)', 'Static']
dem_source = 'USGS 3DEP 1m' #@param ['DeltaDTM', 'FABDEM', 'England 1m DTM', 'DEM France 1m', 'Netherlands 0.5m DTM', 'AUSTRALIA 5M DEM', 'USGS 3DEP 1m', 'NASA', 'COPERNICUS', 'Flat']
meshsize = 5 #@param {type:"number"}
#@title 3.2. Set optional parameters, data sources and meshsize (m)
kwargs = {
# "building_path": 'path_to_building_source_file', #To set path to building base data source when you select 'Local file' as building_source.
# "building_complementary_path": 'path_to_building_complemntary_source_file', #To set path to building complementary data source when you select 'Local file' as building_complementary_source.
"building_complementary_source": building_complementary_source, # # Specify the complementary building data source (e.g. 'Open Building 2.5D Temporal', 'None', etc.)
"building_complement_height": 10, # Default height in meters to use for buildings when height data is missing
# "complement_polygon": True, #Set 'True' if you want to incorporate building footprints from building complementary source. Only building footprints that do not have any intersections with footprints from building source are included.
"output_dir": 'output/test', #To set directory path for output files
# "remove_perimeter_object": 0.1, #Set value more than 0 if you want to remove objects including buildings and trees near to domain boundarys (perimeter). Used mainly for CFD simulation. For instance, when you set 0.1, objects with distances from domain boundaries less than 0.1 * domain width (height) are removed.
# "gridvis": True, #Set 'True' if you want to visualize extracted 2D grid data.
# "mapvis": False, #Set 'True' if you want to visualize extracted 2D grid data on a basemap. Note that it take longer time than "gridvis".
# "voxelvis": False, #Set 'True' if you want to visualize generated voxel 3d city model. Note that this visualiztion takes long time if the number of voxels is huge, e.g., more than one million.
# "voxelvis_img_save_path": None, #Set path to save image file of generated voxel 3d city model.
# "trunk_height_ratio": None, #To set ratio of tree trunk height against tree canopy height. Default: 0.59 (11.76 / 19.98).
# "min_canopy_height": None, #To set minimum canopy height in meters if you want to exclude trees lower than that height.
"dem_interpolation": True, #Set 'True' when mesh size if finer than resolution of dem data source and if you want to use interporation.
# "dynamic_world_date": '2021-04-02', #To set date of Dynamic World.
# "esri_landcover_year": '2023' #To set year of Esri Land Cover.
}
๐๏ธ Step 4: Generate 3D Voxel Modelยถ
This step downloads data from selected sources and creates the 3D voxel city model.
Note: Generation time depends on:
Area size
Voxel resolution (meshsize)
Data source availability and response times
#@title 4.1. Obtain grid data integrate it to create voxel data
from voxcity.generator import get_voxcity
city = get_voxcity(
rectangle_vertices,
meshsize=meshsize,
building_source=building_source,
land_cover_source=land_cover_source,
canopy_height_source=canopy_height_source,
dem_source=dem_source,
**kwargs
)
# Access grids from the VoxCity object
voxcity_grid = city.voxels.classes
building_height_grid = city.buildings.heights
building_min_height_grid = city.buildings.min_heights
building_id_grid = city.buildings.ids
land_cover_grid = city.land_cover.classes
dem_grid = city.dem.elevation
canopy_height_grid = city.tree_canopy.top
building_gdf = city.extras.get('building_gdf', None)
#@title 4.2. Visualize the generated 3D city model (This process may take a long time, so skip this step if you set a large target area)
from voxcity.visualizer import visualize_voxcity
visualize_voxcity(city, mode="static")
๐พ Step 5: Export to Different Formatsยถ
Export your 3D city model to various formats for use in other applications.
Format |
Use Case |
|---|---|
INX |
ENVI-met CFD simulations |
VOX |
MagicaVoxel visualization |
OBJ |
Blender, Rhino, 3D modeling |
#@title 5.1. Export INX file for ENVI-MET
from voxcity.exporter.envimet import export_inx, generate_edb_file
envimet_kwargs = {
"author_name": "enter your name", # Optional. To set author name in INX.
"model_description": "generated and exported using VoxCity", # Optional. To desctibe model in INX.
"domain_building_max_height_ratio": 2, # Optional. To set ratio between domain height (Z) and maximum height (building + terrain). Default: 2.0.
"useTelescoping_grid": True, # Optional. To activate telescoping grid. Default: False.
"verticalStretch": 20, # Optional. To set vertical stretch (%). Default: 0%.
"min_grids_Z": 20, # Optional. To set minimum number of vertical grid cells (Z-axis). Default: 20.
}
export_inx(
city,
output_directory='output/envimet',
file_basename='voxcity',
land_cover_source=land_cover_source,
**envimet_kwargs
)
generate_edb_file(lad=1.0)
#@title 5.2. VOX file for MagicaVoxel
from voxcity.exporter.magicavoxel import export_magicavoxel_vox
output_path = f"output/magicavoxel"
export_magicavoxel_vox(city, output_path)
#@title 5.3. OBJ file
from voxcity.exporter.obj import export_obj
output_directory = './output/obj'
output_file_name = 'voxcity'
export_obj(city, output_directory, output_file_name)
Step6 Urban simulationsยถ
6.1. Solar radiationยถ
#@title 6.1.1. Solar irradiance on building surfaces (instantaneous)
from voxcity.simulator.solar import get_building_global_solar_irradiance_using_epw
from voxcity.visualizer import visualize_voxcity
# Define kwargs dictionary for solar irradiance calculation
irradiance_kwargs = {
"calc_type": "instantaneous", # For single time point calculation
"download_nearest_epw": True,
"rectangle_vertices": rectangle_vertices,
# "epw_file_path": epw_file_path,
"calc_time": "01-01 12:00:00", # Solar noon on June 21 (summer solstice)
}
# Example 1: Instantaneous calculation
instantaneous_irradiance = get_building_global_solar_irradiance_using_epw(
city,
**irradiance_kwargs
)
# Visualize instantaneous results with building mesh overlay
visualize_voxcity(
city,
mode="static",
building_sim_mesh=instantaneous_irradiance,
building_value_name="global",
building_colormap="magma",
building_vmin=0,
voxel_color_map="grayscale"
)
#@title 6.1.2. Solar irradiance on building surfaces (cumulative)
from voxcity.simulator.solar import get_building_global_solar_irradiance_using_epw
from voxcity.visualizer import visualize_voxcity
# Define kwargs dictionary for cumulative calculation
cumulative_kwargs = {
"calc_type": "cumulative", # For cumulative calculations over a period
"download_nearest_epw": True,
"rectangle_vertices": rectangle_vertices,
# "epw_file_path": epw_file_path,
"period_start": "01-01 07:00:00", # June 1st at 7 AM
"period_end": "01-31 19:00:00", # June 30th at 7 PM
}
# Example 2: Cumulative calculation for a month
cumulative_irradiance = get_building_global_solar_irradiance_using_epw(
city,
**cumulative_kwargs
)
# Visualize cumulative results with building mesh overlay
visualize_voxcity(
city,
mode="static",
building_sim_mesh=cumulative_irradiance,
building_value_name="global",
building_colormap="magma",
building_vmin=0,
voxel_color_map="grayscale",
projection_type="perspective",
distance_factor=0.7
)
#@title 6.1.3. An example post analysis for solar irradiance on building surfaces (cumulative)
import pandas as pd
from shapely.geometry import box
from voxcity.visualizer import visualize_numerical_gdf_on_basemap
# Create dataframe with the three columns
cumulative_solar_df = pd.DataFrame({
'area': cumulative_irradiance.area_faces,
'cumulative_solar_irradiance': cumulative_irradiance.metadata['global'],
'building_id': cumulative_irradiance.metadata['building_id']
})
# Filter out rows where cumulative_solar_irradiance is NaN
cumulative_solar_df = cumulative_solar_df.dropna(subset=['cumulative_solar_irradiance'])
# Calculate mean solar irradiance per building
building_solar = cumulative_solar_df.groupby('building_id')['cumulative_solar_irradiance'].mean().reset_index()
# Merge building_gdf with solar data
# First merge solar data
building_gdf_with_solar = building_gdf.merge(
building_solar,
left_on='id',
right_on='building_id',
how='left'
)
# Create a bounding box from rectangle vertices and filter buildings within it
bbox = box(
rectangle_vertices[0][0], # minx
rectangle_vertices[0][1], # miny
rectangle_vertices[2][0], # maxx
rectangle_vertices[1][1] # maxy
)
building_gdf_with_solar = building_gdf_with_solar[building_gdf_with_solar.geometry.within(bbox)]
# Visualize using the utility function
visualize_numerical_gdf_on_basemap(
building_gdf_with_solar,
value_name='cumulative_solar_irradiance',
cmap='magma',
figsize=(20,8),
basemap='CartoDB light',
show_edge=True,
edge_color='black',
edge_width=0.5
)
# Print summary statistics
print("\nSummary Statistics for Building Solar Irradiance:")
print("------------------------------------------------")
print(building_gdf_with_solar['cumulative_solar_irradiance'].describe())
#@title 6.1.4. Ground-level solar irradiance (instantaneous)
from voxcity.simulator.solar import get_global_solar_irradiance_using_epw
solar_kwargs = {
"download_nearest_epw": True, # Whether to automatically download nearest EPW weather file based on location from Climate.OneBuilding.Org
"rectangle_vertices": rectangle_vertices, # Coordinates defining the area of interest for calculation
# "epw_file_path": "./output/new.york-downtown.manhattan.heli_ny_usa_1.epw", # Path to EnergyPlus Weather (EPW) file containing climate data. Set if you already have an EPW file.
"calc_time": "01-01 12:00:00", # Time for instantaneous calculation in format "MM-DD HH:MM:SS"
"view_point_height": 1.5, # Height of view point in meters for calculating solar access. Default: 1.5 m
"tree_k": 0.6, # Static extinction coefficient - controls how much sunlight is blocked by trees (higher = more blocking)
"tree_lad": 0.5, # Leaf area density of trees - density of leaves/branches that affect shading (higher = denser foliage)
"colormap": 'magma', # Matplotlib colormap for visualization. Default: 'viridis'
"obj_export": True, # Whether to export results as 3D OBJ file
"output_directory": 'output/test', # Directory for saving output files
"output_file_name": 'instantaneous_solar_irradiance', # Base filename for outputs (without extension)
"alpha": 1.0, # Transparency of visualization (0.0-1.0)
"vmin": 0, # Minimum value for colormap scaling in visualization
# "vmax": 900, # Maximum value for colormap scaling in visualization
}
# Compute global solar irradiance map (direct + diffuse radiation)
solar_grid = get_global_solar_irradiance_using_epw(
city, # VoxCity object representing the urban environment
calc_type='instantaneous', # Calculate instantaneous irradiance at specified time
direct_normal_irradiance_scaling=1.0, # Scaling factor for direct solar radiation (1.0 = no scaling)
diffuse_irradiance_scaling=1.0, # Scaling factor for diffuse solar radiation (1.0 = no scaling)
**solar_kwargs # Pass all the parameters defined above
)
#@title 6.1.5. Ground-level solar irradiance (cumulative)
from voxcity.simulator.solar import get_global_solar_irradiance_using_epw
# Dictionary containing parameters for solar irradiance calculation
solar_kwargs = {
"download_nearest_epw": True, # Whether to automatically download nearest EPW weather file based on location from Climate.OneBuilding.Org
"rectangle_vertices": rectangle_vertices, # Coordinates defining the area of interest for calculation
# "epw_file_path": "./output/USA_CA_Marina.Muni.AP.690070_TMYx.epw", # Path to EnergyPlus Weather (EPW) file containing climate data. Set if you already have an EPW file.
"start_time": "01-01 01:00:00", # Start time for cumulative calculation
"end_time": "01-31 23:00:00", # End time for cumulative calculation
"view_point_height": 1.5, # Height of view point in meters for calculating solar access. Default: 1.5 m
"tree_k": 0.6, # Static extinction coefficient - controls how much sunlight is blocked by trees (higher = more blocking)
"tree_lad": 0.5, # Leaf area density of trees - density of leaves/branches that affect shading (higher = denser foliage)
"colormap": 'magma', # Matplotlib colormap for visualization. Default: 'viridis'
"obj_export": True, # Whether to export results as 3D OBJ file
"output_directory": 'output/test', # Directory for saving output files
"output_file_name": 'solar_irradiance', # Base filename for outputs (without extension)
"alpha": 1.0, # Transparency of visualization (0.0-1.0)
"vmin": 0, # Minimum value for colormap scaling in visualization (commented out to auto-scale)
# "vmax": 900, # Maximum value for colormap scaling in visualization (commented out to auto-scale)
}
# Calculate cumulative solar irradiance over the specified time period
cum_solar_grid = get_global_solar_irradiance_using_epw(
city, # VoxCity object representing the urban environment
calc_type='cumulative', # Calculate cumulative irradiance over time period instead of instantaneous
direct_normal_irradiance_scaling=1.0, # Scaling factor for direct solar radiation (1.0 = no scaling)
diffuse_irradiance_scaling=1.0, # Scaling factor for diffuse solar radiation (1.0 = no scaling)
**solar_kwargs # Pass all the parameters defined above
)
#6.2. View index
#@title 6.2.1. Green view index, Sky view index
from voxcity.simulator.view import get_view_index
view_kwargs = {
"view_point_height": 1.5, # Height of observer viewpoint in meters
"tree_k": 0.6, # Static extinction coefficient - controls how much sunlight is blocked by trees (higher = more blocking)
"tree_lad": 1.0, # Leaf area density of trees - density of leaves/branches that affect shading (higher = denser foliage)
"colormap": "viridis", # Colormap for visualization
"obj_export": True, # Whether to export as OBJ file
"output_directory": "output", # Directory to save output files
"output_file_name": "gvi" # Base filename for outputs
}
# Compute Green View Index using mode='green'
gvi_grid = get_view_index(city, mode='green', **view_kwargs)
# Adjust parameters for Sky View Index
view_kwargs["colormap"] = "BuPu_r"
view_kwargs["output_file_name"] = "svi"
view_kwargs["elevation_min_degrees"] = 0
# Compute Sky View Index using mode='sky'
svi_grid = get_view_index(city, mode='sky', **view_kwargs)
6.3. Landmark visibilityยถ
#@title 6.3.1. Map for drawing a polygon to specify landmark buildings - draw a polygon that includes all desired landmark building footprints
from voxcity.geoprocessor.draw import display_buildings_and_draw_polygon
# Create interactive map with buildings and get drawn polygon
m_landmark, landmark_polygon = display_buildings_and_draw_polygon(building_gdf)
m_landmark
#@title 6.3.2. Simulate landmark visibility
from voxcity.simulator.view import get_landmark_visibility_map # Import function to analyze landmark visibility from different viewpoints
landmark_kwargs = {
"view_point_height": 1.5, # To set height of view point in meters. Default: 1.5 m.
"landmark_polygon": landmark_polygon,
"colormap": 'cool', # Choose a colormap. Default: 'viridis'.
"obj_export": True, # Set "True" if you export the result in an OBJ file.
"output_directory": 'output/obj', # To set directory path for output files. Default: False.
"output_file_name": 'landmark_visibility', # To set file name excluding extension. Default: 'view_index.
"alpha": 1.0, # Set transparency (0.0 to 1.0)
"vmin": 0.0, # Minimum value for colormap normalization
"vmax": 1.0 # Maximum value for colormap normalization
}
landmark_vis_map, voxcity_grid_landmark = get_landmark_visibility_map( # Calculate visibility map for landmarks
city, # VoxCity object representing the 3D urban environment
building_gdf=building_gdf, # GeoJSON containing building footprint data
**landmark_kwargs # Pass all configuration parameters defined above
)
#@title Data aggregation by road networks
from voxcity.geoprocessor.network import get_network_values
network_kwargs = {
"network_type": "walk",
"colormap": "magma",
"vis_graph": True,
"vmin": 0.0,
# "vmax": 65000,
"edge_width": 2,
"alpha": 0.8,
"zoom": 16
}
G, solar_edge_gdf = get_network_values(
cum_solar_grid,
rectangle_vertices,
meshsize,
value_name='Cumulative Global Solar Irradiance (W/mยฒยทhour), yearly, 7:00-9:00',
**network_kwargs
)
#@title Visualize a simulation result on a 3D city model
from voxcity.visualizer import visualize_voxcity
vis_3d_kwargs = {
"ground_sim_grid": cum_solar_grid,
"ground_dem_grid": dem_grid,
"ground_colormap": 'magma',
"ground_view_point_height": 1.5,
"ground_vmin": 0.0,
# "ground_vmax": 1500000,
"projection_type": "orthographic",
"distance_factor": 1.5
}
visualize_voxcity(
city,
mode="static",
output_directory='output/3d_vis',
**vis_3d_kwargs
)
from voxcity.visualizer import visualize_voxcity
vis_3d_kwargs = {
"ground_sim_grid": landmark_vis_map,
"ground_dem_grid": dem_grid,
"ground_colormap": 'cool',
"ground_view_point_height": 1.5,
"ground_vmin": 0.0,
"ground_vmax": 1.0,
"projection_type": "perspective",
"distance_factor": 0.5
}
visualize_voxcity(
city,
mode="static",
output_directory='output/3d_vis',
**vis_3d_kwargs
)