From 0f12cebf7a6978b919836f67761b098687d16e76 Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 19 Mar 2024 11:24:43 +0100 Subject: [PATCH 1/3] Improve OSM data extraction and OSM data formating --- .../geoclimate/geoindicators/DataUtils.groovy | 20 +-- .../geoindicators/SpatialUnits.groovy | 18 ++- .../geoclimate/osm/InputDataFormatting.groovy | 147 +++++++++++++----- .../geoclimate/osm/WorkflowOSM.groovy | 97 +++++++----- .../osm/InputDataFormattingTest.groovy | 15 +- .../geoclimate/osm/WorflowOSMTest.groovy | 11 +- 6 files changed, 208 insertions(+), 100 deletions(-) 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 f4987d5e36..3ada827c0e 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/DataUtils.groovy @@ -89,14 +89,16 @@ String joinTables(JdbcDataSource datasource, Map inputTableNamesWithId, String o /** * An utility process to save several tables in a folder * + * + * @param datasource connection to the database * @param inputTableNames to be stored in the directory. - * Note : A spatial table is saved in a geojson file and the other in csv + * Note : A spatial table is saved in a flatgeobuffer file and the other in csv + * @param delete true to delete the file is exist * @param directory folder to save the tables - * @param datasource connection to the database * * @return the directory where the tables are saved */ -String saveTablesAsFiles(JdbcDataSource datasource, List inputTableNames, boolean delete = false, String directory) { +String saveTablesAsFiles(JdbcDataSource datasource, List inputTableNames, boolean delete = true, String directory) { if (directory == null) { error "The directory to save the data cannot be null" return @@ -113,7 +115,7 @@ String saveTablesAsFiles(JdbcDataSource datasource, List inputTableNames, boolea inputTableNames.each { tableName -> if (tableName) { def fileToSave = dirFile.absolutePath + File.separator + tableName + - (datasource."$tableName".spatial ? ".geojson" : ".csv") + (datasource."$tableName".spatial ? ".fgb" : ".csv") def table = datasource.getTable(tableName) if (table) { table.save(fileToSave, delete) @@ -164,10 +166,10 @@ static String aliasColumns(JdbcDataSource datasource, String tableName, String a /** * Create the select projection and alias all columns - * @param datasource - * @param tableName - * @param alias - * @param exceptColumns + * @param datasource connection to the database + * @param tableName table name + * @param alias for the columns + * @param exceptColumns columns to remove * @return */ static String aliasColumns(JdbcDataSource datasource, def tableName, def alias, def exceptColumns){ @@ -176,4 +178,4 @@ static String aliasColumns(JdbcDataSource datasource, def tableName, def alias, return columnNames.inject([]) { result, iter -> result += "$alias.$iter" }.join(",") -} +} \ No newline at end of file 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 b5e29b27e2..ddea9619e8 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -123,17 +123,23 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz return null } if (inputzone) { - datasource.createSpatialIndex(inputTableName, "the_geom") + String polygons = postfix("polygons") + datasource.execute("""DROP TABLE IF EXISTS $polygons; + CREATE TABLE $polygons as 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 + FROM $inputTableName)') + """.toString()) + datasource.createSpatialIndex(polygons) + datasource.createSpatialIndex(inputTableName) 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 - FROM ST_EXPLODE('( - SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom - FROM $inputTableName)') AS a, + SELECT a.* + FROM $polygons AS a, $inputzone AS b WHERE a.the_geom && b.the_geom - AND ST_INTERSECTS(ST_POINTONSURFACE(a.THE_GEOM), b.the_geom) and st_area(a.the_geom) > $area + AND ST_INTERSECTS(ST_POINTONSURFACE(a.THE_GEOM), b.the_geom) and st_area(a.the_geom) > $area; + DROP TABLE IF EXISTS $polygons; """.toString() } else { datasource """ diff --git a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy index e9e9c62708..78dbf968ce 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy @@ -23,6 +23,8 @@ import groovy.json.JsonSlurper import groovy.transform.BaseScript import org.locationtech.jts.geom.Geometry import org.locationtech.jts.geom.Polygon +import org.locationtech.jts.geom.prep.PreparedGeometry +import org.locationtech.jts.geom.prep.PreparedGeometryFactory import org.orbisgis.data.jdbc.JdbcDataSource import org.orbisgis.geoclimate.Geoindicators @@ -72,44 +74,39 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone def columnToMap = parametersMap.columns def inputSpatialTable = datasource."$building" if (inputSpatialTable.rowCount > 0) { + def heightPattern = Pattern.compile("((?:\\d+\\/|(?:\\d+|^|\\s)\\.)?\\d+)\\s*([^\\s\\d+\\-.,:;^\\/]+(?:\\^\\d+(?:\$|(?=[\\s:;\\/])))?(?:\\/[^\\s\\d+\\-.,:;^\\/]+(?:\\^\\d+(?:\$|(?=[\\s:;\\/])))?)*)?", Pattern.CASE_INSENSITIVE) def columnNames = inputSpatialTable.columns columnNames.remove("THE_GEOM") queryMapper += columnsMapper(columnNames, columnToMap) + queryMapper += " , st_force2D(a.the_geom) as the_geom FROM $building as a " if (zone) { - datasource.createSpatialIndex(building, "the_geom") - datasource.createSpatialIndex(zone, "the_geom") - queryMapper += " , st_force2D(a.the_geom) as the_geom FROM $building as a, $zone as b WHERE a.the_geom && b.the_geom and st_intersects( " + - "a.the_geom, b.the_geom) " - } else { - queryMapper += " , st_force2D(a.the_geom) as the_geom FROM $building as a " - } - - def heightPattern = Pattern.compile("((?:\\d+\\/|(?:\\d+|^|\\s)\\.)?\\d+)\\s*([^\\s\\d+\\-.,:;^\\/]+(?:\\^\\d+(?:\$|(?=[\\s:;\\/])))?(?:\\/[^\\s\\d+\\-.,:;^\\/]+(?:\\^\\d+(?:\$|(?=[\\s:;\\/])))?)*)?", Pattern.CASE_INSENSITIVE) - def id_build = 1; - datasource.withBatch(100) { stmt -> - datasource.eachRow(queryMapper) { row -> - def typeAndUseValues = getTypeAndUse(row, columnNames, mappingTypeAndUse) - def use = typeAndUseValues[1] - def type = typeAndUseValues[0] - if (type) { - String height = row.height - String roof_height = row.'roof:height' - String b_lev = row.'building:levels' - String roof_lev = row.'roof:levels' - def heightRoof = getHeightRoof(height, heightPattern) - def heightWall = getHeightWall(heightRoof, roof_height) - def nbLevels = getNbLevels(b_lev, roof_lev) - def formatedHeight = Geoindicators.WorkflowGeoIndicators.formatHeightsAndNbLevels(heightWall, heightRoof, nbLevels, h_lev_min, type, typeAndLevel) - def zIndex = getZIndex(row.'layer') - String roof_shape = row.'roof:shape' - - if (formatedHeight.nbLevels > 0 && zIndex >= 0 && type) { - Geometry geom = row.the_geom - def srid = geom.getSRID() - for (int i = 0; i < geom.getNumGeometries(); i++) { - Geometry subGeom = geom.getGeometryN(i) - if (subGeom instanceof Polygon && subGeom.getArea() > 1) { - stmt.addBatch """ + Geometry geomZone = datasource.firstRow("select st_union(st_accum(the_geom)) as the_geom from $zone").the_geom + PreparedGeometry pZone = PreparedGeometryFactory.prepare(geomZone) + def id_build = 1; + datasource.withBatch(200) { stmt -> + datasource.eachRow(queryMapper) { row -> + def typeAndUseValues = getTypeAndUse(row, columnNames, mappingTypeAndUse) + def use = typeAndUseValues[1] + def type = typeAndUseValues[0] + if (type) { + String height = row.height + String roof_height = row.'roof:height' + String b_lev = row.'building:levels' + String roof_lev = row.'roof:levels' + def heightRoof = getHeightRoof(height, heightPattern) + def heightWall = getHeightWall(heightRoof, roof_height) + def nbLevels = getNbLevels(b_lev, roof_lev) + def formatedHeight = Geoindicators.WorkflowGeoIndicators.formatHeightsAndNbLevels(heightWall, heightRoof, nbLevels, h_lev_min, type, typeAndLevel) + def zIndex = getZIndex(row.'layer') + String roof_shape = row.'roof:shape' + if (formatedHeight.nbLevels > 0 && zIndex >= 0 && type) { + Geometry geom = row.the_geom + if(pZone.intersects(geom)){ + def srid = geom.getSRID() + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i) + if (subGeom instanceof Polygon && subGeom.getArea() > 1) { + stmt.addBatch """ INSERT INTO ${outputTableName} values( ST_GEOMFROMTEXT('${subGeom}',$srid), $id_build, @@ -123,31 +120,96 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone ${roof_shape ? "'" + roof_shape + "'" : null}) """.toString() - if (formatedHeight.estimated) { + if (formatedHeight.estimated) { + stmt.addBatch """ + INSERT INTO ${outputEstimateTableName} values( + $id_build, + '${row.id}') + """.toString() + } + id_build++ + } + } + } + } + } + } + } + } else { + def id_build = 1; + datasource.withBatch(200) { stmt -> + datasource.eachRow(queryMapper) { row -> + def typeAndUseValues = getTypeAndUse(row, columnNames, mappingTypeAndUse) + def use = typeAndUseValues[1] + def type = typeAndUseValues[0] + if (type) { + String height = row.height + String roof_height = row.'roof:height' + String b_lev = row.'building:levels' + String roof_lev = row.'roof:levels' + def heightRoof = getHeightRoof(height, heightPattern) + def heightWall = getHeightWall(heightRoof, roof_height) + def nbLevels = getNbLevels(b_lev, roof_lev) + def formatedHeight = Geoindicators.WorkflowGeoIndicators.formatHeightsAndNbLevels(heightWall, heightRoof, nbLevels, h_lev_min, type, typeAndLevel) + def zIndex = getZIndex(row.'layer') + String roof_shape = row.'roof:shape' + + if (formatedHeight.nbLevels > 0 && zIndex >= 0 && type) { + Geometry geom = row.the_geom + def srid = geom.getSRID() + for (int i = 0; i < geom.getNumGeometries(); i++) { + Geometry subGeom = geom.getGeometryN(i) + if (subGeom instanceof Polygon && subGeom.getArea() > 1) { stmt.addBatch """ + INSERT INTO ${outputTableName} values( + ST_GEOMFROMTEXT('${subGeom}',$srid), + $id_build, + '${row.id}', + ${formatedHeight.heightWall}, + ${formatedHeight.heightRoof}, + ${formatedHeight.nbLevels}, + ${singleQuote(type)}, + ${singleQuote(use)}, + ${zIndex}, + ${roof_shape ? "'" + roof_shape + "'" : null}) + """.toString() + + if (formatedHeight.estimated) { + stmt.addBatch """ INSERT INTO ${outputEstimateTableName} values( $id_build, '${row.id}') """.toString() + } + id_build++ } - id_build++ } } } } } } + + //Improve building type using the urban areas table if (urban_areas) { datasource.createSpatialIndex(outputTableName, "the_geom") datasource.createIndex(outputTableName, "id_build") datasource.createIndex(outputTableName, "type") - datasource.createSpatialIndex(urban_areas, "the_geom") - def urbanAreasPart = "urbanAreasPart${UUID.randomUUID().toString().replaceAll("-", "_")}" - def buildinType = "BUILDING_TYPE_${UUID.randomUUID().toString().replaceAll("-", "_")}" + def urbanAreasPart = postfix("urbanAreasPart") + def buildinType = postfix("BUILDING_TYPE") + String triangles = postfix("triangles") + + //Tessellate the urban_areas to perform the intersection on large geometry + datasource.execute("""DROP TABLE IF EXISTS $triangles; + CREATE TABLE $triangles as + SELECT * FROM st_explode('(SELECT CASE WHEN ST_AREA(THE_GEOM) > 100000 + THEN ST_Tessellate(st_buffer(the_geom, -0.001)) ELSE THE_GEOM END AS THE_GEOM, + type FROM $urban_areas)')""".toString()) + datasource.createSpatialIndex(triangles) datasource.execute """DROP TABLE IF EXISTS $urbanAreasPart, $buildinType ; - CREATE TABLE $urbanAreasPart as SELECT b.type, a.id_build, st_area(st_intersection(a.the_geom,b.the_geom))/st_area(a.the_geom) as part FROM $outputTableName a, $urban_areas b + CREATE TABLE $urbanAreasPart as SELECT b.type, a.id_build, st_area(st_intersection(a.the_geom,b.the_geom))/st_area(a.the_geom) as part FROM $outputTableName a, $triangles b WHERE a.the_geom && b.the_geom and st_intersects(a.the_geom, b.the_geom) AND a.TYPE ='building'; CREATE INDEX ON $urbanAreasPart(id_build); create table $buildinType as select @@ -201,7 +263,7 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone stmt.addBatch("UPDATE $outputTableName set type = '${new_type}' , main_use = '${new_type}'where id_build=${id_build_}") } } - datasource.execute("""DROP TABLE IF EXISTS $urbanAreasPart, $buildinType""".toString()) + datasource.execute("""DROP TABLE IF EXISTS $urbanAreasPart, $buildinType, $triangles""".toString()) } } @@ -1012,7 +1074,8 @@ String formatUrbanAreas(JdbcDataSource datasource, String urban_areas, String zo } } } - //Merging landuse + //Merging urban_areas + def mergingUrbanAreas = postfix("merging_urban_areas") datasource.execute(""" CREATE TABLE $mergingUrbanAreas as select CAST((row_number() over()) as Integer) as id_urban,the_geom, type from 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 60fe0fcf95..bde2048844 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy @@ -417,13 +417,13 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d def zones = extractOSMZone(h2gis_datasource, id_zone, processing_parameters.distance, bbox_size) if (zones) { id_zone = id_zone in Collection ? id_zone.join('_') : id_zone - def zone = zones.zone - def zoneEnvelopeTableName = zones.zone_envelope - if (h2gis_datasource.getRowCount(zone) == 0) { + def utm_zone_table = zones.utm_zone_table + def utm_extended_bbox_table = zones.utm_extended_bbox_table + if (h2gis_datasource.getRowCount(utm_zone_table) == 0) { error "Cannot find any geometry to define the zone to extract the OSM data" return } - def srid = h2gis_datasource.getSrid(zone) + def srid = zones.utm_srid def reproject = false if (outputSRID) { if (outputSRID != srid) { @@ -432,13 +432,13 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d } else { outputSRID = srid } - //Prepare OSM extraction + //Prepare OSM extraction from the osm_envelope_extented //TODO set key values ? def osm_date="" if(overpass_date){ osm_date = "[date:\"$overpass_date\"]" } - def query = "[timeout:$overpass_timeout][maxsize:$overpass_maxsize]$osm_date" + OSMTools.Utilities.buildOSMQuery(zones.envelope, null, OSMElement.NODE, OSMElement.WAY, OSMElement.RELATION) + def query = "[timeout:$overpass_timeout][maxsize:$overpass_maxsize]$osm_date" + OSMTools.Utilities.buildOSMQuery(zones.osm_envelope_extented, null, OSMElement.NODE, OSMElement.WAY, OSMElement.RELATION) if (downloadAllOSMData) { //Create a custom OSM query to download all requiered data. It will take more time and resources @@ -447,14 +447,15 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d "leisure", "highway", "natural", "landuse", "landcover", "vegetation", "waterway", "area", "aeroway", "area:aeroway", "tourism", "sport"] - query = "[timeout:$overpass_timeout][maxsize:$overpass_maxsize]$osm_date" + OSMTools.Utilities.buildOSMQueryWithAllData(zones.envelope, keysValues, OSMElement.NODE, OSMElement.WAY, OSMElement.RELATION) + query = "[timeout:$overpass_timeout][maxsize:$overpass_maxsize]$osm_date" + OSMTools.Utilities.buildOSMQueryWithAllData(zones.osm_envelope_extented, keysValues, OSMElement.NODE, OSMElement.WAY, OSMElement.RELATION) } def extract = OSMTools.Loader.extract(query) if (extract) { - Geometry geomArea = h2gis_datasource.getExtent(zones.zone) - geomArea.setSRID(srid) - Map gisLayersResults = OSM.InputDataLoading.createGISLayers(h2gis_datasource, extract, geomArea, srid) + //We must build the GIS layers on the extended bbox area + Geometry utm_extented_geom = h2gis_datasource.getExtent(utm_extended_bbox_table) + utm_extented_geom.setSRID(srid) + Map gisLayersResults = OSM.InputDataLoading.createGISLayers(h2gis_datasource, extract, utm_extented_geom, srid) if (gisLayersResults) { if (deleteOSMFile) { if (new File(extract).delete()) { @@ -466,34 +467,47 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d def road_traffic = processing_parameters.road_traffic def worldpop_indicators = processing_parameters.worldpop_indicators - debug "Formating OSM GIS layers" + info "Formating OSM GIS layers" //Format urban areas - String urbanAreasTable = OSM.InputDataFormatting.formatUrbanAreas(h2gis_datasource, gisLayersResults.urban_areas, zoneEnvelopeTableName) + String urbanAreasTable = OSM.InputDataFormatting.formatUrbanAreas(h2gis_datasource, gisLayersResults.urban_areas, utm_extended_bbox_table) + info "Urban areas formatted" + /* + * Do not filter the data when formatting becausethe job is already done when extracting osm data * + */ Map formatBuilding = OSM.InputDataFormatting.formatBuildingLayer( h2gis_datasource, gisLayersResults.building, - zoneEnvelopeTableName, urbanAreasTable, + null, urbanAreasTable, processing_parameters.hLevMin) - + info "Building formatted" def buildingTableName = formatBuilding.building def buildingEstimateTableName = formatBuilding.building_estimated - String railTableName = OSM.InputDataFormatting.formatRailsLayer(h2gis_datasource, gisLayersResults.rail, zoneEnvelopeTableName) + String railTableName = OSM.InputDataFormatting.formatRailsLayer(h2gis_datasource, gisLayersResults.rail, null) + + info "Rail formatted" - String vegetationTableName = OSM.InputDataFormatting.formatVegetationLayer(h2gis_datasource, gisLayersResults.vegetation, zoneEnvelopeTableName) + String vegetationTableName = OSM.InputDataFormatting.formatVegetationLayer(h2gis_datasource, gisLayersResults.vegetation, utm_extended_bbox_table) + + info "Vegetation formatted" String hydrographicTableName = OSM.InputDataFormatting.formatWaterLayer( h2gis_datasource, gisLayersResults.water, - zoneEnvelopeTableName) + utm_extended_bbox_table) + + info "Water formatted" String imperviousTableName = OSM.InputDataFormatting.formatImperviousLayer( - h2gis_datasource, gisLayersResults.impervious, zoneEnvelopeTableName) + h2gis_datasource, gisLayersResults.impervious, utm_extended_bbox_table) + + info "Impervious formatted" //Sea/Land mask String seaLandMaskTableName = OSM.InputDataFormatting.formatSeaLandMask( - h2gis_datasource, gisLayersResults.coastline, zoneEnvelopeTableName, hydrographicTableName) + h2gis_datasource, gisLayersResults.coastline, utm_extended_bbox_table, hydrographicTableName) + info "Sea/Land formatted" if(h2gis_datasource.getRowCount(seaLandMaskTableName)>0){ //Select the water and sea features @@ -504,7 +518,11 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d //Format road String roadTableName = OSM.InputDataFormatting.formatRoadLayer( h2gis_datasource, gisLayersResults.road, - zoneEnvelopeTableName) + utm_extended_bbox_table) + + info "Road formatted" + + info "All layers have been formatted" //Drop the intermediate GIS layers h2gis_datasource.dropTable(gisLayersResults.values().toArray(new String[0])) @@ -512,7 +530,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d debug "OSM GIS layers formated" //Add the GIS layers to the list of results def results = [:] - results.put("zone", zone) + results.put("zone", utm_zone_table) results.put("road", roadTableName) results.put("rail", railTableName) results.put("water", hydrographicTableName) @@ -535,7 +553,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d rsu_indicators_params.put("utrfModelName", "UTRF_OSM_RF_2_2.model") rsu_indicators_params.put("buildingHeightModelName", estimateHeight) Map geoIndicators = Geoindicators.WorkflowGeoIndicators.computeAllGeoIndicators( - h2gis_datasource, zone, + h2gis_datasource, utm_zone_table, buildingTableName, roadTableName, railTableName, vegetationTableName, hydrographicTableName, imperviousTableName, @@ -601,7 +619,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d def geomEnv; if (noise_indicators) { if (noise_indicators.ground_acoustic) { - geomEnv = h2gis_datasource.getSpatialTable(zone).getExtent() + geomEnv = h2gis_datasource.getSpatialTable(utm_zone_table).getExtent() def outputTable = Geoindicators.SpatialUnits.createGrid(h2gis_datasource, geomEnv, 200, 200) if (outputTable) { String ground_acoustic = Geoindicators.NoiseIndicators.groundAcousticAbsorption(h2gis_datasource, outputTable, "id_grid", @@ -618,7 +636,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d def outputGrid = "geojson" if (grid_indicators_params) { if (!geomEnv) { - geomEnv = h2gis_datasource.getSpatialTable(zone).getExtent() + geomEnv = h2gis_datasource.getSpatialTable(utm_zone_table).getExtent() } outputGrid = grid_indicators_params.output def x_size = grid_indicators_params.x_size @@ -697,6 +715,11 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d * @param zoneToExtract the osm filter : place or bbox * @param distance to expand the OSM bbox * @return + * utm_zone_table the geometries that represents the processed zone in UTM coordinate system + * utm_extended_bbox_table the envelope of the zone to be processed extended by a distance in UTM coordinate system + * osm_envelope_extented the bbox used to query overpass in lat/lon + * osm_geometry the geometry that represents the processed zone in lat/lon + * utm_srid the UTM srid code */ def extractOSMZone(def datasource, def zoneToExtract, def distance, def bbox_size) { def outputZoneTable = "ZONE_${UUID.randomUUID().toString().replaceAll("-", "_")}".toString() @@ -719,38 +742,40 @@ def extractOSMZone(def datasource, def zoneToExtract, def distance, def bbox_siz /** * Create the OSM BBOX to be extracted + * Expand it with a distance in meters */ - def envelope = GeographyUtilities.expandEnvelopeByMeters(geom.getEnvelopeInternal(), distance) + def lat_lon_bbox_extended = GeographyUtilities.expandEnvelopeByMeters(geom.getEnvelopeInternal(), distance) //Find the best utm zone //Reproject the geometry and its envelope to the UTM zone def con = datasource.getConnection() - def interiorPoint = envelope.centre() + def interiorPoint = lat_lon_bbox_extended.centre() def epsg = GeographyUtilities.getSRID(con, interiorPoint.y as float, interiorPoint.x as float) - def geomUTM = ST_Transform.ST_Transform(con, geom, epsg) + def source_geom_utm = ST_Transform.ST_Transform(con, geom, epsg) //Check the size of the bbox - if ((geomUTM.getEnvelopeInternal().getArea() / 1.0E+6) >= bbox_size) { + if ((source_geom_utm.getEnvelopeInternal().getArea() / 1.0E+6) >= bbox_size) { error("The size of the OSM BBOX is greated than the limit : ${bbox_size} in km².\n" + "Please increase the area parameter if you want to skip this limit.") return null } - def tmpGeomEnv = geom.getFactory().toGeometry(envelope) - tmpGeomEnv.setSRID(4326) + def lat_lon_bbox_geom_extended = geom.getFactory().toGeometry(lat_lon_bbox_extended) + lat_lon_bbox_geom_extended.setSRID(4326) datasource.execute """drop table if exists ${outputZoneTable}; create table ${outputZoneTable} (the_geom GEOMETRY(${GEOMETRY_TYPE}, $epsg), ID_ZONE VARCHAR); - INSERT INTO ${outputZoneTable} VALUES (ST_GEOMFROMTEXT('${geomUTM.toString()}', ${epsg}), '${zoneToExtract.toString()}');""".toString() + INSERT INTO ${outputZoneTable} VALUES (ST_GEOMFROMTEXT('${source_geom_utm.toString()}', ${epsg}), '${zoneToExtract.toString()}');""".toString() datasource.execute """drop table if exists ${outputZoneEnvelopeTable}; create table ${outputZoneEnvelopeTable} (the_geom GEOMETRY(POLYGON, $epsg), ID_ZONE VARCHAR); - INSERT INTO ${outputZoneEnvelopeTable} VALUES (ST_GEOMFROMTEXT('${ST_Transform.ST_Transform(con, tmpGeomEnv, epsg).toString()}',${epsg}), '${zoneToExtract.toString()}'); + INSERT INTO ${outputZoneEnvelopeTable} VALUES (ST_GEOMFROMTEXT('${ST_Transform.ST_Transform(con, lat_lon_bbox_geom_extended, epsg).toString()}',${epsg}), '${zoneToExtract.toString()}'); """.toString() - return ["zone" : outputZoneTable, - "zone_envelope": outputZoneEnvelopeTable, - "envelope" : envelope, - "geometry" : geom + return ["utm_zone_table" : outputZoneTable, + "utm_extended_bbox_table": outputZoneEnvelopeTable, + "osm_envelope_extented" : lat_lon_bbox_extended, + "osm_geometry" : geom, + "utm_srid" : epsg ] } else { error "The zone to extract cannot be null or empty" diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy index 401dd1540d..4745e656b3 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy @@ -287,7 +287,7 @@ class InputDataFormattingTest { } def h2GIS = H2GIS.open("${file.absolutePath + File.separator}osm_gislayers;AUTO_SERVER=TRUE".toString()) - def zoneToExtract = "Redon" + def zoneToExtract = "Marseille" Map extractData = OSM.InputDataLoading.extractAndCreateGISLayers(h2GIS, zoneToExtract) @@ -310,27 +310,34 @@ class InputDataFormattingTest { extractData.urban_areas,extractData.zone) h2GIS.save(inputUrbanAreas,"${file.absolutePath + File.separator}osm_urban_areas_${formatedPlaceName}.geojson", true) + println("Urban areas formatted") + //Buildings h2GIS.save(extractData.building,"${file.absolutePath + File.separator}building_${formatedPlaceName}.geojson", true) def inputBuildings = OSM.InputDataFormatting.formatBuildingLayer(h2GIS, - extractData.building,extractData.zone,inputUrbanAreas) + extractData.building,extractData.zone,null) h2GIS.save(inputBuildings.building,"${file.absolutePath + File.separator}osm_building_${formatedPlaceName}.geojson", true) + println("Building formatted") //Roads def inputRoadTableName = OSM.InputDataFormatting.formatRoadLayer( h2GIS,extractData.road, extractData.zone_envelope) h2GIS.save(inputRoadTableName,"${file.absolutePath + File.separator}osm_road_${formatedPlaceName}.geojson", true) + println("Road formatted") + //Rails def inputRailTableName = OSM.InputDataFormatting.formatRailsLayer( h2GIS,extractData.rail, extractData.zone_envelope) h2GIS.save(inputRailTableName,"${file.absolutePath + File.separator}osm_rail_${formatedPlaceName}.geojson", true) + println("Rail formatted") //Vegetation def inputVegetationTableName = OSM.InputDataFormatting.formatVegetationLayer( h2GIS,extractData.vegetation,extractData.zone_envelope) h2GIS.save(inputVegetationTableName,"${file.absolutePath + File.separator}osm_vegetation_${formatedPlaceName}.geojson", true) + println("Vegetation formatted") //Hydrography def inputWaterTableName = OSM.InputDataFormatting.formatWaterLayer(h2GIS, extractData.water, extractData.zone_envelope) @@ -340,6 +347,8 @@ class InputDataFormattingTest { extractData.zone_envelope) h2GIS.save(imperviousTable,"${file.absolutePath + File.separator}osm_impervious_${formatedPlaceName}.geojson", true) + println("Impervious formatted") + //Save coastlines to debug h2GIS.save(extractData.coastline,"${file.absolutePath + File.separator}osm_coastlines_${formatedPlaceName}.geojson", true) @@ -349,6 +358,8 @@ class InputDataFormattingTest { extractData.zone_envelope, inputWaterTableName) h2GIS.save(inputSeaLandTableName,"${file.absolutePath + File.separator}osm_sea_land_${formatedPlaceName}.geojson", true) + println("Sea land mask formatted") + //Save it after sea/land mask because the water table can be modified h2GIS.save(inputWaterTableName,"${file.absolutePath + File.separator}osm_water_${formatedPlaceName}.geojson", true) 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 4b10af7d9c..ae3f365de0 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -649,11 +649,11 @@ class WorflowOSMTest extends WorkflowAbstractTest { File dirFile = new File(directory) dirFile.delete() dirFile.mkdir() - def location = "Dijon" + def location = "Marseille" def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location) def grid_size = 100 - location = nominatim.bbox - //location=[47.4, -4.8, 47.6, -4.6] + //location = nominatim.bbox + //location=[43.214935,5.336351,43.244890,5.383558] def osm_parmeters = [ "description" : "Example of configuration file to run the OSM workflow and store the result in a folder", "geoclimatedb": [ @@ -677,7 +677,7 @@ class WorflowOSMTest extends WorkflowAbstractTest { "indicatorUse": ["LCZ", "UTRF", "TEB"] - ],"grid_indicators": [ + ],/*"grid_indicators": [ "x_size": grid_size, "y_size": grid_size, //"rowCol": true, @@ -702,7 +702,8 @@ class WorflowOSMTest extends WorkflowAbstractTest { H2GIS h2gis = H2GIS.open("${directory + File.separator}geoclimate_test_integration;AUTO_SERVER=TRUE") def tableNames = results.values() def gridTable = tableNames.grid_indicators[0] - String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2gis, gridTable,grid_size, grid_size*grid_size) + String multiscaleGrid = Geoindicators.GridIndicators.multiscaleLCZGrid(h2gis, gridTable, "id_grid", 1) + String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2gis, multiscaleGrid,grid_size) def folder_save =location in Collection ? location.join("_") : location def path = directory + File.separator + "osm_$folder_save" + File.separator path = "/tmp/" From 7d1b7a68b8c3976881834964fefb7f283c303513 Mon Sep 17 00:00:00 2001 From: ebocher Date: Tue, 19 Mar 2024 14:15:15 +0100 Subject: [PATCH 2/3] Improve tsu, remove small gaps --- .../bdtopo/AbstractBDTopoWorkflow.groovy | 2 +- .../bdtopo/WorkflowAbstractTest.groovy | 4 +-- .../geoindicators/GenericIndicators.groovy | 2 +- .../geoindicators/SpatialUnits.groovy | 4 +-- .../geoindicators/DataUtilsTests.groovy | 2 +- .../GeoIndicatorsExtensionModuleTests.groovy | 4 +-- .../geoindicators/GridIndicatorsTests.groovy | 2 +- .../PopulationIndicatorsTests.groovy | 2 +- .../geoindicators/SpatialUnitsTests.groovy | 17 +++++------- .../WorkflowGeoIndicatorsTest.groovy | 13 +++++----- .../geoclimate/osm/InputDataFormatting.groovy | 1 - .../geoclimate/osm/WorkflowOSM.groovy | 2 +- .../osm/InputDataFormattingTest.groovy | 26 +++++++++---------- .../geoclimate/osm/WorflowOSMTest.groovy | 8 +++--- .../osm/WorkflowAbstractTest.groovy | 6 ++--- .../geoclimate/osmtools/TransformTest.groovy | 4 +-- 16 files changed, 48 insertions(+), 51 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 73b03945d0..7bb16c6efc 100644 --- a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy +++ b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy @@ -1030,7 +1030,7 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { outputFiles.each { if (it == "grid_indicators") { if (outputGrid == "geojson") { - Geoindicators.WorkflowUtilities.saveToGeojson(results."$it", "${outputFolder + File.separator + it}.geojson", h2gis_datasource, outputSRID, reproject, deleteOutputData) + Geoindicators.WorkflowUtilities.saveToGeojson(results."$it", "${outputFolder + File.separator + it}.fgb", h2gis_datasource, outputSRID, reproject, deleteOutputData) } else if (outputGrid == "asc") { Geoindicators.WorkflowUtilities.saveToAscGrid(results."$it", outputFolder, it, h2gis_datasource, outputSRID, reproject, deleteOutputData) } diff --git a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowAbstractTest.groovy b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowAbstractTest.groovy index 569bcd8ec3..5e76443916 100644 --- a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowAbstractTest.groovy +++ b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/WorkflowAbstractTest.groovy @@ -66,7 +66,7 @@ abstract class WorkflowAbstractTest { if(dataFolder){ def files = [:] new File(dataFolder).eachFileRecurse groovy.io.FileType.FILES, { file -> - if (file.name.toLowerCase().endsWith(".geojson")) { + if (file.name.toLowerCase().endsWith(".fgb")) { files.put((file.name.take(file.name.lastIndexOf('.'))), file.getAbsolutePath()) } } @@ -210,7 +210,7 @@ abstract class WorkflowAbstractTest { H2GIS h2gis = H2GIS.open("${folder.absolutePath + File.separator}testFullWorflowSRID;AUTO_SERVER=TRUE") assertTrue h2gis.firstRow("select count(*) as count from $grid_table".toString()).count == 100 assertTrue h2gis.firstRow("select count(*) as count from $grid_table where BUILDING_FRACTION>0".toString()).count > 0 - File grid_file = new File(folder.absolutePath + File.separator + "bdtopo_" + getVersion() + "_" + location.join("_") + File.separator + "grid_indicators.geojson") + File grid_file = new File(folder.absolutePath + File.separator + "bdtopo_" + getVersion() + "_" + location.join("_") + File.separator + "grid_indicators.fgb") assertTrue(grid_file.exists()) h2gis.load(grid_file.absolutePath, "grid_indicators_file", true) assertEquals(4326, h2gis.getSpatialTable("grid_indicators_file").srid) 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 88bfed584b..c689bed99b 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/GenericIndicators.groovy @@ -108,7 +108,7 @@ String unweightedOperationFromLowerScale(JdbcDataSource datasource, String input } } } else { - warn """ The column $var doesn't exist or should be numeric""" + debug("""The column $var doesn't exist or should be numeric""") } } query += "b.$inputIdUp FROM $inputLowerScaleTableName a RIGHT JOIN $inputUpperScaleTableName b " + 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 ddea9619e8..55334b6f6c 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -125,7 +125,7 @@ String createTSU(JdbcDataSource datasource, String inputTableName, String inputz if (inputzone) { String polygons = postfix("polygons") datasource.execute("""DROP TABLE IF EXISTS $polygons; - CREATE TABLE $polygons as SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(ST_BUFFER(the_geom,0.01), $epsg) AS the_geom FROM + CREATE TABLE $polygons as SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(ST_BUFFER(ST_BUFFER(the_geom,-0.01, 'join=mitre'), 0.01, 'join=mitre'), $epsg) AS the_geom FROM ST_EXPLODE('(SELECT ST_POLYGONIZE(ST_UNION(ST_NODE(ST_ACCUM(the_geom)))) AS the_geom FROM $inputTableName)') """.toString()) @@ -145,7 +145,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_buffer(the_geom, -0.01), $epsg) AS the_geom + SELECT EXPLOD_ID AS $COLUMN_ID_NAME, ST_SETSRID(ST_BUFFER(ST_BUFFER(the_geom,-0.01, 'join=mitre'), 0.01, 'join=mitre'), $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() 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 c637bd8ffe..2c41f9d229 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/DataUtilsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/DataUtilsTests.groovy @@ -89,7 +89,7 @@ class DataUtilsTests { directory) assert p - assert 1 == h2GIS.table(h2GIS.load(directory + File.separator + "tablegeom.geojson", true)).rowCount + assert 1 == h2GIS.table(h2GIS.load(directory + File.separator + "tablegeom.fgb", true)).rowCount assert 1 == h2GIS.table(h2GIS.load(directory + File.separator + "tablea.csv", true)).rowCount } } diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GeoIndicatorsExtensionModuleTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GeoIndicatorsExtensionModuleTests.groovy index f3a8fb9c95..d62ce2816f 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GeoIndicatorsExtensionModuleTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GeoIndicatorsExtensionModuleTests.groovy @@ -45,7 +45,7 @@ class GeoIndicatorsExtensionModuleTests { @Test void saveGeometryToFile() { Geometry geom = new GeometryFactory().createPoint(new Coordinate(10, 10)) - File outputFile = new File(folder, "geometry.geojson") + File outputFile = new File(folder, "geometry.fgb") geom.save(h2GIS, outputFile.getAbsolutePath()) h2GIS.load(outputFile.getAbsolutePath(), "mygeom", true) ISpatialTable table = h2GIS.getSpatialTable("mygeom") @@ -58,7 +58,7 @@ class GeoIndicatorsExtensionModuleTests { void saveGeometrySridToFile() { Geometry geom = new GeometryFactory().createPoint(new Coordinate(10, 10)) geom.setSRID(4326) - File outputFile = new File(folder, "geometry.geojson") + File outputFile = new File(folder, "geometry.fgb") geom.save(h2GIS, outputFile.getAbsolutePath()) h2GIS.load(outputFile.getAbsolutePath(), "mygeom", true) ISpatialTable table = h2GIS.getSpatialTable("mygeom") diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy index a3ec4cc5c4..52a6484765 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/GridIndicatorsTests.groovy @@ -89,7 +89,7 @@ class GridIndicatorsTests { from $grid_scale group by ID_ROW_LOD_${i}, ID_COL_LOD_${i}; """.toString()) - h2GIS.save(grid_lod, "/tmp/grid_lod_${i}.geojson", true) + h2GIS.save(grid_lod, "/tmp/grid_lod_${i}.fgb", true) h2GIS.dropTable(grid_lod) } } diff --git a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/PopulationIndicatorsTests.groovy b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/PopulationIndicatorsTests.groovy index d5f0926463..c221de3c32 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/PopulationIndicatorsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/PopulationIndicatorsTests.groovy @@ -91,7 +91,7 @@ class PopulationIndicatorsTests { "building", "rsu", "grid") results.each { it -> - h2GIS.save(it.value, "./target/${it.value}.geojson", true) + h2GIS.save(it.value, "./target/${it.value}.fgb", true) } def rows = h2GIS.rows("SELECT id_build, pop from ${results.buildingTable} order by id_build".toString()) 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 2478dd81d3..402a9048a2 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnitsTests.groovy @@ -310,8 +310,6 @@ class SpatialUnitsTests { """.toString()) String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,"grid","id_grid",1) String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 0) - h2GIS.save(sprawl_areas, "/tmp/sprawl_areas.geojson", true) - h2GIS.save(grid_scales, "/tmp/grid_scales.geojson", true) assertEquals(1, h2GIS.firstRow("select count(*) as count from $sprawl_areas".toString()).count) assertEquals(9, h2GIS.firstRow("select st_area(st_accum(the_geom)) as area from $sprawl_areas".toString()).area, 0.0001) } @@ -339,7 +337,6 @@ class SpatialUnitsTests { ((240 150, 320 150, 320 80, 240 80, 240 150)))'::GEOMETRY as the_geom; """.toString()) String inverse = Geoindicators.SpatialUnits.inversePolygonsLayer(h2GIS, "polygons") - h2GIS.save(inverse, "/tmp/inverse.geojson", true) assertEquals(1, h2GIS.firstRow("select count(*) as count from $inverse".toString()).count) assertTrue(expectedGeom.equals(h2GIS.firstRow("select the_geom from $inverse".toString()).the_geom)) } @@ -364,20 +361,20 @@ class SpatialUnitsTests { void sprawlAreasTestIntegration() { //Data for test String path = "/home/ebocher/Autres/data/geoclimate/uhi_lcz/Angers/" - String data = h2GIS.load("${path}grid_indicators.geojson") + String data = h2GIS.load("${path}grid_indicators.fgb") String grid_scales = Geoindicators.GridIndicators.multiscaleLCZGrid(h2GIS,data,"id_grid", 1) String sprawl_areas = Geoindicators.SpatialUnits.computeSprawlAreas(h2GIS, grid_scales, 100) - h2GIS.save(sprawl_areas, "/tmp/sprawl_areas_indic.geojson", true) - h2GIS.save(grid_scales, "/tmp/grid_indicators.geojson", true) + h2GIS.save(sprawl_areas, "/tmp/sprawl_areas_indic.fgb", true) + h2GIS.save(grid_scales, "/tmp/grid_indicators.fgb", true) String distances = Geoindicators.GridIndicators.gridDistances(h2GIS, sprawl_areas, data, "id_grid") - h2GIS.save(distances, "/tmp/distances.geojson", true) + h2GIS.save(distances, "/tmp/distances.fgb", true) //Method to compute the cool areas distances String cool_areas = Geoindicators.SpatialUnits.extractCoolAreas(h2GIS, grid_scales) - h2GIS.save(cool_areas, "/tmp/cool_areas.geojson", true) + h2GIS.save(cool_areas, "/tmp/cool_areas.fgb", true) String inverse_cool_areas = Geoindicators.SpatialUnits.inversePolygonsLayer(h2GIS,cool_areas) - h2GIS.save(inverse_cool_areas, "/tmp/inverse_cool_areas.geojson", true) + h2GIS.save(inverse_cool_areas, "/tmp/inverse_cool_areas.fgb", true) distances = Geoindicators.GridIndicators.gridDistances(h2GIS, inverse_cool_areas, data, "id_grid") - h2GIS.save(distances, "/tmp/cool_inverse_distances.geojson", true) + h2GIS.save(distances, "/tmp/cool_inverse_distances.fgb", true) } } \ No newline at end of file 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 b3b1749238..e7801eed56 100644 --- a/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy +++ b/geoindicators/src/test/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowGeoIndicatorsTest.groovy @@ -133,6 +133,8 @@ class WorkflowGeoIndicatorsTest { inputTableNames.railTable, inputTableNames.vegetationTable, inputTableNames.hydrographicTable, "", "", "", "","", ["indicatorUse": indicatorUse, svfSimplified: false], prefixName) + datasource.save(inputTableNames.vegetationTable, "/tmp/vegetation.fgb", true) + datasource.save(inputTableNames.roadTable, "/tmp/road.fgb", true) assertNotNull(geoIndicatorsCompute_i) checkRSUIndicators(datasource, geoIndicatorsCompute_i.rsu_indicators) assertEquals(listUrbTyp.Bu.sort(), datasource.getTable(geoIndicatorsCompute_i.building_indicators).columns.sort()) @@ -479,7 +481,6 @@ class WorkflowGeoIndicatorsTest { assertNull(rows.ASPECT_RATIO) assertTrue(0.5 -rows.SVF<0.1) assertEquals(1d, rows.TYPE_OFFICE) - } @@ -490,15 +491,15 @@ class WorkflowGeoIndicatorsTest { def checkRSUIndicators(def datasource, def rsuIndicatorsTableName) { //Check road_fraction > 0 def countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE ROAD_FRACTION>0".toString()) - assertEquals(195, countResult.count) + assertEquals(184, countResult.count) //Check building_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE BUILDING_FRACTION>0".toString()) - assertEquals(73, countResult.count) + assertEquals(70, countResult.count) //Check high_vegetation_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE high_vegetation_fraction>0".toString()) - assertEquals(43, countResult.count) + assertEquals(18, countResult.count) //Check high_vegetation_water_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE high_vegetation_water_fraction>0".toString()) @@ -514,7 +515,7 @@ class WorkflowGeoIndicatorsTest { //Check high_vegetation_road_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE high_vegetation_road_fraction>0".toString()) - assertEquals(43, countResult.count) + assertEquals(20, countResult.count) //Check high_vegetation_impervious_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE high_vegetation_impervious_fraction>0".toString()) @@ -526,7 +527,7 @@ class WorkflowGeoIndicatorsTest { //Check low_vegetation_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE low_vegetation_fraction>0".toString()) - assertEquals(76, countResult.count) + assertEquals(50, countResult.count) //Check low_vegetation_fraction > 0 countResult = datasource.firstRow("select count(*) as count from ${rsuIndicatorsTableName} WHERE impervious_fraction>0".toString()) diff --git a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy index 78dbf968ce..2f20b9e90a 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy @@ -49,7 +49,6 @@ Map formatBuildingLayer(JdbcDataSource datasource, String building, String zone if (!h_lev_min) { h_lev_min = 3 } - def outputTableName = postfix "INPUT_BUILDING" debug 'Formating building layer' def outputEstimateTableName = "EST_${outputTableName}" 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 bde2048844..042b10b11f 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy @@ -364,7 +364,7 @@ Map workflow(def input) { outputDatasource, outputTables, outputSRID, downloadAllOSMData, deleteOutputData, deleteOSMFile, logTableZones, osm_size_area, overpass_timeout, overpass_maxsize, osm_date) if (!osmprocessing) { - h2gis_datasource.save(logTableZones,"${file_outputFolder.getAbsolutePath() + File.separator}logzones.geojson", true) + h2gis_datasource.save(logTableZones,"${file_outputFolder.getAbsolutePath() + File.separator}logzones.fgb", true) return null } if (delete_h2gis) { diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy index 4745e656b3..3e67d9631b 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy @@ -218,7 +218,7 @@ class InputDataFormattingTest { String inputSeaLandTableName = OSM.InputDataFormatting.formatSeaLandMask(h2GIS, extractData.coastline, zoneEnvelopeTableName) assertEquals(2, h2GIS.getTable(inputSeaLandTableName).getRowCount()) assertTrue h2GIS.firstRow("select count(*) as count from ${inputSeaLandTableName} where type='land'").count == 1 - h2GIS.getTable(inputSeaLandTableName).save(new File(folder, "osm_sea_land.geojson").getAbsolutePath(), true) + h2GIS.getTable(inputSeaLandTableName).save(new File(folder, "osm_sea_land.fgb").getAbsolutePath(), true) } @Test @@ -300,42 +300,42 @@ class InputDataFormattingTest { if (extractData.zone != null) { //Zone def epsg = h2GIS.getSpatialTable(extractData.zone).srid - h2GIS.getTable(extractData.zone).save("${file.absolutePath + File.separator}osm_zone_${formatedPlaceName}.geojson", true) + h2GIS.getTable(extractData.zone).save("${file.absolutePath + File.separator}osm_zone_${formatedPlaceName}.fgb", true) //Zone envelope - h2GIS.getTable(extractData.zone_envelope).save("${file.absolutePath + File.separator}osm_zone_envelope_${formatedPlaceName}.geojson", true) + h2GIS.getTable(extractData.zone_envelope).save("${file.absolutePath + File.separator}osm_zone_envelope_${formatedPlaceName}.fgb", true) //Urban Areas def inputUrbanAreas = OSM.InputDataFormatting.formatUrbanAreas(h2GIS, extractData.urban_areas,extractData.zone) - h2GIS.save(inputUrbanAreas,"${file.absolutePath + File.separator}osm_urban_areas_${formatedPlaceName}.geojson", true) + h2GIS.save(inputUrbanAreas,"${file.absolutePath + File.separator}osm_urban_areas_${formatedPlaceName}.fgb", true) println("Urban areas formatted") //Buildings - h2GIS.save(extractData.building,"${file.absolutePath + File.separator}building_${formatedPlaceName}.geojson", true) + h2GIS.save(extractData.building,"${file.absolutePath + File.separator}building_${formatedPlaceName}.fgb", true) def inputBuildings = OSM.InputDataFormatting.formatBuildingLayer(h2GIS, extractData.building,extractData.zone,null) - h2GIS.save(inputBuildings.building,"${file.absolutePath + File.separator}osm_building_${formatedPlaceName}.geojson", true) + h2GIS.save(inputBuildings.building,"${file.absolutePath + File.separator}osm_building_${formatedPlaceName}.fgb", true) println("Building formatted") //Roads def inputRoadTableName = OSM.InputDataFormatting.formatRoadLayer( h2GIS,extractData.road, extractData.zone_envelope) - h2GIS.save(inputRoadTableName,"${file.absolutePath + File.separator}osm_road_${formatedPlaceName}.geojson", true) + h2GIS.save(inputRoadTableName,"${file.absolutePath + File.separator}osm_road_${formatedPlaceName}.fgb", true) println("Road formatted") //Rails def inputRailTableName = OSM.InputDataFormatting.formatRailsLayer( h2GIS,extractData.rail, extractData.zone_envelope) - h2GIS.save(inputRailTableName,"${file.absolutePath + File.separator}osm_rail_${formatedPlaceName}.geojson", true) + h2GIS.save(inputRailTableName,"${file.absolutePath + File.separator}osm_rail_${formatedPlaceName}.fgb", true) println("Rail formatted") //Vegetation def inputVegetationTableName = OSM.InputDataFormatting.formatVegetationLayer( h2GIS,extractData.vegetation,extractData.zone_envelope) - h2GIS.save(inputVegetationTableName,"${file.absolutePath + File.separator}osm_vegetation_${formatedPlaceName}.geojson", true) + h2GIS.save(inputVegetationTableName,"${file.absolutePath + File.separator}osm_vegetation_${formatedPlaceName}.fgb", true) println("Vegetation formatted") @@ -345,23 +345,23 @@ class InputDataFormattingTest { //Impervious String imperviousTable = OSM.InputDataFormatting.formatImperviousLayer(h2GIS, extractData.impervious, extractData.zone_envelope) - h2GIS.save(imperviousTable,"${file.absolutePath + File.separator}osm_impervious_${formatedPlaceName}.geojson", true) + h2GIS.save(imperviousTable,"${file.absolutePath + File.separator}osm_impervious_${formatedPlaceName}.fgb", true) println("Impervious formatted") //Save coastlines to debug - h2GIS.save(extractData.coastline,"${file.absolutePath + File.separator}osm_coastlines_${formatedPlaceName}.geojson", true) + h2GIS.save(extractData.coastline,"${file.absolutePath + File.separator}osm_coastlines_${formatedPlaceName}.fgb", true) //Sea/Land mask def inputSeaLandTableName = OSM.InputDataFormatting.formatSeaLandMask(h2GIS, extractData.coastline, extractData.zone_envelope, inputWaterTableName) - h2GIS.save(inputSeaLandTableName,"${file.absolutePath + File.separator}osm_sea_land_${formatedPlaceName}.geojson", true) + h2GIS.save(inputSeaLandTableName,"${file.absolutePath + File.separator}osm_sea_land_${formatedPlaceName}.fgb", true) println("Sea land mask formatted") //Save it after sea/land mask because the water table can be modified - h2GIS.save(inputWaterTableName,"${file.absolutePath + File.separator}osm_water_${formatedPlaceName}.geojson", true) + h2GIS.save(inputWaterTableName,"${file.absolutePath + File.separator}osm_water_${formatedPlaceName}.fgb", true) } else { assertTrue(false) 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 ae3f365de0..606d79416e 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -336,7 +336,7 @@ class WorflowOSMTest extends WorkflowAbstractTest { def geoFiles = [] def folder = new File("${directory + File.separator}osm_Pont-de-Veyle") folder.eachFileRecurse groovy.io.FileType.FILES, { file -> - if (file.name.toLowerCase().endsWith(".geojson")) { + if (file.name.toLowerCase().endsWith(".fgb")) { geoFiles << file.getAbsolutePath() } } @@ -370,7 +370,7 @@ class WorflowOSMTest extends WorkflowAbstractTest { def folder = new File(directory + File.separator + "osm_" + bbox.join("_")) def countFiles = 0; folder.eachFileRecurse groovy.io.FileType.FILES, { file -> - if (file.name.toLowerCase().endsWith(".geojson")) { + if (file.name.toLowerCase().endsWith(".fgb")) { countFiles++ } } @@ -402,12 +402,12 @@ class WorflowOSMTest extends WorkflowAbstractTest { def folder = new File(directory + File.separator + "osm_" + bbox.join("_")) def resultFiles = [] folder.eachFileRecurse groovy.io.FileType.FILES, { file -> - if (file.name.toLowerCase().endsWith(".geojson")) { + if (file.name.toLowerCase().endsWith(".fgb")) { resultFiles << file.getAbsolutePath() } } assertTrue(resultFiles.size() == 1) - assertTrue(resultFiles.get(0) == folder.absolutePath + File.separator + "zone.geojson") + assertTrue(resultFiles.get(0) == folder.absolutePath + File.separator + "zone.fgb") } @Test diff --git a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorkflowAbstractTest.groovy b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorkflowAbstractTest.groovy index 11bb31ee3b..2aa847762b 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorkflowAbstractTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorkflowAbstractTest.groovy @@ -84,7 +84,7 @@ class WorkflowAbstractTest { } if (saveResults) { logger.debug("Saving building indicators") - datasource.getSpatialTable(buildingIndicators).save(directory + File.separator + "${buildingIndicators}.geojson", true) + datasource.getSpatialTable(buildingIndicators).save(directory + File.separator + "${buildingIndicators}.fgb", true) } //Check we have the same number of buildings @@ -101,7 +101,7 @@ class WorkflowAbstractTest { assertNotNull(blockIndicators) if (saveResults) { logger.debug("Saving block indicators") - datasource.getSpatialTable(blockIndicators).save(directory + File.separator + "${blockIndicators}.geojson", true) + datasource.getSpatialTable(blockIndicators).save(directory + File.separator + "${blockIndicators}.fgb", true) } //Check if we have the same number of blocks def countRelationBlocks = datasource.firstRow("select count(*) as count from ${relationBlocks}".toString()) @@ -125,7 +125,7 @@ class WorkflowAbstractTest { assert rsuIndicators if (saveResults) { logger.debug("Saving RSU indicators") - datasource.getSpatialTable(rsuIndicators).save(directory + File.separator + "${rsuIndicators}.geojson", true) + datasource.getSpatialTable(rsuIndicators).save(directory + File.separator + "${rsuIndicators}.fgb", true) } //Check if we have the same number of RSU diff --git a/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy b/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy index 55ef11a540..14e9ae1943 100644 --- a/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy +++ b/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy @@ -848,7 +848,7 @@ class TransformTest extends AbstractOSMToolsTest { //Create building layer def tags = ["amenity", "landuse", "railway", "water"] String outputTableName = OSMTools.Transform.toPolygons(ds, prefix, 4326, tags) - ds.save(outputTableName, "/tmp/polygons.geojson", true) + ds.save(outputTableName, "/tmp/polygons.fgb", true) } } } @@ -863,7 +863,7 @@ class TransformTest extends AbstractOSMToolsTest { Map r = OSMTools.Loader.fromArea(h2GIS, [48.733493,-3.076869,48.733995,-3.075829]) println(r) def lines = OSMTools.Transform.toPolygons(h2GIS, r.prefix,4326, [], []) - h2GIS.save(lines, "/tmp/building.geojson") + h2GIS.save(lines, "/tmp/building.fgb") } From 8f8651b4bcb3d258bb56351531a26c7e757ee979 Mon Sep 17 00:00:00 2001 From: ebocher Date: Wed, 20 Mar 2024 10:12:43 +0100 Subject: [PATCH 3/3] Fix test --- .../bdtopo/AbstractBDTopoWorkflow.groovy | 48 +++++++++---------- .../bdtopo/WorkflowAbstractTest.groovy | 3 +- .../bdtopo/v2/WorkflowBDTopoV2Test.groovy | 27 ++--------- .../bdtopo/v3/WorkflowBDTopoV3Test.groovy | 10 ++-- .../geoindicators/WorkflowUtilities.groovy | 4 +- .../geoclimate/osm/WorkflowOSM.groovy | 16 +++---- 6 files changed, 44 insertions(+), 64 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 7bb16c6efc..cf36e6b7aa 100644 --- a/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy +++ b/bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/AbstractBDTopoWorkflow.groovy @@ -447,20 +447,20 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { def defaultParameters = [distance : 1000f, distance_buffer: 500f, prefixName: "", hLevMin : 3] - def rsu_indicators_default = [indicatorUse : [], - svfSimplified : true, - surface_vegetation: 10000f, - surface_hydro : 2500f, + def rsu_indicators_default = [indicatorUse : [], + svfSimplified : true, + surface_vegetation : 10000f, + surface_hydro : 2500f, surface_urban_areas: 10000f, - snappingTolerance : 0.01f, - mapOfWeights : ["sky_view_factor" : 4, - "aspect_ratio" : 3, - "building_surface_fraction" : 8, - "impervious_surface_fraction" : 0, - "pervious_surface_fraction" : 0, - "height_of_roughness_elements": 6, - "terrain_roughness_length" : 0.5], - utrfModelName : "UTRF_BDTOPO_V2_RF_2_2.model"] + snappingTolerance : 0.01f, + mapOfWeights : ["sky_view_factor" : 4, + "aspect_ratio" : 3, + "building_surface_fraction" : 8, + "impervious_surface_fraction" : 0, + "pervious_surface_fraction" : 0, + "height_of_roughness_elements": 6, + "terrain_roughness_length" : 0.5], + utrfModelName : "UTRF_BDTOPO_V2_RF_2_2.model"] defaultParameters.put("rsu_indicators", rsu_indicators_default) if (processing_parameters) { @@ -564,13 +564,13 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { def grid_indicators_tmp = [ "x_size" : x_size, "y_size" : y_size, - "output" : "geojson", + "output" : "fgb", "rowCol" : false, "indicators": allowedOutputIndicators ] def grid_output = grid_indicators.output if (grid_output) { - if (grid_output.toLowerCase() in ["asc", "geojson"]) { + if (grid_output.toLowerCase() in ["asc", "fgb"]) { grid_indicators_tmp.output = grid_output.toLowerCase() } } @@ -675,7 +675,7 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { //Add the GIS layers to the list of results def outputTableNamesResult = [:] def grid_indicators_params = processing_parameters.grid_indicators - def outputGrid = "geojson" + def outputGrid = "fgb" if (grid_indicators_params) { outputGrid = grid_indicators_params.output } @@ -936,7 +936,7 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { Map geoIndicators = Geoindicators.WorkflowGeoIndicators.computeAllGeoIndicators(h2gis_datasource, zone, building, road, rail, vegetation, - water, impervious, "", "", urban_areas,"", + water, impervious, "", "", urban_areas, "", rsu_indicators_params, processing_parameters.prefixName) if (!geoIndicators) { error "Cannot build the geoindicators for the zone $id_zone" @@ -961,7 +961,7 @@ abstract class AbstractBDTopoWorkflow extends BDTopoUtils { if (!building) { info "Cannot compute any population data at building level" } - tablesToDrop<