๐Ÿ™๏ธ VoxCity Complete Demoยถ

VoxCity is a Python package for generating 3D voxel-based urban models from open geospatial data.

What This Notebook Coversยถ

  1. Environment Setup - Install VoxCity and authenticate Google Earth Engine

  2. Area Selection - Define your study area using coordinates or interactive map

  3. Configuration - Select data sources and parameters

  4. Model Generation - Create 3D voxel city models

  5. 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
)