Deltares / ra2ce

RA2CE helps to quantify resilience of critical infrastructure networks, prioritize interventions and adaptation measures and select the most appropriate action perspective to increase resilience considering future conditions.
https://deltares.github.io/ra2ce/
Other
9 stars 2 forks source link

Downloading a network from a polygon crashes if the polygon is in a different projection (not in EPSG:4326) #458

Open Cham8920 opened 2 months ago

Cham8920 commented 2 months ago

To DO:

Ra2ce version checks

Reproducible example

Network from a shapefile which is an OSM-download

_network_section = NetworkSection(
    network_type=NetworkTypeEnum.DRIVE,
    road_types=[RoadTypeEnum.MOTORWAY,
                RoadTypeEnum.MOTORWAY_LINK,
                RoadTypeEnum.TRUNK, 
                RoadTypeEnum.TRUNK_LINK,
                RoadTypeEnum.PRIMARY, 
                RoadTypeEnum.PRIMARY_LINK,
                RoadTypeEnum.SECONDARY,
                RoadTypeEnum.SECONDARY_LINK,
                RoadTypeEnum.TERTIARY,
                RoadTypeEnum.TERTIARY_LINK,
                RoadTypeEnum.ROAD
                ], 
)

#pass the specified sections as arguments for configuration
_network_config_data = NetworkConfigData(
    root_path= root_dir,
    output_path= root_dir/"output",
    network=_network_section,
    static_path=root_dir.joinpath('static'),
    )

#download the network based on the polygon extent and the specified road characteristics
_graph,_gdf = OsmNetworkWrapper.get_network_from_polygon(_network_config_data, bbox_polygon)

Here the bbox coordinates are not in EPSG:4326 - The notebook crashes

Current behaviour

When the polygon is in another projection the notebook crashes and throws this error:

File ~\RA2CE\RA2CE_Rep\ra2ce\ra2ce\network\network_wrappers\osm_network_wrapper\osm_network_wrapper.py:133, in OsmNetworkWrapper.get_network_from_geojson(config_data)
    [125](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:125) if (
    [126](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:126)     not isinstance(config_data.network.polygon, Path)
    [127](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:127)     or not config_data.network.polygon.is_file()
    [128](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:128) ):
    [129](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:129)     raise ValueError(
    [130](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:130)         "A valid network polygon (.geojson) file path needs to be provided."
    [131](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:131)     )
--> [133](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:133) return OsmNetworkWrapper(config_data).get_network()

File ~\RA2CE\RA2CE_Rep\ra2ce\ra2ce\network\network_wrappers\osm_network_wrapper\osm_network_wrapper.py:64, in OsmNetworkWrapper.__init__(self, config_data)
     [62](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:62) self.network_type = config_data.network.network_type
     [63](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/network_wrappers/osm_network_wrapper/osm_network_wrapper.py:63) self.road_types = config_data.network.road_types
...
   [1419](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/networks_utils.py:1419)         raise ValueError(
   [1420](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/networks_utils.py:1420)             "Latitude is out of bounds, check your JSON format or data. The Coordinate Reference System should be in EPSG:4326."
   [1421](/ra2ce/examples/~/RA2CE/RA2CE_Rep/ra2ce/ra2ce/network/networks_utils.py:1421)         )

ValueError: Longitude is out of bounds, check your JSON format or data. The Coordinate Reference System should be in EPSG:4326.

Desired behaviour

01) Here the ability to specify a crs and the polygon gets reprojected through RA2CE to download the network through osm

class NetworkConfigData(ConfigDataProtocol):
    root_path: Optional[Path] = None
    input_path: Optional[Path] = None
    output_path: Optional[Path] = None
    static_path: Optional[Path] = None
    # CRS is not yet supported in the ini file, it might be relocated to a subsection.
    **crs: CRS = field(default_factory=lambda: CRS.from_user_input(4326))**
    project: ProjectSection = field(default_factory=ProjectSection)
    network: NetworkSection = field(default_factory=NetworkSection)
    origins_destinations: OriginsDestinationsSection = field(
        default_factory=OriginsDestinationsSection
    )

Additional context

This becomes an issue when

For b, I reproject the hazard map if its not in EPSG:4326 so I can get the correct projection: (If there's a better workflow that would also be great)

hazard_map = Path(hazard_map)
output_map = hazard_map.parent / f'reprojected_{hazard_map.stem}.tif'

# Open the source raster file
with rasterio.open(hazard_map) as src:
    src_crs = src.crs
    target_crs = 'EPSG:4326'

    # Check if the source CRS is already EPSG:4326
    if src_crs.to_string() == target_crs:
        # If the source raster is already in the target CRS, copy the file
        shutil.copy(hazard_map, output_map)
        print(f"The source raster is already in {target_crs}. The file has been copied to {output_map}.")
    else:
        src_dtype = src.dtypes[0]

        transform, width, height = calculate_default_transform(
            src_crs, target_crs, src.width, src.height, *src.bounds)

        with rasterio.open(output_map,
                           'w',
                           driver='GTiff',
                           height=height,
                           width=width,
                           count=src.count,
                           dtype=src_dtype,
                           crs=target_crs,
                           transform=transform) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src_crs,
                    dst_transform=transform,
                    dst_crs=target_crs,
                    resampling=Resampling.nearest)
        print(f"Reprojection complete. Output saved to {output_map}")