From 8d750badc30552a133dad7689ee4700fc6fd9803 Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 26 Sep 2023 15:12:22 +0200 Subject: [PATCH 01/12] Add new exact indicators at grid scale --- .../geoindicators/GenericIndicators.groovy | 5 ++++ .../geoindicators/RsuIndicators.groovy | 3 ++- .../WorkflowGeoIndicators.groovy | 24 +++++++++++++++---- .../geoindicators/SpatialUnitsTests.groovy | 1 + .../geoclimate/osm/WorkflowOSM.groovy | 2 +- .../geoclimate/osm/WorflowOSMTest.groovy | 1 + 6 files changed, 30 insertions(+), 6 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy index 41c28a5d5e..fef695a134 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy @@ -455,6 +455,11 @@ String distributionCharacterization(JdbcDataSource datasource, String distribTab def distribColumns = allColumns.minus(inputId.toUpperCase()) def nbDistCol = distribColumns.size() + if(distribColumns.size==0){ + error("Any columns to compute the distribution from the table $distribTableName".toString()) + return + } + def idxExtrem = nbDistCol - 1 def idxExtrem_1 = nbDistCol - 2 if (extremum.toUpperCase() == "LEAST") { diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index da1c592f7a..d0a976efba 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -838,7 +838,8 @@ String roofAreaDistribution(JdbcDataSource datasource, String rsu, String buildi String effectiveTerrainRoughnessLength(JdbcDataSource datasource, String rsuTable, String projectedFacadeAreaName, String geometricMeanBuildingHeightName, - List listLayersBottom = [0, 10, 20, 30, 40, 50], int numberOfDirection = 12, String prefixName) { + List listLayersBottom = [0, 10, 20, 30, 40, 50], + int numberOfDirection = 12, String prefixName) { def GEOMETRIC_COLUMN = "the_geom" def ID_COLUMN_RSU = "id_rsu" def BASE_NAME = "effective_terrain_roughness_length" diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index 53ddb2ddc4..5627839429 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -357,7 +357,8 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, [facadeDensListLayersBottom : [0, 10, 20, 30, 40, 50], facadeDensNumberOfDirection : 12, svfPointDensity : 0.008, - svfRayLength : 100, svfNumberOfDirection: 60, + svfRayLength : 100, + svfNumberOfDirection : 60, heightColumnName : "height_roof", inputFields : ["id_build", "the_geom"], levelForRoads : [0], @@ -1724,7 +1725,9 @@ String rasterizeIndicators(JdbcDataSource datasource, def unweightedBuildingIndicators = [:] def weightedBuildingIndicators = [:] list_indicators.each { - if (it.equalsIgnoreCase("BUILDING_FRACTION") || it.equalsIgnoreCase("BUILDING_SURFACE_DENSITY")) { + if (it.equalsIgnoreCase("BUILDING_FRACTION") + || it.equalsIgnoreCase("BUILDING_SURFACE_DENSITY") || + it.equalsIgnoreCase("ASPECT_RATIO")) { columnFractionsList.put(priorities.indexOf("building"), "building") } else if (it.equalsIgnoreCase("WATER_FRACTION")) { columnFractionsList.put(priorities.indexOf("water"), "water") @@ -1831,7 +1834,9 @@ String rasterizeIndicators(JdbcDataSource datasource, } } - if ((list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY") && building) || + def freeFacadeDensityExact + + if ((list_indicators*.toUpperCase().containsAll(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO"]) && building) || (list_indicators*.toUpperCase().contains("BUILDING_SURFACE_DENSITY") && building)) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings @@ -1843,7 +1848,7 @@ String rasterizeIndicators(JdbcDataSource datasource, return } } - def freeFacadeDensityExact = Geoindicators.RsuIndicators.freeExternalFacadeDensityExact(datasource, createScalesRelationsGridBl, + freeFacadeDensityExact = Geoindicators.RsuIndicators.freeExternalFacadeDensityExact(datasource, createScalesRelationsGridBl, grid, grid_column_identifier, prefixName) if (freeFacadeDensityExact) { if (list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY")) { @@ -1866,6 +1871,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } } + if (list_indicators*.toUpperCase().contains("BUILDING_HEIGHT_DIST") && building) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings @@ -1958,6 +1964,15 @@ String rasterizeIndicators(JdbcDataSource datasource, info "Cannot merge all indicators in grid table $grid_indicators_table." return } + + //Compute the aspect_ratio + if (list_indicators*.toUpperCase().contains("ASPECT_RATIO") && building){ + //Because all other indicators have been yet computed we just alter the table with the aspect_ratio column + datasource.execute(""" + ALTER TABLE $grid_indicators_table ADD COLUMN ASPECT_RATIO DOUBLE PRECISION; + UPDATE $grid_indicators_table SET ASPECT_RATIO = free_external_facade_density / (1 - building_fraction) + """.toString()) + } datasource.dropTable(indicatorTablesToJoin.keySet().toArray(new String[0])) // Remove temporary tables @@ -1972,6 +1987,7 @@ String rasterizeIndicators(JdbcDataSource datasource, * @param zoneEnvelopeTableName * @param x_size x size of the grid * @param y_size y size of the grid + * @param srid used to reproject the grid * @param outputTableName the name of grid table * @return */ diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy index fb0248aea4..69f301f2b3 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -177,6 +177,7 @@ class SpatialUnitsTests { assert h2GIS.getTable(pPointOnSurface).getRowCount() == 2 } + @Test void prepareGeometriesForRSUWithFilterTest() { h2GIS.load(SpatialUnitsTests.class.class.getResource("road_test.geojson"), true) diff --git a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy index 63629f5ad8..8221aacbaf 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy @@ -445,7 +445,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d def extract = OSMTools.Loader.extract(query) if (extract) { - Geometry geomArea = h2gis_datasource.firstRow("select st_extent(the_geom) as the_geom from ${zones.zone}").the_geom + Geometry geomArea = h2gis_datasource.firstRow("select st_extent(the_geom) as the_geom from ${zones.zone}".toString()).the_geom geomArea.setSRID(srid) Map gisLayersResults = OSM.InputDataLoading.createGISLayers(h2gis_datasource, extract, geomArea, srid) if (gisLayersResults) { diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy index 682332af30..613a0d3a55 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -652,6 +652,7 @@ class WorflowOSMTest extends WorkflowAbstractTest { dirFile.delete() dirFile.mkdir() + //def nominatim = OSMTools.Utilities.getNominatimData("Lorient") def osm_parmeters = [ From 73c024efb52571dfe3b3f47df2db3bb9ac5abdbb Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 26 Sep 2023 16:36:00 +0200 Subject: [PATCH 02/12] Add SVF_FRACTION on grid --- .../bdtopo/AbstractBDTopoWorkflow.groovy | 3 +- .../geoclimate/geoindicators/DataUtils.groovy | 2 +- .../geoindicators/RsuIndicators.groovy | 6 ++-- .../WorkflowGeoIndicators.groovy | 28 +++++++++++++++++-- .../geoindicators/RsuIndicatorsTests.groovy | 2 +- .../geoclimate/osm/WorkflowOSM.groovy | 6 ++-- .../geoclimate/osm/WorflowOSMTest.groovy | 14 +++++----- 7 files changed, 44 insertions(+), 17 deletions(-) diff --git a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy index 5fbac7e8f1..16f1f58bda 100644 --- a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy +++ b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy @@ -542,7 +542,8 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { } def allowed_grid_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION", "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION", "UTRF_FLOOR_AREA_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY", - "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION"] + "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", + "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF_FRACTION"] def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase()) if (allowedOutputIndicators) { //Update the RSU indicators list according the grid indicators diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy index 9e046f840f..66b73dcbf8 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy @@ -73,7 +73,7 @@ String joinTables(JdbcDataSource datasource, Map inputTableNamesWithId, String o } leftQuery += " LEFT JOIN $key as $alias ON $alias.$value = $columnKey " } - indexes += "CREATE INDEX IF NOT EXISTS ${key}_ids ON $key USING BTREE($value);" + indexes += "CREATE INDEX IF NOT EXISTS ${key}_ids ON $key ($value);" alias++ } diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index d0a976efba..25865a6955 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -225,11 +225,11 @@ String freeExternalFacadeDensityExact(JdbcDataSource datasource, String building * @author Jérémy Bernard * @author Erwan Bocher */ -String groundSkyViewFactor(JdbcDataSource datasource, String rsu, String correlationBuildingTable, float pointDensity, +String groundSkyViewFactor(JdbcDataSource datasource, String rsu, String id_rsu, String correlationBuildingTable, float pointDensity, float rayLength, int numberOfDirection, String prefixName) { def GEOMETRIC_COLUMN_RSU = "the_geom" def GEOMETRIC_COLUMN_BU = "the_geom" - def ID_COLUMN_RSU = "id_rsu" + def ID_COLUMN_RSU = id_rsu def HEIGHT_WALL = "height_wall" def BASE_NAME = "ground_sky_view_factor" @@ -318,7 +318,7 @@ String groundSkyViewFactor(JdbcDataSource datasource, String rsu, String correla // The result of the SVF calculation is averaged at RSU scale datasource """ - CREATE TABLE $outputTableName(id_rsu integer, $BASE_NAME double) + CREATE TABLE $outputTableName($ID_COLUMN_RSU integer, $BASE_NAME double) AS (SELECT a.$ID_COLUMN_RSU, CASE WHEN AVG(b.SVF) is not null THEN AVG(b.SVF) ELSE 1 END FROM $rsu a LEFT JOIN $svfPts b diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index 5627839429..c626982fe0 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -692,7 +692,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, AS SELECT 1-extended_free_facade_fraction AS GROUND_SKY_VIEW_FACTOR, $columnIdRsu FROM ${computeExtFF}; DROP TABLE ${computeExtFF}""".toString() } else { - def computeSVF = Geoindicators.RsuIndicators.groundSkyViewFactor(datasource, intermediateJoinTable, + def computeSVF = Geoindicators.RsuIndicators.groundSkyViewFactor(datasource, intermediateJoinTable,"id_rsu", buildingTable, parameters.svfPointDensity, parameters.svfRayLength, parameters.svfNumberOfDirection, temporaryPrefName) @@ -1915,7 +1915,6 @@ String rasterizeIndicators(JdbcDataSource datasource, } } } - tablesToDrop< Date: Fri, 29 Sep 2023 18:36:36 +0200 Subject: [PATCH 03/12] Prepare HEIGHT_OF_ROUGHNESS_ELEMENTS at grid scale Improve some grid indicators to compute the exact part of geometry inside a grid cell Update Groovy to 3.0.19 because 3.0.11 is no longer available on the groovy download page --- .../bdtopo/AbstractBDTopoWorkflow.groovy | 2 +- .../geoindicators/RsuIndicators.groovy | 82 ++++++---- .../WorkflowGeoIndicators.groovy | 152 ++++++++++++++---- .../geoindicators/RsuIndicatorsTests.groovy | 10 +- .../geoclimate/osm/WorkflowOSM.groovy | 2 +- .../geoclimate/osm/WorflowOSMTest.groovy | 9 +- pom.xml | 2 +- 7 files changed, 184 insertions(+), 75 deletions(-) diff --git a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy index 16f1f58bda..c4be56987d 100644 --- a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy +++ b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy @@ -543,7 +543,7 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { def allowed_grid_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION", "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION", "UTRF_FLOOR_AREA_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY", "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", - "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF_FRACTION"] + "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS"] def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase()) if (allowedOutputIndicators) { //Update the RSU indicators list according the grid indicators diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index 25865a6955..b33512da61 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -207,6 +207,7 @@ String freeExternalFacadeDensityExact(JdbcDataSource datasource, String building * @param datasource A connexion to a database (H2GIS, PostGIS, ...) where are stored the input Table and in which * the resulting database will be stored * @param rsu The name of the input ITable where are stored the RSU + * @param id_rsu Unique identifier column name * @param correlationBuildingTable The name of the input ITable where are stored the buildings and the relationships * between buildings and RSU * @param pointDensity The density of points (nb / free m²) used to calculate the spatial average SVF. Use 0.008f @@ -392,13 +393,13 @@ String aspectRatio(JdbcDataSource datasource, String rsuTable, String rsuFreeExt * @return A database table name. * @author Jérémy Bernard */ -String projectedFacadeAreaDistribution(JdbcDataSource datasource, String building, String rsu, +String projectedFacadeAreaDistribution(JdbcDataSource datasource, String building, String rsu,String id_rsu, List listLayersBottom = [0, 10, 20, 30, 40, 50], int numberOfDirection = 12, String prefixName) { def BASE_NAME = "projected_facade_area_distribution" def GEOMETRIC_COLUMN_RSU = "the_geom" def GEOMETRIC_COLUMN_BU = "the_geom" - def ID_COLUMN_RSU = "id_rsu" + def ID_COLUMN_RSU = id_rsu def ID_COLUMN_BU = "id_build" def HEIGHT_WALL = "height_wall" @@ -410,7 +411,7 @@ String projectedFacadeAreaDistribution(JdbcDataSource datasource, String buildin datasource.createSpatialIndex(building,"the_geom") datasource.createSpatialIndex(rsu,"the_geom") datasource.createIndex(building,"id_build") - datasource.createIndex(rsu,"id_rsu") + datasource.createIndex(rsu,ID_COLUMN_RSU) if (360 % numberOfDirection == 0 && numberOfDirection % 2 == 0) { @@ -504,7 +505,7 @@ String projectedFacadeAreaDistribution(JdbcDataSource datasource, String buildin // Intersections between free facades and rsu geometries are calculated datasource """ DROP TABLE IF EXISTS $buildingFreeExpl; - CREATE TABLE $buildingFreeExpl(id_rsu INTEGER, the_geom GEOMETRY, $namesAndType) AS + CREATE TABLE $buildingFreeExpl($ID_COLUMN_RSU INTEGER, the_geom GEOMETRY, $namesAndType) AS (SELECT a.$ID_COLUMN_RSU, ST_INTERSECTION(a.$GEOMETRIC_COLUMN_RSU, ST_TOMULTILINE(b.the_geom)), ${onlyNamesB} FROM $rsu a, $buildingLayer b WHERE a.$GEOMETRIC_COLUMN_RSU && b.the_geom @@ -513,8 +514,8 @@ String projectedFacadeAreaDistribution(JdbcDataSource datasource, String buildin // Intersections facades are exploded to multisegments datasource """DROP TABLE IF EXISTS $rsuInter; - CREATE TABLE $rsuInter(id_rsu INTEGER, the_geom GEOMETRY, $namesAndType) - AS (SELECT id_rsu, the_geom, ${onlyNames} FROM ST_EXPLODE('$buildingFreeExpl'))""".toString() + CREATE TABLE $rsuInter($ID_COLUMN_RSU INTEGER, the_geom GEOMETRY, $namesAndType) + AS (SELECT $ID_COLUMN_RSU, the_geom, ${onlyNames} FROM ST_EXPLODE('$buildingFreeExpl'))""".toString() // The analysis is then performed for each direction ('numberOfDirection' / 2 because calculation @@ -550,8 +551,8 @@ String projectedFacadeAreaDistribution(JdbcDataSource datasource, String buildin sumNamesDir = sumNamesDir.join(",") def query = "DROP TABLE IF EXISTS $finalIndicator; " + - "CREATE TABLE $finalIndicator AS SELECT a.id_rsu," + queryColumns + - " FROM (SELECT id_rsu, CASE WHEN ST_AZIMUTH(ST_STARTPOINT(the_geom), ST_ENDPOINT(the_geom)) >= PI()" + + "CREATE TABLE $finalIndicator AS SELECT a.$ID_COLUMN_RSU," + queryColumns + + " FROM (SELECT $ID_COLUMN_RSU, CASE WHEN ST_AZIMUTH(ST_STARTPOINT(the_geom), ST_ENDPOINT(the_geom)) >= PI()" + "THEN ST_AZIMUTH(ST_STARTPOINT(the_geom), ST_ENDPOINT(the_geom)) - PI() " + "ELSE ST_AZIMUTH(ST_STARTPOINT(the_geom), ST_ENDPOINT(the_geom)) END AS azimuth," + " ST_LENGTH(the_geom) AS length, ${onlyNames} FROM $rsuInter) a" @@ -559,15 +560,15 @@ String projectedFacadeAreaDistribution(JdbcDataSource datasource, String buildin datasource query.toString() - datasource.createIndex(finalIndicator,"id_rsu") + datasource.createIndex(finalIndicator,ID_COLUMN_RSU) // Sum area at RSU scale and fill null values with 0 datasource """ DROP TABLE IF EXISTS $outputTableName; CREATE TABLE ${outputTableName} - AS SELECT a.id_rsu, ${sumNamesDir} + AS SELECT a.$ID_COLUMN_RSU, ${sumNamesDir} FROM $rsu a LEFT JOIN $finalIndicator b - ON a.id_rsu = b.id_rsu - GROUP BY a.id_rsu""".toString() + ON a.$ID_COLUMN_RSU = b.$ID_COLUMN_RSU + GROUP BY a.$ID_COLUMN_RSU""".toString() // Remove all temporary tables created datasource """DROP TABLE IF EXISTS $buildingIntersection, $buildingIntersectionExpl, @@ -820,6 +821,7 @@ String roofAreaDistribution(JdbcDataSource datasource, String rsu, String buildi * the resulting database will be stored * @param rsuTable the name of the input ITable where are stored the RSU geometries, the rsu_geometric_mean_height and * the rsu_projected_facade_area_distribution fields + * @param id_rsu Unique identifier column name * @param projectedFacadeAreaName the name base used for naming the projected facade area field within the * inputA rsuTable (default rsu_projected_facade_area_distribution - note that the field is also constructed * with the direction and vertical height informations) @@ -835,13 +837,13 @@ String roofAreaDistribution(JdbcDataSource datasource, String rsu, String buildi * * @author Jérémy Bernard */ -String effectiveTerrainRoughnessLength(JdbcDataSource datasource, String rsuTable, +String effectiveTerrainRoughnessLength(JdbcDataSource datasource, String rsuTable, String id_rsu, String projectedFacadeAreaName, String geometricMeanBuildingHeightName, List listLayersBottom = [0, 10, 20, 30, 40, 50], int numberOfDirection = 12, String prefixName) { def GEOMETRIC_COLUMN = "the_geom" - def ID_COLUMN_RSU = "id_rsu" + def ID_COLUMN_RSU = id_rsu def BASE_NAME = "effective_terrain_roughness_length" debug "Executing RSU effective terrain roughness length" @@ -1682,7 +1684,7 @@ String buildingSurfaceDensity(JdbcDataSource datasource, String facadeDensityTab } /** - * Script to compute the distributio of (only) horizontal roof area fraction for each layer of the canopy. If the height + * Script to compute the distribution of (only) horizontal roof area fraction for each layer of the canopy. If the height * roof and height wall differ, take the average value at building height. Note that this process first cut the buildings * according to RSU in order to calculate the exact distribution. * @@ -1695,6 +1697,7 @@ String buildingSurfaceDensity(JdbcDataSource datasource, String facadeDensityTab * @param prefixName String use as prefix to name the output table * @param listLayersBottom the list of height corresponding to the bottom of each vertical layers (default * [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]) + * @param cutBuilding set to true if the building layer must be cuted by the rsu layer * * @return Table name in which the rsu id and their corresponding indicator value are stored * @@ -1702,7 +1705,7 @@ String buildingSurfaceDensity(JdbcDataSource datasource, String facadeDensityTab */ String roofFractionDistributionExact(JdbcDataSource datasource, String rsu, String building, String idRsu, List listLayersBottom = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], - boolean density = true, String prefixName) { + boolean cutBuilding=true, String prefixName) { def GEOMETRIC_COLUMN_RSU = "the_geom" def GEOMETRIC_COLUMN_BU = "the_geom" def ID_COLUMN_BU = "id_build" @@ -1719,15 +1722,17 @@ String roofFractionDistributionExact(JdbcDataSource datasource, String rsu, Stri // To avoid overwriting the output files of this step, a unique identifier is created // Temporary table names - def buildInter = postfix "build_inter" + def buildInter =building def rsuBuildingArea = postfix "rsu_building_area" def buildFracH = postfix "roof_frac_H" def bufferTable = postfix "buffer_table" - // 1. Create the intersection between buildings and RSU polygons - datasource."$building"."$ID_COLUMN_BU".createIndex() - datasource."$rsu"."$idRsu".createIndex() - datasource """ + if(cutBuilding) { + buildInter = postfix "build_inter" + // 1. Create the intersection between buildings and RSU polygons + datasource."$building"."$ID_COLUMN_BU".createIndex() + datasource."$rsu"."$idRsu".createIndex() + datasource """ DROP TABLE IF EXISTS $buildInter; CREATE TABLE $buildInter AS SELECT a.$ID_COLUMN_BU, a.$idRsu, @@ -1735,9 +1740,10 @@ String roofFractionDistributionExact(JdbcDataSource datasource, String rsu, Stri (a.$HEIGHT_WALL + a.$HEIGHT_ROOF) / 2 AS $BUILD_HEIGHT FROM $building AS a LEFT JOIN $rsu AS b ON a.$idRsu = b.$idRsu""".toString() + } // 2. Calculate the total building roof area within each RSU - datasource."$buildInter"."$idRsu".createIndex() + datasource.createIndex(buildInter, idRsu) datasource """ DROP TABLE IF EXISTS $rsuBuildingArea; CREATE TABLE $rsuBuildingArea @@ -1746,9 +1752,8 @@ String roofFractionDistributionExact(JdbcDataSource datasource, String rsu, Stri GROUP BY $idRsu""".toString() // 3. Calculate the fraction of roof for each level of the canopy (defined by 'listLayersBottom') except the last - datasource."$buildInter"."$BUILD_HEIGHT".createIndex() - datasource."$buildInter"."$idRsu".createIndex() - datasource."$rsuBuildingArea"."$idRsu".createIndex() + datasource.createIndex(buildInter,BUILD_HEIGHT) + datasource.createIndex(rsuBuildingArea,idRsu) def tab_H = [:] def indicToJoin = [:] for (i in 1..(listLayersBottom.size() - 1)) { @@ -1841,14 +1846,17 @@ String roofFractionDistributionExact(JdbcDataSource datasource, String rsu, Stri * [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50]) * @param numberOfDirection the number of directions used for the calculation - according to the method used it should * be divisor of 360 AND a multiple of 2 (default 12) + * @param distributionAsIndex set the value to false to avoid normalizing the distribution by the area of the input rsu + * and the range of the bottom layers * @param prefixName String use as prefix to name the output table * * @return A database table name. * @author Jérémy Bernard + * TODO : https://github.com/orbisgis/geoclimate/issues/848 */ String frontalAreaIndexDistribution(JdbcDataSource datasource, String building, String rsu, String idRsu, List listLayersBottom = [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], - int numberOfDirection = 12, String prefixName) { + int numberOfDirection = 12, boolean distributionAsIndex = true, String prefixName) { def GEOMETRIC_FIELD_RSU = "the_geom" def GEOMETRIC_FIELD_BU = "the_geom" def ID_FIELD_BU = "id_build" @@ -1940,7 +1948,8 @@ String frontalAreaIndexDistribution(JdbcDataSource datasource, String building, // Indicator name def indicName = "FRONTAL_AREA_INDEX_H${layer_bottom}_${layer_top}_D${k * angleRangeDeg}_${(k + 1) * angleRangeDeg}" // Define query to sum the projected facade for buildings and shared facades - dirQueryVertFrac[k] = """ + if(distributionAsIndex) { + dirQueryVertFrac[k] = """ CASE WHEN $v > AZIMUTH AND $v-AZIMUTH < PI() THEN CASE WHEN $HEIGHT_WALL >= $layer_top THEN LENGTH*SIN($v-AZIMUTH) @@ -1954,7 +1963,24 @@ String frontalAreaIndexDistribution(JdbcDataSource datasource, String building, ELSE 0 END END AS $indicName""" - dirQueryDiv[k] = """COALESCE(SUM(b.$indicName)/ST_AREA(a.$GEOMETRIC_FIELD_RSU), 0) AS $indicName""" + dirQueryDiv[k] = """COALESCE(SUM(b.$indicName)/ST_AREA(a.$GEOMETRIC_FIELD_RSU), 0) AS $indicName""" + }else{ + dirQueryVertFrac[k] = """ + CASE WHEN $v > AZIMUTH AND $v-AZIMUTH < PI() + THEN CASE WHEN $HEIGHT_WALL >= $layer_top + THEN LENGTH*SIN($v-AZIMUTH) + ELSE LENGTH*SIN($v-AZIMUTH)*($HEIGHT_WALL-$layer_bottom) + END + ELSE CASE WHEN $v - AZIMUTH < -PI() + THEN CASE WHEN $HEIGHT_WALL >= $layer_top + THEN LENGTH*SIN($v+2*PI()-AZIMUTH) + ELSE LENGTH*SIN($v+2*PI()-AZIMUTH)*($HEIGHT_WALL-$layer_bottom) + END + ELSE 0 + END + END AS $indicName""" + dirQueryDiv[k] = """COALESCE(SUM(b.$indicName), 0) AS $indicName""" + } } // Calculates projected surfaces for buildings and shared facades datasource """ diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index c626982fe0..4906a63474 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -26,6 +26,8 @@ import org.locationtech.jts.geom.Geometry import org.orbisgis.data.jdbc.JdbcDataSource import org.orbisgis.geoclimate.Geoindicators +import java.sql.SQLException + @BaseScript Geoindicators geoindicators /** @@ -638,7 +640,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, facadeDensNumberOfDirection: 8 } def projFacadeDist = Geoindicators.RsuIndicators.projectedFacadeAreaDistribution(datasource, buildingTable, - rsu, facadeDensListLayersBottom, + rsu, "id_rsu", facadeDensListLayersBottom, facadeDensNumberOfDirection, temporaryPrefName) if (!projFacadeDist) { info "Cannot compute the projected facade distribution. " @@ -708,7 +710,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, // rsu_effective_terrain_roughness if (indicatorUse*.toUpperCase().contains("LCZ") || indicatorUse*.toUpperCase().contains("TEB")) { // Create the join tables to have all needed input fields for aspect ratio computation - def effRoughHeight = Geoindicators.RsuIndicators.effectiveTerrainRoughnessLength(datasource, intermediateJoinTable, + def effRoughHeight = Geoindicators.RsuIndicators.effectiveTerrainRoughnessLength(datasource, intermediateJoinTable,"id_rsu", "projected_facade_area_distribution", "geom_avg_$heightColumnName", facadeDensListLayersBottom, facadeDensNumberOfDirection, temporaryPrefName) @@ -722,7 +724,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, if (indicatorUse*.toUpperCase().contains("LCZ")) { def roughClass = Geoindicators.RsuIndicators.effectiveTerrainRoughnessClass(datasource, effRoughHeight, "effective_terrain_roughness_length", temporaryPrefName) if (!roughClass) { - info "Cannot compute the SVF calculation." + info "Cannot compute the roughness class." return } finalTablesToJoin.put(roughClass, columnIdRsu) @@ -1724,6 +1726,7 @@ String rasterizeIndicators(JdbcDataSource datasource, def unweightedBuildingIndicators = [:] def weightedBuildingIndicators = [:] + def height_roof_unweighted_list =[] list_indicators.each { if (it.equalsIgnoreCase("BUILDING_FRACTION") || it.equalsIgnoreCase("BUILDING_SURFACE_DENSITY") || @@ -1739,12 +1742,18 @@ String rasterizeIndicators(JdbcDataSource datasource, } else if (it.equalsIgnoreCase("IMPERVIOUS_FRACTION")) { columnFractionsList.put(priorities.indexOf("impervious"), "impervious") } else if (it.equalsIgnoreCase("BUILDING_HEIGHT") && building) { - unweightedBuildingIndicators.put("height_roof", ["AVG", "STD"]) + height_roof_unweighted_list.addAll( ["AVG", "STD"]) } else if (it.equalsIgnoreCase("BUILDING_HEIGHT_WEIGHTED") && building) { weightedBuildingIndicators["height_roof"] = ["area": ["AVG", "STD"]] } else if (it.equalsIgnoreCase("BUILDING_POP") && building) { unweightedBuildingIndicators.put("pop", ["SUM"]) } + else if (it.equalsIgnoreCase("HEIGHT_OF_ROUGHNESS_ELEMENTS") && building) { + height_roof_unweighted_list.add("GEOM_AVG") + } + } + if(height_roof_unweighted_list){ + unweightedBuildingIndicators.put("height_roof",height_roof_unweighted_list) } // Calculate all surface fractions indicators on the GRID cell @@ -1768,6 +1777,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } } String createScalesRelationsGridBl + String computeBuildingStats if (unweightedBuildingIndicators) { // Create the relations between grid cells and buildings createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, building, grid, @@ -1777,7 +1787,7 @@ String rasterizeIndicators(JdbcDataSource datasource, info "Cannot compute the scales relations between buildings and grid cells." return } - def computeBuildingStats = Geoindicators.GenericIndicators.unweightedOperationFromLowerScale(datasource, createScalesRelationsGridBl, + computeBuildingStats = Geoindicators.GenericIndicators.unweightedOperationFromLowerScale(datasource, createScalesRelationsGridBl, grid, grid_column_identifier, grid_column_identifier, unweightedBuildingIndicators, prefixName) @@ -1789,17 +1799,15 @@ String rasterizeIndicators(JdbcDataSource datasource, } + String buildingCutted if (weightedBuildingIndicators) { - if (!createScalesRelationsGridBl) { - // Create the relations between grid cells and buildings - createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, building, - grid, grid_column_identifier, null, prefixName) - if (!createScalesRelationsGridBl) { - info "Cannot compute the scales relations between buildings and grid cells." - return - } + //Cut the building to compute exact fractions + buildingCutted = cutBuilding(datasource, grid, building) + if(!buildingCutted){ + info "Cannot split the building with the grid to compute the weighted statistics" + return } - def computeWeightedAggregStat = Geoindicators.GenericIndicators.weightedAggregatedStatistics(datasource, createScalesRelationsGridBl, + def computeWeightedAggregStat = Geoindicators.GenericIndicators.weightedAggregatedStatistics(datasource, buildingCutted, grid, grid_column_identifier, weightedBuildingIndicators, prefixName) if (!computeWeightedAggregStat) { @@ -1811,20 +1819,17 @@ String rasterizeIndicators(JdbcDataSource datasource, } if (list_indicators*.toUpperCase().contains("BUILDING_TYPE") && building) { - if (!createScalesRelationsGridBl) { - // Create the relations between grid cells and buildings - createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, building, grid, - grid_column_identifier, null, - prefixName) - if (!createScalesRelationsGridBl) { - info "Cannot compute the scales relations between buildings and grid cells." + if (!buildingCutted) { + buildingCutted = cutBuilding(datasource, grid, building) + if(!buildingCutted){ + info "Cannot split the building with the grid to compute the weighted statistics" return } } def indicatorName = "TYPE" def upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics(datasource, grid, grid_column_identifier, - createScalesRelationsGridBl, + buildingCutted, indicatorName, indicatorName,false, "building_type_fraction") if (upperScaleAreaStatistics) { @@ -1873,19 +1878,16 @@ String rasterizeIndicators(JdbcDataSource datasource, if (list_indicators*.toUpperCase().contains("BUILDING_HEIGHT_DIST") && building) { - if (!createScalesRelationsGridBl) { - // Create the relations between grid cells and buildings - createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, building, - grid, grid_column_identifier, null, - prefixName) - if (!createScalesRelationsGridBl) { - info "Cannot compute the scales relations between buildings and grid cells." + if (!buildingCutted) { + buildingCutted = cutBuilding(datasource, grid, building) + if(!buildingCutted){ + info "Cannot split the building with the grid to compute the building height distance" return } } def roofFractionDistributionExact = Geoindicators.RsuIndicators.roofFractionDistributionExact(datasource, - grid, createScalesRelationsGridBl,grid_column_identifier, - [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], prefixName) + grid, buildingCutted,grid_column_identifier, + [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], false, prefixName) if (roofFractionDistributionExact) { indicatorTablesToJoin.put(roofFractionDistributionExact, grid_column_identifier) } else { @@ -1907,7 +1909,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } def frontalAreaIndexDistribution = Geoindicators.RsuIndicators.frontalAreaIndexDistribution(datasource, createScalesRelationsGridBl, grid, grid_column_identifier, - [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], 12, prefixName,) + [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], 12, true, prefixName) if (frontalAreaIndexDistribution) { indicatorTablesToJoin.put(frontalAreaIndexDistribution, grid_column_identifier) } else { @@ -1958,7 +1960,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if (list_indicators*.toUpperCase().contains("SVF_FRACTION") && building) { + if (list_indicators*.toUpperCase().contains("SVF") && building) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, building, @@ -1980,7 +1982,58 @@ String rasterizeIndicators(JdbcDataSource datasource, tablesToDrop<2.1.0-SNAPSHOT 2.0.9 1.4.11 - 3.0.11 + 3.0.19 4.6.3 2.11.0 1.21 From edccc5423d5aee349f2b41d122fd44089f561e6e Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 3 Oct 2023 17:32:30 +0200 Subject: [PATCH 04/12] Add HEIGHT_OF_ROUGHNESS_ELEMENTS at grid scale --- .../geoindicators/RsuIndicators.groovy | 13 ++-- .../WorkflowGeoIndicators.groovy | 68 ++++++++++--------- .../geoclimate/osm/WorflowOSMTest.groovy | 4 +- 3 files changed, 47 insertions(+), 38 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index b33512da61..ba151cb151 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -861,14 +861,16 @@ String effectiveTerrainRoughnessLength(JdbcDataSource datasource, String rsuTabl // The name of the outputTableName is constructed def outputTableName = prefix prefixName, "rsu_" + BASE_NAME + def layerSize = listLayersBottom.size() // The lambda_f indicator is first calculated def lambdaQuery = "DROP TABLE IF EXISTS $lambdaTable;" + "CREATE TABLE $lambdaTable AS SELECT $ID_COLUMN_RSU, $geometricMeanBuildingHeightName, (" - for (int i in 1..listLayersBottom.size()) { + for (int i in 1..layerSize) { + //TODO : why an array here and not a variable names[i - 1] = "${projectedFacadeAreaName}_H${listLayersBottom[i - 1]}_${listLayersBottom[i]}" - if (i == listLayersBottom.size()) { - names[listLayersBottom.size() - 1] = - "${projectedFacadeAreaName}_H${listLayersBottom[listLayersBottom.size() - 1]}" + if (i == layerSize) { + names[layerSize - 1] = + "${projectedFacadeAreaName}_H${listLayersBottom[layerSize - 1]}" } for (int d = 0; d < numberOfDirection / 2; d++) { def dirDeg = d * 360 / numberOfDirection @@ -2017,6 +2019,9 @@ String frontalAreaIndexDistribution(JdbcDataSource datasource, String building, dirList.each { k, v -> // Indicator name def indicName = "FRONTAL_AREA_INDEX_H${layer_bottom}_${layer_top}_D${k * angleRangeDeg}_${(k + 1) * angleRangeDeg}" + if(!distributionAsIndex){ + indicName = "FRONTAL_AREA_INDEX_H${layer_bottom}_D${k * angleRangeDeg}_${(k + 1) * angleRangeDeg}" + } // Define query to calculate the vertical fraction of projected facade for buildings and shared facades dirQueryVertFrac[k] = """ CASE WHEN $v-AZIMUTH > 0 AND $v-AZIMUTH < PI() diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index 4906a63474..79e947845a 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -1622,6 +1622,11 @@ String rasterizeIndicators(JdbcDataSource datasource, info "The list of indicator names cannot be null or empty" return } + + //Concert the list of indicators to upper case + + def list_indicators_upper = list_indicators.collect{ it.toUpperCase() } + def tablesToDrop=[] // Temporary (and output tables) are created def tesselatedSeaLandTab = postfix "TESSELATED_SEA_LAND" @@ -1635,7 +1640,7 @@ String rasterizeIndicators(JdbcDataSource datasource, /* * Make aggregation process with previous grid and current rsu lcz */ - if (list_indicators.findAll { it.toUpperCase() in ["LCZ_FRACTION", "LCZ_PRIMARY"] } && rsu_lcz) { + if (list_indicators_upper.intersect(["LCZ_FRACTION", "LCZ_PRIMARY"]) && rsu_lcz) { def indicatorName = "LCZ_PRIMARY" String upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics( datasource, grid, grid_column_identifier, @@ -1644,7 +1649,7 @@ String rasterizeIndicators(JdbcDataSource datasource, if (upperScaleAreaStatistics) { indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier) - if (list_indicators*.toUpperCase().contains("LCZ_PRIMARY")) { + if (list_indicators_upper.contains("LCZ_PRIMARY")) { def resultsDistrib = Geoindicators.GenericIndicators.distributionCharacterization(datasource, upperScaleAreaStatistics, upperScaleAreaStatistics, grid_column_identifier, ["equality", "uniqueness"], @@ -1693,7 +1698,7 @@ String rasterizeIndicators(JdbcDataSource datasource, /* * Make aggregation process with previous grid and current rsu urban typo area */ - if (list_indicators*.toUpperCase().contains("UTRF_AREA_FRACTION") && rsu_utrf_area) { + if (list_indicators_upper.contains("UTRF_AREA_FRACTION") && rsu_utrf_area) { String indicatorName = "TYPO_MAJ" String upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics(datasource, grid, grid_column_identifier, rsu_utrf_area, @@ -1705,7 +1710,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } } - if (list_indicators*.toUpperCase().contains("UTRF_FLOOR_AREA_FRACTION") && rsu_utrf_floor_area) { + if (list_indicators_upper.contains("UTRF_FLOOR_AREA_FRACTION") && rsu_utrf_floor_area) { def indicatorName = "TYPO_MAJ" def upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics(datasource, grid, grid_column_identifier, rsu_utrf_floor_area, indicatorName, "FLOOR_AREA_TYPO_MAJ", false, @@ -1727,28 +1732,28 @@ String rasterizeIndicators(JdbcDataSource datasource, def unweightedBuildingIndicators = [:] def weightedBuildingIndicators = [:] def height_roof_unweighted_list =[] - list_indicators.each { - if (it.equalsIgnoreCase("BUILDING_FRACTION") - || it.equalsIgnoreCase("BUILDING_SURFACE_DENSITY") || - it.equalsIgnoreCase("ASPECT_RATIO")) { + list_indicators_upper.each { + if (it == "BUILDING_FRACTION" + || it=="BUILDING_SURFACE_DENSITY" || + it=="ASPECT_RATIO") { columnFractionsList.put(priorities.indexOf("building"), "building") - } else if (it.equalsIgnoreCase("WATER_FRACTION")) { + } else if (it=="WATER_FRACTION") { columnFractionsList.put(priorities.indexOf("water"), "water") - } else if (it.equalsIgnoreCase("VEGETATION_FRACTION")) { + } else if (it=="VEGETATION_FRACTION") { columnFractionsList.put(priorities.indexOf("high_vegetation"), "high_vegetation") columnFractionsList.put(priorities.indexOf("low_vegetation"), "low_vegetation") - } else if (it.equalsIgnoreCase("ROAD_FRACTION")) { + } else if (it=="ROAD_FRACTION") { columnFractionsList.put(priorities.indexOf("road"), "road") - } else if (it.equalsIgnoreCase("IMPERVIOUS_FRACTION")) { + } else if (it=="IMPERVIOUS_FRACTION") { columnFractionsList.put(priorities.indexOf("impervious"), "impervious") - } else if (it.equalsIgnoreCase("BUILDING_HEIGHT") && building) { + } else if (it=="BUILDING_HEIGHT" && building) { height_roof_unweighted_list.addAll( ["AVG", "STD"]) - } else if (it.equalsIgnoreCase("BUILDING_HEIGHT_WEIGHTED") && building) { + } else if (it=="BUILDING_HEIGHT_WEIGHTED" && building) { weightedBuildingIndicators["height_roof"] = ["area": ["AVG", "STD"]] - } else if (it.equalsIgnoreCase("BUILDING_POP") && building) { + } else if (it=="BUILDING_POP" && building) { unweightedBuildingIndicators.put("pop", ["SUM"]) } - else if (it.equalsIgnoreCase("HEIGHT_OF_ROUGHNESS_ELEMENTS") && building) { + else if (it=="HEIGHT_OF_ROUGHNESS_ELEMENTS" && building) { height_roof_unweighted_list.add("GEOM_AVG") } } @@ -1818,7 +1823,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if (list_indicators*.toUpperCase().contains("BUILDING_TYPE") && building) { + if (list_indicators_upper.contains("BUILDING_TYPE") && building) { if (!buildingCutted) { buildingCutted = cutBuilding(datasource, grid, building) if(!buildingCutted){ @@ -1839,10 +1844,9 @@ String rasterizeIndicators(JdbcDataSource datasource, } } - def freeFacadeDensityExact - if ((list_indicators*.toUpperCase().containsAll(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO"]) && building) || - (list_indicators*.toUpperCase().contains("BUILDING_SURFACE_DENSITY") && building)) { + + if ((list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO", "BUILDING_SURFACE_DENSITY"]) && building)) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, @@ -1853,20 +1857,20 @@ String rasterizeIndicators(JdbcDataSource datasource, return } } - freeFacadeDensityExact = Geoindicators.RsuIndicators.freeExternalFacadeDensityExact(datasource, createScalesRelationsGridBl, + def freeFacadeDensityExact = Geoindicators.RsuIndicators.freeExternalFacadeDensityExact(datasource, createScalesRelationsGridBl, grid, grid_column_identifier, prefixName) if (freeFacadeDensityExact) { - if (list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY")) { + if (list_indicators_upper.intersect(["FREE_EXTERNAL_FACADE_DENSITY", "ASPECT_RATIO"])) { indicatorTablesToJoin.put(freeFacadeDensityExact, grid_column_identifier) } - if (list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY")) { + if (list_indicators_upper.contains("FREE_EXTERNAL_FACADE_DENSITY")) { def buildingSurfDensity = Geoindicators.RsuIndicators.buildingSurfaceDensity( datasource, freeFacadeDensityExact, surfaceFractionsProcess, "FREE_EXTERNAL_FACADE_DENSITY", "building_fraction", grid_column_identifier, prefixName) if (buildingSurfDensity) { - if (list_indicators*.toUpperCase().contains("BUILDING_SURFACE_DENSITY")) { + if (list_indicators_upper.contains("BUILDING_SURFACE_DENSITY")) { indicatorTablesToJoin.put(buildingSurfDensity, grid_column_identifier) } } @@ -1877,7 +1881,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if (list_indicators*.toUpperCase().contains("BUILDING_HEIGHT_DIST") && building) { + if (list_indicators_upper.contains("BUILDING_HEIGHT_DIST") && building) { if (!buildingCutted) { buildingCutted = cutBuilding(datasource, grid, building) if(!buildingCutted){ @@ -1895,7 +1899,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } } - if (list_indicators*.toUpperCase().contains("FRONTAL_AREA_INDEX") && building) { + if (list_indicators_upper.contains("FRONTAL_AREA_INDEX") && building) { if (!datasource.getTable(building).isEmpty()) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings @@ -1918,11 +1922,11 @@ String rasterizeIndicators(JdbcDataSource datasource, } } - if (list_indicators*.toUpperCase().contains("SEA_LAND_FRACTION") && sea_land_mask) { + if (list_indicators_upper.contains("SEA_LAND_FRACTION") && sea_land_mask) { // If only one type of surface (land or sea) is in the zone, no need for computational fraction calculation def sea_land_type_rows = datasource.rows("""SELECT $seaLandTypeField, COUNT(*) AS NB_TYPES FROM $sea_land_mask - GROUP BY $seaLandTypeField""") + GROUP BY $seaLandTypeField""".toString()) if (sea_land_type_rows[seaLandTypeField].count("sea") == 0) { datasource """ DROP TABLE IF EXISTS $seaLandFractionTab; @@ -1960,7 +1964,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if (list_indicators*.toUpperCase().contains("SVF") && building) { + if (list_indicators_upper.contains("SVF") && building) { if (!createScalesRelationsGridBl) { // Create the relations between grid cells and buildings createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin(datasource, building, @@ -1982,7 +1986,7 @@ String rasterizeIndicators(JdbcDataSource datasource, tablesToDrop< Date: Tue, 3 Oct 2023 18:02:40 +0200 Subject: [PATCH 05/12] Add TERRAIN_ROUGHNESS_CLASS at grid scale --- .../bdtopo/AbstractBDTopoWorkflow.groovy | 3 +- .../geoindicators/RsuIndicators.groovy | 6 +- .../WorkflowGeoIndicators.groovy | 190 +++++++++--------- .../geoindicators/RsuIndicatorsTests.groovy | 2 +- .../geoclimate/osm/WorkflowOSM.groovy | 3 +- .../geoclimate/osm/WorflowOSMTest.groovy | 6 +- 6 files changed, 109 insertions(+), 101 deletions(-) diff --git a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy index c4be56987d..2393453556 100644 --- a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy +++ b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy @@ -543,7 +543,8 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { def allowed_grid_indicators = ["BUILDING_FRACTION", "BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION", "WATER_FRACTION", "VEGETATION_FRACTION", "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION", "UTRF_FLOOR_AREA_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY", "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", - "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO", "SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS"] + "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_FRACTION", "ASPECT_RATIO", + "SVF", "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS"] def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase()) if (allowedOutputIndicators) { //Update the RSU indicators list according the grid indicators diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index ba151cb151..62919539b4 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -1094,6 +1094,7 @@ String linearRoadOperations(JdbcDataSource datasource, String rsuTable, String r * @param datasource A connexion to a database (H2GIS, PostGIS, ...) where are stored the input Table and in which * the resulting database will be stored * @param rsuTable the name of the input ITable where are stored the effectiveTerrainRoughnessHeight values + * @param id_rsu Unique identifier column name * @param effectiveTerrainRoughnessLength the field name corresponding to the RSU effective terrain roughness class due * to roughness elements (buildings, trees, etc.) (in the rsuTable) * @param prefixName String use as prefix to name the output table @@ -1102,8 +1103,7 @@ String linearRoadOperations(JdbcDataSource datasource, String rsuTable, String r * * @author Jérémy Bernard */ -String effectiveTerrainRoughnessClass(JdbcDataSource datasource, String rsu, String effectiveTerrainRoughnessLength, String prefixName) { - def ID_COLUMN_RSU = "id_rsu" +String effectiveTerrainRoughnessClass(JdbcDataSource datasource, String rsu, String id_rsu, String effectiveTerrainRoughnessLength, String prefixName) { def BASE_NAME = "effective_terrain_roughness_class" debug "Executing RSU effective terrain roughness class" @@ -1113,7 +1113,7 @@ String effectiveTerrainRoughnessClass(JdbcDataSource datasource, String rsu, Str // Based on the lookup Table of Davenport datasource """DROP TABLE IF EXISTS $outputTableName; - CREATE TABLE $outputTableName AS SELECT $ID_COLUMN_RSU, + CREATE TABLE $outputTableName AS SELECT $id_rsu, CASEWHEN($effectiveTerrainRoughnessLength<0.0 OR $effectiveTerrainRoughnessLength IS NULL, null, CASEWHEN($effectiveTerrainRoughnessLength<0.00035, 1, CASEWHEN($effectiveTerrainRoughnessLength<0.01525, 2, diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index 79e947845a..c483927913 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -367,7 +367,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, angleRangeSizeBuDirection : 30, svfSimplified : true, indicatorUse : ["LCZ", "UTRF", "TEB"], - surfSuperpositions : ["high_vegetation": ["water", "building", "low_vegetation", "rail","road", "impervious"]], + surfSuperpositions : ["high_vegetation": ["water", "building", "low_vegetation", "rail", "road", "impervious"]], surfPriorities : ["water", "building", "high_vegetation", "low_vegetation", "rail", "road", "impervious"], buildingAreaTypeAndComposition: ["light_industry": ["light_industry"], "heavy_industry": ["heavy_industry"], @@ -422,7 +422,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, def columnIdRsu = "id_rsu" def columnIdBuild = "id_build" def BASE_NAME = "rsu_indicators" - def tablesToDrop=[] + def tablesToDrop = [] // Maps for intermediate or final joins def finalTablesToJoin = [:] def intermediateJoin = [:] @@ -450,9 +450,9 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, if (!surfaceFractions) { // Calculate all surface fractions indicators // Need to create the smallest geometries used as input of the surface fraction process - def superpositionsTable = Geoindicators.RsuIndicators.smallestCommunGeometry(datasource, + def superpositionsTable = Geoindicators.RsuIndicators.smallestCommunGeometry(datasource, rsu, "id_rsu", buildingTable, road, water, vegetation, - impervious,rail, temporaryPrefName) + impervious, rail, temporaryPrefName) if (!superpositionsTable) { info "Cannot compute the smallest commun geometries" return @@ -467,7 +467,7 @@ String computeRSUIndicators(JdbcDataSource datasource, String buildingTable, return } - tablesToDrop< urbTypoCorrespondenceTabInverted[ini] = fin } - datasource.createIndex(utrfBuild,"I_TYPO") + datasource.createIndex(utrfBuild, "I_TYPO") def queryDistinct = """SELECT DISTINCT I_TYPO AS I_TYPO FROM $utrfBuild""" def mapTypos = datasource.rows(queryDistinct.toString()) def listTypos = [] @@ -912,8 +912,8 @@ Map computeTypologyIndicators(JdbcDataSource datasource, String building_indicat } queryCaseWhenReplace = queryCaseWhenReplace + " 'unknown' " + endCaseWhen utrfBuilding = postfix "UTRF_BUILDING" - datasource.createIndex(utrfBuild,COLUMN_ID_BUILD) - datasource.createIndex(building_indicators,COLUMN_ID_BUILD) + datasource.createIndex(utrfBuild, COLUMN_ID_BUILD) + datasource.createIndex(building_indicators, COLUMN_ID_BUILD) datasource """ DROP TABLE IF EXISTS $utrfBuilding; CREATE TABLE $utrfBuilding AS SELECT a.$COLUMN_ID_BUILD, a.$COLUMN_ID_RSU, a.THE_GEOM, @@ -926,9 +926,9 @@ Map computeTypologyIndicators(JdbcDataSource datasource, String building_indicat def queryCasewhen = [:] queryCasewhen["AREA"] = "" queryCasewhen["FLOOR_AREA"] = "" - datasource.createIndex(rsu_indicators,COLUMN_ID_RSU) - datasource.createIndex(building_indicators,COLUMN_ID_RSU) - datasource.createIndex(utrfBuilding,COLUMN_ID_BUILD) + datasource.createIndex(rsu_indicators, COLUMN_ID_RSU) + datasource.createIndex(building_indicators, COLUMN_ID_RSU) + datasource.createIndex(utrfBuilding, COLUMN_ID_BUILD) queryCasewhen.keySet().each { ind -> def querySum = "" @@ -946,26 +946,26 @@ Map computeTypologyIndicators(JdbcDataSource datasource, String building_indicat WHERE b.$COLUMN_ID_RSU IS NOT NULL GROUP BY b.$COLUMN_ID_RSU """.toString() - tablesToDrop< Date: Tue, 3 Oct 2023 20:33:25 +0200 Subject: [PATCH 06/12] Fix #831 --- .../org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy index 46107ee6da..8280267cb0 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -124,7 +124,7 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz datasource """ DROP TABLE IF EXISTS $outputTableName; CREATE TABLE $outputTableName AS - SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(a.the_geom, $epsg) AS the_geom + SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(st_buffer(a.the_geom, -0.01), $epsg) AS the_geom FROM ST_EXPLODE('( SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom FROM $inputTableName)') AS a, @@ -136,7 +136,7 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz datasource """ DROP TABLE IF EXISTS $outputTableName; CREATE TABLE $outputTableName AS - SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(ST_FORCE2D(the_geom), $epsg) AS the_geom + SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(st_buffer(a.the_geom, -0.01), $epsg) AS the_geom FROM ST_EXPLODE('( SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom FROM $inputTableName)') where st_area(the_geom) > $area""".toString() From c8e51b4ce5903fb3404a418fb288043c480e1cf7 Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 3 Oct 2023 21:03:12 +0200 Subject: [PATCH 07/12] Fix #831 --- .../geoclimate/geoindicators/GenericIndicators.groovy | 4 ++-- .../orbisgis/geoclimate/geoindicators/SpatialUnits.groovy | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy index fef695a134..21b9025751 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy @@ -447,7 +447,7 @@ String distributionCharacterization(JdbcDataSource datasource, String distribTab def outputTableName = prefix prefixName, BASENAME // Get all columns from the distribution table and remove the geometry column if exists - def allColumns = datasource."$distribTableName".columns + def allColumns = datasource.getColumnNames(distribTableName) if (allColumns.contains(GEOMETRY_FIELD)) { allColumns -= GEOMETRY_FIELD } @@ -455,7 +455,7 @@ String distributionCharacterization(JdbcDataSource datasource, String distribTab def distribColumns = allColumns.minus(inputId.toUpperCase()) def nbDistCol = distribColumns.size() - if(distribColumns.size==0){ + if(distribColumns.size()==0){ error("Any columns to compute the distribution from the table $distribTableName".toString()) return } diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy index 8280267cb0..48d4fa77fa 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -124,9 +124,9 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz datasource """ DROP TABLE IF EXISTS $outputTableName; CREATE TABLE $outputTableName AS - SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(st_buffer(a.the_geom, -0.01), $epsg) AS the_geom + SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(a.the_geom, $epsg) AS the_geom FROM ST_EXPLODE('( - SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom + SELECT ST_BUFFER(ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))), -0.01) AS the_geom FROM $inputTableName)') AS a, $inputzone AS b WHERE a.the_geom && b.the_geom @@ -136,9 +136,9 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz datasource """ DROP TABLE IF EXISTS $outputTableName; CREATE TABLE $outputTableName AS - SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(st_buffer(a.the_geom, -0.01), $epsg) AS the_geom + SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(st_buffer(the_geom, -0.01), $epsg) AS the_geom FROM ST_EXPLODE('( - SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom + SELECT ST_BUFFER(ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))),-0.01) AS the_geom FROM $inputTableName)') where st_area(the_geom) > $area""".toString() } From 0efec82a9209d7fbe5bf21b5da9dbebb73ec6b75 Mon Sep 17 00:00:00 2001 From: ebocher Date: Wed, 4 Oct 2023 13:38:30 +0200 Subject: [PATCH 08/12] Fix free external density indicators when any building share a facade Add rasterizeIndicators test --- .../geoindicators/RsuIndicators.groovy | 15 +++++- .../WorkflowGeoIndicators.groovy | 12 +++-- .../geoindicators/RsuIndicatorsTests.groovy | 24 +++++++++ .../WorkflowGeoIndicatorsTest.groovy | 50 +++++++++++++++++++ 4 files changed, 95 insertions(+), 6 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index 62919539b4..81701554a1 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -170,13 +170,24 @@ String freeExternalFacadeDensityExact(JdbcDataSource datasource, String building // 4. Calculates the free facade density by RSU (subtract 3 and 2 and divide by RSU area) datasource.createIndex(buildLineRsu,idRsu) datasource.createIndex(sharedLineRsu,idRsu) - datasource """ + if(datasource.getRowCount(sharedLineRsu)>0) { + datasource """ DROP TABLE IF EXISTS $onlyBuildRsu; CREATE TABLE $onlyBuildRsu AS SELECT a.$idRsu, (a.$FACADE_AREA-b.$FACADE_AREA)/a.$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY FROM $buildLineRsu AS a LEFT JOIN $sharedLineRsu AS b - ON a.$idRsu = b.$idRsu""".toString() + ON a.$idRsu = b.$idRsu + """.toString() + }else{ + datasource """ + DROP TABLE IF EXISTS $onlyBuildRsu; + CREATE TABLE $onlyBuildRsu + AS SELECT $idRsu, + $FACADE_AREA/$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY + FROM $buildLineRsu + """.toString() + } // 5. Join RSU having no buildings and set their value to 0 datasource.createIndex(onlyBuildRsu,idRsu) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index c483927913..298d4f12f5 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -1613,7 +1613,7 @@ Map computeGeoclimateIndicators(JdbcDataSource datasource, String zone, String b * @return */ String rasterizeIndicators(JdbcDataSource datasource, - String grid, List list_indicators = [], + String grid, List list_indicators, String building, String road, String vegetation, String water, String impervious, String rsu_lcz, String rsu_utrf_area, String rsu_utrf_floor_area, String sea_land_mask, @@ -1822,7 +1822,7 @@ String rasterizeIndicators(JdbcDataSource datasource, } - if (list_indicators_upper.contains("BUILDING_TYPE") && building) { + if (list_indicators_upper.contains("BUILDING_TYPE_FRACTION") && building) { if (!buildingCutted) { buildingCutted = cutBuilding(datasource, grid, building) if (!buildingCutted) { @@ -1838,6 +1838,7 @@ String rasterizeIndicators(JdbcDataSource datasource, "building_type_fraction") if (upperScaleAreaStatistics) { indicatorTablesToJoin.put(upperScaleAreaStatistics, grid_column_identifier) + tablesToDrop< Date: Wed, 4 Oct 2023 13:57:40 +0200 Subject: [PATCH 09/12] Fix test --- .../geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy index d8474d3b64..c6d5b2877e 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy @@ -453,7 +453,6 @@ class WorkflowGeoIndicatorsTest { assertNotNull(grid_indicators) assertEquals(1, datasource.getRowCount(grid_indicators)) def rows = datasource.firstRow("select * from $grid_indicators".toString()) - println(rows) assertEquals(1d, rows.BUILDING_FRACTION) assertEquals(0d, rows.HIGH_VEGETATION_FRACTION) assertEquals(0d, rows.IMPERVIOUS_FRACTION) @@ -473,7 +472,7 @@ class WorkflowGeoIndicatorsTest { assertEquals(8, rows.EFFECTIVE_TERRAIN_ROUGHNESS_CLASS) assertEquals(2d, rows.ASPECT_RATIO) assertTrue(0.5 -rows.SVF<0.1) - assertEquals(1, rows.TYPE_OFFICE) + assertEquals(1d, rows.TYPE_OFFICE) } From b3242599e146dcfef1282474ffba4ccdbf1dcf72 Mon Sep 17 00:00:00 2001 From: ebocher Date: Wed, 4 Oct 2023 14:13:13 +0200 Subject: [PATCH 10/12] Fix aspect ratio must be null when building_fraction = 1 --- .../geoclimate/geoindicators/WorkflowGeoIndicators.groovy | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy index 298d4f12f5..db1146f88e 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -2062,7 +2062,7 @@ String rasterizeIndicators(JdbcDataSource datasource, //Because all other indicators have been yet computed we just alter the table with the aspect_ratio column datasource.execute(""" ALTER TABLE $grid_indicators_table ADD COLUMN ASPECT_RATIO DOUBLE PRECISION; - UPDATE $grid_indicators_table SET ASPECT_RATIO = case when building_fraction = 1 then free_external_facade_density else free_external_facade_density / (1 - building_fraction) end + UPDATE $grid_indicators_table SET ASPECT_RATIO = case when building_fraction = 1 then null else free_external_facade_density / (1 - building_fraction) end """.toString()) } tablesToDrop << createScalesRelationsGridBl From 318f8cc302f0f04cf172563d9dbf3470de3a71e2 Mon Sep 17 00:00:00 2001 From: ebocher Date: Wed, 4 Oct 2023 15:10:50 +0200 Subject: [PATCH 11/12] Last fix on join --- .../geoindicators/RsuIndicators.groovy | 20 ++++++++----------- .../geoindicators/DataUtilsTests.groovy | 3 ++- .../geoindicators/RsuIndicatorsTests.groovy | 7 +++++-- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy index 81701554a1..e6d8a1497d 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -170,24 +170,20 @@ String freeExternalFacadeDensityExact(JdbcDataSource datasource, String building // 4. Calculates the free facade density by RSU (subtract 3 and 2 and divide by RSU area) datasource.createIndex(buildLineRsu,idRsu) datasource.createIndex(sharedLineRsu,idRsu) - if(datasource.getRowCount(sharedLineRsu)>0) { - datasource """ + datasource """ DROP TABLE IF EXISTS $onlyBuildRsu; CREATE TABLE $onlyBuildRsu AS SELECT a.$idRsu, - (a.$FACADE_AREA-b.$FACADE_AREA)/a.$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY + a.$FACADE_AREA/a.$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY FROM $buildLineRsu AS a LEFT JOIN $sharedLineRsu AS b + ON a.$idRsu = b.$idRsu WHERE b.$idRsu IS NULL + union all + SELECT a.$idRsu, + (a.$FACADE_AREA-b.$FACADE_AREA)/a.$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY + FROM $buildLineRsu AS a right JOIN $sharedLineRsu AS b ON a.$idRsu = b.$idRsu """.toString() - }else{ - datasource """ - DROP TABLE IF EXISTS $onlyBuildRsu; - CREATE TABLE $onlyBuildRsu - AS SELECT $idRsu, - $FACADE_AREA/$RSU_AREA AS FREE_EXTERNAL_FACADE_DENSITY - FROM $buildLineRsu - """.toString() - } + // 5. Join RSU having no buildings and set their value to 0 datasource.createIndex(onlyBuildRsu,idRsu) diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/DataUtilsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/DataUtilsTests.groovy index 8934d3c31f..c637bd8ffe 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/DataUtilsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/DataUtilsTests.groovy @@ -23,6 +23,7 @@ import org.junit.jupiter.api.BeforeAll import org.junit.jupiter.api.BeforeEach import org.junit.jupiter.api.Test import org.junit.jupiter.api.io.TempDir +import org.orbisgis.data.H2GIS import org.orbisgis.geoclimate.Geoindicators import static org.orbisgis.data.H2GIS.open @@ -31,7 +32,7 @@ class DataUtilsTests { @TempDir static File folder - private static def h2GIS + private static H2GIS h2GIS @BeforeAll static void beforeAll() { diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy index 3edbe7cb6f..8690354c73 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy @@ -112,9 +112,11 @@ class RsuIndicatorsTests { CREATE TABLE tempo_build(id_build int, the_geom geometry, height_wall double); INSERT INTO tempo_build VALUES (1, 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'::GEOMETRY, 10), (2, 'POLYGON ((10 0, 20 0, 20 20, 10 20, 10 0))'::GEOMETRY, 10), - (3, 'POLYGON ((30 30, 50 30, 50 50, 30 50, 30 30))'::GEOMETRY, 10); + (3, 'POLYGON ((30 30, 50 30, 50 50, 30 50, 30 30))'::GEOMETRY, 10), + (4, 'POLYGON ((120 60, 130 60, 130 50, 120 50, 120 60))'::GEOMETRY, 10); CREATE TABLE tempo_rsu(id_rsu int, the_geom geometry); - INSERT INTO tempo_rsu VALUES (1, 'POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))'::GEOMETRY); + INSERT INTO tempo_rsu VALUES (1, 'POLYGON((0 0, 100 0, 100 100, 0 100, 0 0))'::GEOMETRY), + (2, 'POLYGON((100 100, 200 100, 200 0, 100 0, 100 100))'::GEOMETRY) ; """ // First calculate the correlation table between buildings and rsu def buildingTableRelation = Geoindicators.SpatialUnits.spatialJoin(h2GIS, @@ -126,6 +128,7 @@ class RsuIndicatorsTests { "id_rsu", "test") assertNotNull(p) assertEquals 0.16, h2GIS.firstRow("SELECT * FROM ${p} WHERE id_rsu = 1").FREE_EXTERNAL_FACADE_DENSITY + assertEquals(0.04, h2GIS.firstRow("SELECT * FROM ${p} WHERE id_rsu = 2").FREE_EXTERNAL_FACADE_DENSITY) } @Test From 6fd9b902bd92674239840be53554b4bfbd96bba6 Mon Sep 17 00:00:00 2001 From: ebocher Date: Wed, 4 Oct 2023 15:25:45 +0200 Subject: [PATCH 12/12] Last fix on join --- .../geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy index c6d5b2877e..6d607008cb 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy @@ -470,7 +470,7 @@ class WorkflowGeoIndicatorsTest { assertEquals(3d, rows.BUILDING_SURFACE_DENSITY) assertTrue(1.5 -rows.EFFECTIVE_TERRAIN_ROUGHNESS_LENGTH<0.001) assertEquals(8, rows.EFFECTIVE_TERRAIN_ROUGHNESS_CLASS) - assertEquals(2d, rows.ASPECT_RATIO) + assertNull(rows.ASPECT_RATIO) assertTrue(0.5 -rows.SVF<0.1) assertEquals(1d, rows.TYPE_OFFICE)