diff --git a/bdtopo_v2/src/main/groovy/org/orbisgis/geoclimate/bdtopo_v2/WorkflowBDTopo_V2.groovy b/bdtopo_v2/src/main/groovy/org/orbisgis/geoclimate/bdtopo_v2/WorkflowBDTopo_V2.groovy index 91f3557201..e61ec08738 100644 --- a/bdtopo_v2/src/main/groovy/org/orbisgis/geoclimate/bdtopo_v2/WorkflowBDTopo_V2.groovy +++ b/bdtopo_v2/src/main/groovy/org/orbisgis/geoclimate/bdtopo_v2/WorkflowBDTopo_V2.groovy @@ -1250,7 +1250,7 @@ def extractProcessingParameters(def processing_parameters){ } def allowed_grid_indicators=["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP","BUILDING_TYPE_FRACTION","WATER_FRACTION","VEGETATION_FRACTION", "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY","FREE_EXTERNAL_FACADE_DENSITY", - "BUILDING_HEIGHT_WEIGHTED"] + "BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", "BUILDING_HEIGHT_DIST", "FRONTAL_AREA_INDEX", "SEA_LAND_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/GenericIndicators.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy index e23604e483..e3fb139d6d 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy @@ -1073,7 +1073,7 @@ IProcess upperScaleAreaStatistics() { """ }else{ query += """ - , CASE WHEN $lowerColumnName=${it.val} + , CASE WHEN $lowerColumnName='${it.val}' THEN SUM(area) ELSE 0 END AS $aliasColumn """ 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 bac53d4d13..1f33699cf5 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicators.groovy @@ -77,11 +77,15 @@ IProcess freeExternalFacadeDensity() { * in the freeExternalFacadeDensity IProcess since in this process only the part of building being in the RSU * is considered in the calculation of the indicator. * + * WARNING: WITH THE CURRENT METHOD, IF A BUILDING FACADE FOLLOW THE LINE SEPARATING TWO RSU, IT WILL BE COUNTED FOR BOTH RSU + * (EVEN THOUGH THIS SITUATION IS PROBABLY ALMOST IMPOSSIBLE WITH REGULAR SQUARE GRID)... + * * @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 buildingTable The name of the input ITable where are stored the buildings geometry, the building height, * the building and the rsu id - * @param rsuTable The name of the input ITable where are stored the rsu geometries and the id_rsu + * @param rsuTable The name of the input ITable where are stored the rsu geometries and the 'idRsu' + * @param idRsu the name of the id of the RSU table * @param prefixName String use as prefix to name the output table * * @return A database table name. @@ -101,7 +105,7 @@ IProcess freeExternalFacadeDensityExact() { def HEIGHT_WALL = "height_wall" def FACADE_AREA = "facade_area" def RSU_AREA = "rsu_area" - def BASE_NAME = "_exact_free_external_facade_density" + def BASE_NAME = "exact_free_external_facade_density" debug "Executing RSU free external facade density (exact version)" @@ -1686,3 +1690,430 @@ IProcess surfaceFractions() { } } } + +/** + * Process used to compute the sum of all building facades (free facades and roofs) included in a + * Reference Spatial Unit (RSU) divided by the RSU area. The calculation is actually simply the sum + * of two other indicators (building fraction and free external facade density)... + * + * @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 facadeDensityTable The name of the input ITable where are stored the facade density and the rsu id + * @param buildingFractionTable The name of the input ITable where are stored the building fraction and the rsu id + * @param facDensityColumn The name of the column where are stored the facade density values (within the + * 'facadeDensityTable' Table) + * @param buFractionColumn The name of the column where are stored the building fraction values (within the + * 'buildingFractionTable' Table) + * @param idRsu the name of the id of the RSU table + * @param prefixName String use as prefix to name the output table + * + * @return A database table name. + * @author Jérémy Bernard + */ +IProcess buildingSurfaceDensity() { + return create { + title "RSU building surface density" + id "buildingSurfaceDensity" + inputs facadeDensityTable: String, buildingFractionTable: String, facDensityColumn: String, + buFractionColumn: String, idRsu: String, prefixName: String, datasource: JdbcDataSource + outputs outputTableName: String + run { facadeDensityTable, buildingFractionTable, facDensityColumn, buFractionColumn, idRsu, prefixName, datasource -> + + def BASE_NAME = "building_surface_fraction" + + debug "Executing building surface density" + + // The name of the outputTableName is constructed + def outputTableName = prefix prefixName, BASE_NAME + + // Sum free facade density and building fraction... + datasource."$facadeDensityTable"."$idRsu".createIndex() + datasource."$buildingFractionTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS $outputTableName; + CREATE TABLE $outputTableName + AS SELECT a.$idRsu, + a.$buFractionColumn + b.$facDensityColumn AS BUILDING_SURFACE_DENSITY + FROM $buildingFractionTable AS a LEFT JOIN $facadeDensityTable AS b + ON a.$idRsu = b.$idRsu""".toString() + + [outputTableName: outputTableName] + } + } +} + +/** + * Script to compute the distributio 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. + * + * @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 buildingTable the name of the input ITable where are stored the buildings and the relationships + * between buildings and RSU + * @param rsuTable the name of the input ITable where are stored the RSU + * @param idRsu the name of the id of the RSU table + * @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]) + * + * @return outputTableName Table name in which the rsu id and their corresponding indicator value are stored + * + * @author Jérémy Bernard + */ +IProcess roofFractionDistributionExact() { + return create { + title "RSU exact horizontal roof area fraction distribution" + id "roofFractionDistributionExact" + inputs rsuTable: String, buildingTable: String, idRsu: String, listLayersBottom: [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], + prefixName: String, datasource: JdbcDataSource, density: true + outputs outputTableName: String + run { rsuTable, buildingTable, idRsu, listLayersBottom, prefixName, datasource, density -> + + def GEOMETRIC_COLUMN_RSU = "the_geom" + def GEOMETRIC_COLUMN_BU = "the_geom" + def ID_COLUMN_BU = "id_build" + def HEIGHT_WALL = "height_wall" + def HEIGHT_ROOF = "height_roof" + def BUILDING_AREA = "building_area" + def BUILD_HEIGHT = "BUILD_HEIGHT" + def BASE_NAME = "roof_fraction_distribution" + + // The name of the outputTableName is constructed + def outputTableName = prefix prefixName, BASE_NAME + + debug "Executing RSU roof area distribution (and optionally roof density)" + + // To avoid overwriting the output files of this step, a unique identifier is created + // Temporary table names + def buildInter = postfix "build_inter" + 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."$buildingTable"."$ID_COLUMN_BU".createIndex() + datasource."$rsuTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS $buildInter; + CREATE TABLE $buildInter + AS SELECT a.$ID_COLUMN_BU, a.$idRsu, + ST_INTERSECTION(a.$GEOMETRIC_COLUMN_BU, b.$GEOMETRIC_COLUMN_RSU) AS $GEOMETRIC_COLUMN_BU, + (a.$HEIGHT_WALL + a.$HEIGHT_ROOF) / 2 AS $BUILD_HEIGHT + FROM $buildingTable AS a LEFT JOIN $rsuTable AS b + ON a.$idRsu = b.$idRsu""".toString() + + // 2. Calculate the total building roof area within each RSU + datasource."$buildInter"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS $rsuBuildingArea; + CREATE TABLE $rsuBuildingArea + AS SELECT $idRsu, SUM(ST_AREA($GEOMETRIC_COLUMN_BU)) AS $BUILDING_AREA + FROM $buildInter + 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() + def tab_H = [:] + def indicToJoin = [:] + for (i in 1..(listLayersBottom.size() - 1)) { + def layer_top = listLayersBottom[i] + def layer_bottom = listLayersBottom[i-1] + def indicNameH = "${BASE_NAME}_${layer_bottom}_$layer_top".toString() + tab_H[i-1] = "${buildFracH}_$layer_bottom".toString() + datasource """ + DROP TABLE IF EXISTS $bufferTable; + CREATE TABLE $bufferTable + AS SELECT a.$idRsu, + CASE WHEN a.$BUILDING_AREA = 0 + THEN 0 + ELSE SUM(ST_AREA(b.$GEOMETRIC_COLUMN_BU)) / a.$BUILDING_AREA + END AS $indicNameH + FROM $rsuBuildingArea AS a LEFT JOIN $buildInter AS b + ON a.$idRsu = b.$idRsu + WHERE b.$BUILD_HEIGHT >= $layer_bottom AND b.$BUILD_HEIGHT < $layer_top + GROUP BY b.$idRsu""".toString() + // Fill missing values with 0 + datasource."$bufferTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS ${tab_H[i-1]}; + CREATE TABLE ${tab_H[i-1]} + AS SELECT a.$idRsu, + COALESCE(b.$indicNameH, 0) AS $indicNameH + FROM $rsuTable AS a LEFT JOIN $bufferTable AS b + ON a.$idRsu = b.$idRsu""".toString() + // Declare this layer to the layer to join at the end + indicToJoin.put(tab_H[i-1], idRsu) + } + + // 4. Calculate the fraction of roof for the last level of the canopy + def layer_bottom = listLayersBottom[listLayersBottom.size() - 1] + def indicNameH = "${BASE_NAME}_${layer_bottom}_inf".toString() + tab_H[listLayersBottom.size() - 1] = "${buildFracH}_$layer_bottom".toString() + datasource """ + DROP TABLE IF EXISTS $bufferTable; + CREATE TABLE $bufferTable + AS SELECT a.$idRsu, + CASE WHEN a.$BUILDING_AREA = 0 + THEN 0 + ELSE SUM(ST_AREA(b.$GEOMETRIC_COLUMN_BU)) / a.$BUILDING_AREA + END AS $indicNameH + FROM $rsuBuildingArea AS a LEFT JOIN $buildInter AS b + ON a.$idRsu = b.$idRsu + WHERE b.$BUILD_HEIGHT >= $layer_bottom + GROUP BY b.$idRsu""".toString() + // Fill missing values with 0 + datasource."$bufferTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS ${tab_H[listLayersBottom.size() - 1]}; + CREATE TABLE ${tab_H[listLayersBottom.size() - 1]} + AS SELECT a.$idRsu, + COALESCE(b.$indicNameH, 0) AS $indicNameH + FROM $rsuTable AS a LEFT JOIN $bufferTable AS b + ON a.$idRsu = b.$idRsu""".toString() + // Declare this layer to the layer to join at the end + indicToJoin.put(tab_H[listLayersBottom.size() - 1], idRsu) + + // 5. Join all layers in one table + def joinGrids = Geoindicators.DataUtils.joinTables() + if (!joinGrids([inputTableNamesWithId: indicToJoin, + outputTableName : outputTableName, + datasource : datasource])) { + info "Cannot merge all indicators in RSU table $outputTableName." + return + } + + datasource """DROP TABLE IF EXISTS $buildInter, $rsuBuildingArea, $bufferTable, + ${tab_H.values().join(",")}""".toString() + + [outputTableName: outputTableName] + } + } +} + +/** + * Process used to compute the frontal area index for different combinations of direction / canopy levels. The frontal + * area index for a given set of direction / level is defined as the projected area of facade facing the direction + * of analysis divided by the horizontal surface of analysis and divided by the height of the layer + * + * WARNING: WITH THE CURRENT METHOD, IF A BUILDING FACADE FOLLOW THE LINE SEPARATING TWO RSU, IT WILL BE COUNTED FOR BOTH RSU + * (EVEN THOUGH THIS SITUATION IS PROBABLY ALMOST IMPOSSIBLE WITH REGULAR SQUARE GRID)... + * + * @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 buildingTable The name of the input ITable where are stored the buildings geometry, the building height, + * the building and the rsu id + * @param rsuTable The name of the input ITable where are stored the rsu geometries and the 'idRsu' + * @param idRsu the name of the id of the RSU 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 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 prefixName String use as prefix to name the output table + * + * @return A database table name. + * @author Jérémy Bernard + */ +IProcess frontalAreaIndexDistribution() { + return create { + title "RSU frontal area index distribution" + id "frontalAreaIndexDistribution" + inputs buildingTable: String, rsuTable: String, idRsu: String, listLayersBottom: [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], + numberOfDirection: 12, prefixName: String, datasource: JdbcDataSource + outputs outputTableName: String + run { buildingTable, rsuTable, idRsu, listLayersBottom, numberOfDirection, prefixName, datasource -> + + def GEOMETRIC_FIELD_RSU = "the_geom" + def GEOMETRIC_FIELD_BU = "the_geom" + def ID_FIELD_BU = "id_build" + def HEIGHT_WALL = "height_wall" + def BASE_NAME = "frontal_area_index_distribution" + + debug "Executing RSU frontal area index distribution" + + // The name of the outputTableName is constructed + def outputTableName = prefix prefixName, BASE_NAME + + if (360 % numberOfDirection == 0 && numberOfDirection % 2 == 0) { + + // Temporary table names + def buildLine = postfix "buildLine" + def allLinesRsu = postfix "all_lines_rsu" + def buildFracH = postfix "build_frac_H" + def bufferTable = postfix "buffer_table" + + // Consider facades as touching each other within a snap tolerance + def snap_tolerance = 0.01 + + // 1. Convert the building polygons into lines and create the intersection with RSU polygons + datasource."$buildingTable"."$idRsu".createIndex() + datasource."$rsuTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS $buildLine; + CREATE TABLE $buildLine + AS SELECT a.$ID_FIELD_BU, a.$idRsu, + ST_INTERSECTION(ST_TOMULTILINE(a.$GEOMETRIC_FIELD_BU), b.$GEOMETRIC_FIELD_RSU) AS $GEOMETRIC_FIELD_BU, + a.$HEIGHT_WALL + FROM $buildingTable AS a LEFT JOIN $rsuTable AS b + ON a.$idRsu = b.$idRsu""".toString() + + // 2. Keep only intersected facades within a given distance and calculate their length, height and azimuth + datasource."$buildLine"."$GEOMETRIC_FIELD_BU".createSpatialIndex() + datasource."$buildLine"."$idRsu".createIndex() + datasource."$buildLine"."$ID_FIELD_BU".createIndex() + datasource """ + DROP TABLE IF EXISTS $allLinesRsu; + CREATE TABLE $allLinesRsu + AS SELECT -ST_LENGTH($GEOMETRIC_FIELD_BU) AS LENGTH, + ST_AZIMUTH(ST_STARTPOINT($GEOMETRIC_FIELD_BU), ST_ENDPOINT($GEOMETRIC_FIELD_BU)) AS AZIMUTH, + $idRsu, + $HEIGHT_WALL, + $ID_FIELD_BU + FROM ST_EXPLODE('(SELECT ST_TOMULTISEGMENTS(ST_INTERSECTION(a.$GEOMETRIC_FIELD_BU, + ST_SNAP(b.$GEOMETRIC_FIELD_BU, + a.$GEOMETRIC_FIELD_BU, + $snap_tolerance))) AS $GEOMETRIC_FIELD_BU, + LEAST(a.$HEIGHT_WALL, b.$HEIGHT_WALL) AS $HEIGHT_WALL, + a.$idRsu, a.$ID_FIELD_BU + FROM $buildLine AS a LEFT JOIN $buildLine AS b + ON a.$idRsu = b.$idRsu + WHERE a.$GEOMETRIC_FIELD_BU && b.$GEOMETRIC_FIELD_BU AND ST_INTERSECTS(a.$GEOMETRIC_FIELD_BU, + ST_SNAP(b.$GEOMETRIC_FIELD_BU, a.$GEOMETRIC_FIELD_BU, $snap_tolerance)) AND + a.$ID_FIELD_BU <> b.$ID_FIELD_BU)') + WHERE ST_DIMENSION($GEOMETRIC_FIELD_BU) = 1 + UNION ALL + SELECT ST_LENGTH($GEOMETRIC_FIELD_BU) AS LENGTH, + ST_AZIMUTH(ST_STARTPOINT($GEOMETRIC_FIELD_BU), ST_ENDPOINT($GEOMETRIC_FIELD_BU)) AS AZIMUTH, + $idRsu, + $HEIGHT_WALL, + $ID_FIELD_BU + FROM ST_EXPLODE('(SELECT ST_TOMULTISEGMENTS($GEOMETRIC_FIELD_BU) AS $GEOMETRIC_FIELD_BU, + $HEIGHT_WALL, + $idRsu, $ID_FIELD_BU + FROM $buildLine)') + WHERE ST_LENGTH($GEOMETRIC_FIELD_BU) > 0;""".toString() + + // 3. Make the calculations for all directions of each level except the highest one + def dirQueryVertFrac = [:] + def dirQueryDiv = [:] + def angleRangeRad = 2 * Math.PI / numberOfDirection + def angleRangeDeg = 360 / numberOfDirection + def tab_H = [:] + def indicToJoin = [:] + datasource."$rsuTable"."$idRsu".createIndex() + for (i in 1..(listLayersBottom.size() - 1)) { + def layer_top = listLayersBottom[i] + def layer_bottom = listLayersBottom[i-1] + def deltaH = layer_top - layer_bottom + tab_H[i-1] = "${buildFracH}_$layer_bottom".toString() + + // Define queries and indic names + def dirList = [:] + (0..numberOfDirection-1).each{ dirList[it] = (it+0.5)*angleRangeRad} + dirList.each{k, v -> + // 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] = """ + 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)/$deltaH + 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)/$deltaH + END + ELSE 0 + END + END AS $indicName""" + dirQueryDiv[k] = """COALESCE(SUM(b.$indicName)/ST_AREA(a.$GEOMETRIC_FIELD_RSU), 0) AS $indicName""" + } + // Calculates projected surfaces for buildings and shared facades + datasource """ + DROP TABLE IF EXISTS $bufferTable; + CREATE TABLE $bufferTable + AS SELECT $idRsu, + ${dirQueryVertFrac.values().join(",")} + FROM $allLinesRsu + WHERE $HEIGHT_WALL > $layer_bottom""".toString() + // Fill missing values with 0 + datasource."$bufferTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS ${tab_H[i-1]}; + CREATE TABLE ${tab_H[i-1]} + AS SELECT a.$idRsu, + ${dirQueryDiv.values().join(",")} + FROM $rsuTable AS a LEFT JOIN $bufferTable AS b + ON a.$idRsu = b.$idRsu + GROUP BY a.$idRsu""".toString() + // Declare this layer to the layer to join at the end + indicToJoin.put(tab_H[i-1], idRsu) + } + + // 4. Make the calculations for the last level + def layer_bottom = listLayersBottom[listLayersBottom.size() - 1] + // Get the maximum building height + def layer_top = (datasource.firstRow("SELECT MAX($HEIGHT_WALL) AS MAXH FROM $buildingTable").MAXH).trunc()+1 as int + def deltaH = layer_top - layer_bottom + tab_H[listLayersBottom.size() - 1] = "${buildFracH}_$layer_bottom".toString() + + // Define queries and indic names + def dirList = [:] + (0..numberOfDirection-1).each{ dirList[it] = (it+0.5)*angleRangeRad} + dirList.each{k, v -> + // Indicator name + def indicName = "FRONTAL_AREA_INDEX_H${layer_bottom}_${layer_top}_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() + THEN LENGTH*SIN($v-AZIMUTH)*($HEIGHT_WALL-$layer_bottom)/$deltaH + ELSE CASE WHEN $v - AZIMUTH < -PI() + THEN LENGTH*SIN($v+2*PI()-AZIMUTH)*($HEIGHT_WALL-$layer_bottom)/$deltaH + ELSE 0 + END + END AS $indicName""" + dirQueryDiv[k] = """COALESCE(SUM(b.$indicName)/ST_AREA(a.$GEOMETRIC_FIELD_RSU), 0) AS $indicName""" + } + // Calculates projected surfaces for buildings and shared facades + datasource """ + DROP TABLE IF EXISTS $bufferTable; + CREATE TABLE $bufferTable + AS SELECT $idRsu, + ${dirQueryVertFrac.values().join(",")} + FROM $allLinesRsu + WHERE $HEIGHT_WALL > $layer_bottom""".toString() + // Fill missing values with 0 + datasource."$bufferTable"."$idRsu".createIndex() + datasource """ + DROP TABLE IF EXISTS ${tab_H[listLayersBottom.size()-1]}; + CREATE TABLE ${tab_H[listLayersBottom.size()-1]} + AS SELECT a.$idRsu, + ${dirQueryDiv.values().join(",")} + FROM $rsuTable AS a LEFT JOIN $bufferTable AS b + ON a.$idRsu = b.$idRsu + GROUP BY a.$idRsu""".toString() + // Declare this layer to the layer to join at the end + indicToJoin.put(tab_H[listLayersBottom.size()-1], idRsu) + + // 5. Join all layers in one table + def joinGrids = Geoindicators.DataUtils.joinTables() + if (!joinGrids([inputTableNamesWithId: indicToJoin, + outputTableName : outputTableName, + datasource : datasource])) { + info "Cannot merge all indicators in RSU table $outputTableName." + return + } + + // The temporary tables are deleted + //datasource """DROP TABLE IF EXISTS $buildLine, $allLinesRsu, + // $bufferTable, ${tab_H.values().join(",")}""".toString() + } + + [outputTableName: outputTableName] + } + } +} \ No newline at end of file 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 7c6821b72c..e25969e891 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicators.groovy @@ -1851,15 +1851,21 @@ IProcess rasterizeIndicators() { gridTableName: String, list_indicators :[], buildingTable: "", roadTable: "", vegetationTable: "", hydrographicTable: "", imperviousTable: "", rsu_lcz:"", - rsu_utrf_area:"",rsu_utrf_floor_area:"", + rsu_utrf_area:"",rsu_utrf_floor_area:"", seaLandMaskTableName: "", prefixName: String outputs outputTableName: String run { datasource, gridTableName, list_indicators,buildingTable, roadTable, vegetationTable, - hydrographicTable, imperviousTable, rsu_lcz,rsu_utrf_area,rsu_utrf_floor_area, prefixName -> + hydrographicTable, imperviousTable, rsu_lcz,rsu_utrf_area,rsu_utrf_floor_area, seaLandMaskTableName, + prefixName -> if(!list_indicators){ info "The list of indicator names cannot be null or empty" return } + // Temporary (and output tables) are created + def tesselatedSeaLandTab = postfix "TESSELATED_SEA_LAND" + def seaLandFractionTab = postfix "SEA_LAND_FRACTION" + + def seaLandTypeField = "TYPE" def grid_indicators_table = "grid_indicators" def grid_column_identifier ="id" def indicatorTablesToJoin = [:] @@ -1956,13 +1962,14 @@ IProcess rasterizeIndicators() { // If any surface fraction calculation is needed, create the priority list containing only needed fractions // and also set which type of statistics is needed if "BUILDING_HEIGHT" is activated + def surfaceFractionsProcess def columnFractionsList = [:] def priorities = ["water", "building", "high_vegetation", "low_vegetation", "road", "impervious"] def unweightedBuildingIndicators = [:] def weightedBuildingIndicators = [:] list_indicators.each{ - if(it.equalsIgnoreCase("BUILDING_FRACTION")){ + if(it.equalsIgnoreCase("BUILDING_FRACTION") || it.equalsIgnoreCase("BUILDING_SURFACE_DENSITY")){ columnFractionsList.put( priorities.indexOf("building"),"building") } else if(it.equalsIgnoreCase("WATER_FRACTION")){ @@ -1998,7 +2005,7 @@ IProcess rasterizeIndicators() { imperviousTable: imperviousTable, prefixName : prefixName, datasource: datasource])) { def superpositionsTableGrid = computeSmallestGeom.results.outputTableName - def surfaceFractionsProcess = Geoindicators.RsuIndicators.surfaceFractions() + surfaceFractionsProcess = Geoindicators.RsuIndicators.surfaceFractions() def superpositions = [] if(surfaceFractionsProcess.execute([ rsuTable: gridTableName, spatialRelationsTable: superpositionsTableGrid, @@ -2098,7 +2105,8 @@ IProcess rasterizeIndicators() { } } - if(list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY") && buildingTable){ + if((list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY") && buildingTable) || + (list_indicators*.toUpperCase().contains("BUILDING_SURFACE_DENSITY") && buildingTable)){ if(!createScalesRelationsGridBl){ // Create the relations between grid cells and buildings createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin() @@ -2119,12 +2127,127 @@ IProcess rasterizeIndicators() { idRsu : grid_column_identifier, prefixName : prefixName, datasource : datasource])) { - indicatorTablesToJoin.put(freeFacadeDensityExact.results.outputTableName, grid_column_identifier) + if(list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY")){ + indicatorTablesToJoin.put(freeFacadeDensityExact.results.outputTableName, grid_column_identifier) + } + if(list_indicators*.toUpperCase().contains("FREE_EXTERNAL_FACADE_DENSITY")){ + def buildingSurfDensity = Geoindicators.RsuIndicators.buildingSurfaceDensity() + if (buildingSurfDensity.execute([ + facadeDensityTable : freeFacadeDensityExact.results.outputTableName, + buildingFractionTable : surfaceFractionsProcess.results.outputTableName, + facDensityColumn : "FREE_EXTERNAL_FACADE_DENSITY", + buFractionColumn : "building_fraction", + idRsu : grid_column_identifier, + prefixName : prefixName, + datasource : datasource])) { + if (list_indicators*.toUpperCase().contains("BUILDING_SURFACE_DENSITY")) { + indicatorTablesToJoin.put(buildingSurfDensity.results.outputTableName, grid_column_identifier) + } + } + } } else { info "Cannot calculate the exact free external facade density" } } + if(list_indicators*.toUpperCase().contains("BUILDING_HEIGHT_DIST") && buildingTable){ + if(!createScalesRelationsGridBl){ + // Create the relations between grid cells and buildings + createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin() + if (!createScalesRelationsGridBl([datasource : datasource, + sourceTable : buildingTable, + targetTable : gridTableName, + idColumnTarget : grid_column_identifier, + prefixName : prefixName, + nbRelations : null])) { + info "Cannot compute the scales relations between buildings and grid cells." + return + } + } + def roofFractionDistributionExact = Geoindicators.RsuIndicators.roofFractionDistributionExact() + if (roofFractionDistributionExact( + [buildingTable : createScalesRelationsGridBl.results.outputTableName, + rsuTable : gridTableName, + idRsu : grid_column_identifier, + listLayersBottom : [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], + prefixName : prefixName, + datasource : datasource])) { + indicatorTablesToJoin.put(roofFractionDistributionExact.results.outputTableName, grid_column_identifier) + } else { + info "Cannot compute the roof fraction distribution." + } + } + + if(list_indicators*.toUpperCase().contains("FRONTAL_AREA_INDEX") && buildingTable){ + if(!createScalesRelationsGridBl){ + // Create the relations between grid cells and buildings + createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin() + if (!createScalesRelationsGridBl([datasource : datasource, + sourceTable : buildingTable, + targetTable : gridTableName, + idColumnTarget : grid_column_identifier, + prefixName : prefixName, + nbRelations : null])) { + info "Cannot compute the scales relations between buildings and grid cells." + return + } + } + def frontalAreaIndexDistribution = Geoindicators.RsuIndicators.frontalAreaIndexDistribution() + if (frontalAreaIndexDistribution( + [buildingTable : createScalesRelationsGridBl.results.outputTableName, + rsuTable : gridTableName, + idRsu : grid_column_identifier, + listLayersBottom : [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], + numberOfDirection : 12, + prefixName : prefixName, + datasource : datasource])) { + indicatorTablesToJoin.put(frontalAreaIndexDistribution.results.outputTableName, grid_column_identifier) + } else { + info "Cannot compute the frontal area index." + } + } + + if(list_indicators*.toUpperCase().contains("SEA_LAND_FRACTION") && seaLandMaskTableName) { + // 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 $seaLandMaskTableName + GROUP BY $seaLandTypeField""") + if (! sea_land_type_rows[seaLandTypeField].get("SEA")) { + datasource """ + DROP TABLE IF EXISTS $seaLandFractionTab; + CREATE TABLE $seaLandFractionTab + AS SELECT $grid_column_identifier, 1 AS LAND_FRACTION + FROM $grid_indicators_table""" + indicatorTablesToJoin.put(seaLandFractionTab, grid_column_identifier) + } else { + // Split the potentially big complex seaLand geometries into smaller triangles in order to makes calculation faster + datasource """ + DROP TABLE IF EXISTS $tesselatedSeaLandTab; + CREATE TABLE $tesselatedSeaLandTab(id_tesselate serial, the_geom geometry, $seaLandTypeField VARCHAR) + AS SELECT null, the_geom, $seaLandTypeField + FROM ST_EXPLODE('select st_tesselate(the_geom) AS the_geom, $seaLandTypeField FROM $seaLandMaskTableName')""" + + def upperScaleAreaStatistics = Geoindicators.GenericIndicators.upperScaleAreaStatistics() + if (upperScaleAreaStatistics( + [upperTableName: grid_indicators_table, + upperColumnId: grid_column_identifier, + lowerTableName: tesselatedSeaLandTab, + lowerColumnName: seaLandTypeField, + prefixName: prefixName, + datasource: datasource])) { + // Modify columns name to postfix with "_FRACTION" + datasource """ + ALTER TABLE ${upperScaleAreaStatistics.results.outputTableName} RENAME COLUMN LAND TO LAND_FRACTION; + ALTER TABLE ${ + upperScaleAreaStatistics.results.outputTableName + } RENAME COLUMN SEA TO SEA_FRACTION;""" + indicatorTablesToJoin.put(upperScaleAreaStatistics.results.outputTableName, grid_column_identifier) + } else { + info "Cannot compute the frontal area index." + } + } + } + //Join all indicators at grid scale def joinGrids = Geoindicators.DataUtils.joinTables() if (!joinGrids([inputTableNamesWithId: indicatorTablesToJoin, @@ -2133,6 +2256,9 @@ IProcess rasterizeIndicators() { info "Cannot merge all indicators in grid table $grid_indicators_table." return } + + // Remove temporary tables + datasource """DROP TABLE IF EXISTS $tesselatedSeaLandTab, $seaLandFractionTab""" [outputTableName: grid_indicators_table] } } 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 a17433268f..14335052f8 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/RsuIndicatorsTests.groovy @@ -641,4 +641,119 @@ class RsuIndicatorsTests { assertEquals(3.0/10, result0["building_fraction"]) assertEquals(0d, result0["undefined_fraction"]) } + + @Test + void buildingSurfaceDensityTest() { + h2GIS """ + DROP TABLE IF EXISTS facade_density_tab, building_fraction_tab; + CREATE TABLE facade_density_tab(id_rsu int, facade_density double); + INSERT INTO facade_density_tab VALUES (1, 1.2), + (2, 0); + CREATE TABLE building_fraction_tab(id_rsu int, building_fraction double); + INSERT INTO building_fraction_tab VALUES (1, 0.5), + (2, 1); + """ + + def p0 = Geoindicators.RsuIndicators.buildingSurfaceDensity() + assertTrue p0.execute([ + facadeDensityTable: "facade_density_tab", buildingFractionTable: "building_fraction_tab", + facDensityColumn: "facade_density", + buFractionColumn: "building_fraction", + idRsu: "id_rsu", + prefixName: "test", datasource: h2GIS]) + def result1 = h2GIS.firstRow("SELECT * FROM ${p0.results.outputTableName} WHERE id_rsu=1") + assertEquals(1.7, result1["building_surface_density"]) + def result2 = h2GIS.firstRow("SELECT * FROM ${p0.results.outputTableName} WHERE id_rsu=2") + assertEquals(1.0, result2["building_surface_density"]) + } + + @Test + void roofFractionDistributionExactTest() { + // Only the first 1 first created buildings are selected for the tests + h2GIS """ + DROP TABLE IF EXISTS tempo_build, tempo_rsu; + CREATE TABLE tempo_build(id_build int, the_geom geometry, height_wall double, height_roof double); + INSERT INTO tempo_build VALUES (1, 'POLYGON((-50 -50, 50 -50, 50 50, -50 50, -50 -50))'::GEOMETRY, 2, 4), + (2, 'POLYGON((50 -50, 150 -50, 150 50, 50 50, 50 -50))'::GEOMETRY, 18, 24), + (3, 'POLYGON((50 50, 100 50, 100 150, 50 150, 50 50))'::GEOMETRY, 60, 60); + 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), + (2, 'POLYGON((100 0, 200 0, 200 100, 100 100, 100 0))'::GEOMETRY), + (3, 'POLYGON((0 100, 100 100, 100 200, 0 200, 0 100))'::GEOMETRY), + (4, 'POLYGON((100 100, 200 100, 200 200, 100 200, 100 100))'::GEOMETRY); + """ + // First calculate the correlation table between buildings and rsu + def createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin() + createScalesRelationsGridBl([datasource : h2GIS, + sourceTable : "tempo_build", + targetTable : "tempo_rsu", + idColumnTarget : "id_rsu", + prefixName : "test", + nbRelations : null]) + + def p = Geoindicators.RsuIndicators.roofFractionDistributionExact() + assertTrue p([ + buildingTable : createScalesRelationsGridBl.results.outputTableName, + rsuTable : "tempo_rsu", + idRsu : "id_rsu", + prefixName : "test", + listLayersBottom : [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], + datasource : h2GIS]) + assertEquals 1.0/3, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 1").ROOF_FRACTION_DISTRIBUTION_0_5, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 1").ROOF_FRACTION_DISTRIBUTION_15_20, 0.00001 + assertEquals 1.0/3, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 1").ROOF_FRACTION_DISTRIBUTION_20_25, 0.00001 + assertEquals 1.0/3, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 1").ROOF_FRACTION_DISTRIBUTION_50_INF, 0.00001 + assertEquals 1.0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 2").ROOF_FRACTION_DISTRIBUTION_20_25, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 2").ROOF_FRACTION_DISTRIBUTION_0_5, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 2").ROOF_FRACTION_DISTRIBUTION_50_INF, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 3").ROOF_FRACTION_DISTRIBUTION_20_25, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 3").ROOF_FRACTION_DISTRIBUTION_0_5, 0.00001 + assertEquals 1.0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 3").ROOF_FRACTION_DISTRIBUTION_50_INF, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 4").ROOF_FRACTION_DISTRIBUTION_20_25, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 4").ROOF_FRACTION_DISTRIBUTION_0_5, 0.00001 + assertEquals 0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 4").ROOF_FRACTION_DISTRIBUTION_50_INF, 0.00001 + } + + @Test + void frontalAreaIndexDistributionTest() { + // Only the first 1 first created buildings are selected for the tests + h2GIS """ + DROP TABLE IF EXISTS tempo_build, tempo_rsu; + CREATE TABLE tempo_build(id_build int, the_geom geometry, height_wall double); + INSERT INTO tempo_build VALUES (1, 'POLYGON((-50 -50, 50 -50, 50 50, -50 50, -50 -50))'::GEOMETRY, 3), + (2, 'POLYGON((50 -50, 150 -50, 150 50, 50 50, 50 -50))'::GEOMETRY, 21), + (3, 'POLYGON((50 50, 100 50, 100 150, 50 150, 50 50))'::GEOMETRY, 60); + 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), + (2, 'POLYGON((100 0, 200 0, 200 100, 100 100, 100 0))'::GEOMETRY), + (3, 'POLYGON((0 100, 100 100, 100 200, 0 200, 0 100))'::GEOMETRY), + (4, 'POLYGON((100 100, 200 100, 200 200, 100 200, 100 100))'::GEOMETRY), + (5, 'POLYGON((200 200, 300 200, 300 300, 200 300, 200 200))'::GEOMETRY); + """ + // First calculate the correlation table between buildings and rsu + def createScalesRelationsGridBl = Geoindicators.SpatialUnits.spatialJoin() + createScalesRelationsGridBl([datasource : h2GIS, + sourceTable : "tempo_build", + targetTable : "tempo_rsu", + idColumnTarget : "id_rsu", + prefixName : "test", + nbRelations : null]) + + + def p = Geoindicators.RsuIndicators.frontalAreaIndexDistribution() + assertTrue p([ + buildingTable : createScalesRelationsGridBl.results.outputTableName, + rsuTable : "tempo_rsu", + idRsu : "id_rsu", + listLayersBottom : [0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50], + numberOfDirection : 12, + prefixName : "test", + datasource : h2GIS]) + assertEquals 0.00566, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 1").FRONTAL_AREA_INDEX_H0_5_D30_60, 0.00001 + assertEquals 0.00321, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 1").FRONTAL_AREA_INDEX_H50_61_D30_60, 0.00001 + assertEquals 0.00321, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 4").FRONTAL_AREA_INDEX_H50_61_D30_60, 0.00001 + assertEquals 0.0, h2GIS.firstRow("SELECT * FROM ${p.results.outputTableName} WHERE id_rsu = 5").FRONTAL_AREA_INDEX_H0_5_D30_60, 0.00001 + } + } + 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 7ae6c888a6..e30f98a727 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy @@ -1106,8 +1106,8 @@ def extractProcessingParameters(def processing_parameters){ return } def allowed_grid_indicators=["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION","WATER_FRACTION","VEGETATION_FRACTION", - "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_AREA_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY", "FREE_EXTERNAL_FACADE_DENSITY", - "BUILDING_HEIGHT_WEIGHTED"] + "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_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"] def allowedOutputIndicators = allowed_grid_indicators.intersect(list_indicators*.toUpperCase()) if(allowedOutputIndicators){ //Update the RSU indicators list according the grid indicators 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 f6370a94ff..4562487c4e 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -582,8 +582,10 @@ class WorflowOSMTest extends WorkflowAbstractTest { //"rowCol": true, "indicators": ["BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP", "BUILDING_TYPE_FRACTION","WATER_FRACTION","VEGETATION_FRACTION", - "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "LCZ_FRACTION", "LCZ_PRIMARY", - "FREE_EXTERNAL_FACADE_DENSITY", "BUILDING_HEIGHT_WEIGHTED"] + "ROAD_FRACTION", "IMPERVIOUS_FRACTION", "UTRF_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"] ], "worldpop_indicators" : true ] ]