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..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< - if (file.name.toLowerCase().endsWith(".geojson")) { + if (file.name.toLowerCase().endsWith(".fgb")) { files.put((file.name.take(file.name.lastIndexOf('.'))), file.getAbsolutePath()) } } @@ -210,7 +209,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/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v2/WorkflowBDTopoV2Test.groovy b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v2/WorkflowBDTopoV2Test.groovy index f37b7630f7..3aec0944b1 100644 --- a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v2/WorkflowBDTopoV2Test.groovy +++ b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v2/WorkflowBDTopoV2Test.groovy @@ -19,36 +19,17 @@ */ package org.orbisgis.geoclimate.bdtopo.v2 +import org.junit.jupiter.api.io.CleanupMode +import org.junit.jupiter.api.io.TempDir import org.orbisgis.data.H2GIS import org.orbisgis.geoclimate.bdtopo.WorkflowAbstractTest - -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertEquals -import static org.junit.jupiter.api.Assertions.assertTrue -import static org.junit.jupiter.api.Assertions.assertTrue -import static org.junit.jupiter.api.Assertions.assertTrue -import static org.junit.jupiter.api.Assertions.assertTrue import static org.junit.jupiter.api.Assertions.assertTrue class WorkflowBDTopoV2Test extends WorkflowAbstractTest { + @TempDir(cleanup = CleanupMode.ON_SUCCESS) + static File folder @Override ArrayList getFileNames() { diff --git a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v3/WorkflowBDTopoV3Test.groovy b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v3/WorkflowBDTopoV3Test.groovy index 2bdc37c10f..a6ee12691c 100644 --- a/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v3/WorkflowBDTopoV3Test.groovy +++ b/bdtopo/src/test/groovy/org/orbisgis/geoclimate/bdtopo/v3/WorkflowBDTopoV3Test.groovy @@ -20,7 +20,8 @@ package org.orbisgis.geoclimate.bdtopo.v3 -import org.junit.jupiter.api.* +import org.junit.jupiter.api.io.CleanupMode +import org.junit.jupiter.api.io.TempDir import org.orbisgis.data.H2GIS import org.orbisgis.geoclimate.bdtopo.WorkflowAbstractTest @@ -29,11 +30,10 @@ import static org.junit.jupiter.api.Assertions.assertTrue class WorkflowBDTopoV3Test extends WorkflowAbstractTest { + @TempDir(cleanup = CleanupMode.ON_SUCCESS) + static File folder + - public WorkflowBDTopoV3Test(){ - new File("/tmp/test_bd").mkdir() - folder = new File("/tmp/test_bd") - } @Override int getVersion() { return 3 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/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 b5e29b27e2..55334b6f6c 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/SpatialUnits.groovy @@ -123,23 +123,29 @@ 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(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()) + 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 """ 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/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowUtilities.groovy b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowUtilities.groovy index f678650bcf..efb4391879 100644 --- a/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowUtilities.groovy +++ b/geoindicators/src/main/groovy/org/orbisgis/geoclimate/geoindicators/WorkflowUtilities.groovy @@ -212,7 +212,7 @@ def saveToAscGrid(def outputTable, def subFolder, def filePrefix, JdbcDataSource } /** - * Method to save a table into a geojson file + * Method to save a table into a file * @param outputTable name of the table to export * @param filePath path to save the table * @param h2gis_datasource connection to the database @@ -220,7 +220,7 @@ def saveToAscGrid(def outputTable, def subFolder, def filePrefix, JdbcDataSource * @param reproject true if the file must be reprojected * @param deleteOutputData true to delete the file if exists */ -def saveToGeojson(def outputTable, def filePath, H2GIS h2gis_datasource, def outputSRID, def reproject, def deleteOutputData) { +def saveInFile(def outputTable, def filePath, H2GIS h2gis_datasource, def outputSRID, def reproject, def deleteOutputData) { if (outputTable && h2gis_datasource.hasTable(outputTable)) { if (!reproject) { h2gis_datasource.save(outputTable, filePath, deleteOutputData) 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 e9e9c62708..2f20b9e90a 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 @@ -47,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}" @@ -72,44 +73,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 +119,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 +262,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 +1073,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..0152ed1343 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy @@ -71,7 +71,7 @@ import java.sql.SQLException * * , * * [OPTIONAL ENTRY] "output" :{ //If not ouput is set the results are keep in the local database * * "srid" : //optional value to reproject the data - * * "folder" : "/tmp/myResultFolder" //tmp folder to store the computed layers in a geojson format, + * * "folder" : "/tmp/myResultFolder" //tmp folder to store the computed layers in a fgb format, * * "database": { //database parameters to store the computed layers. * * "user": "-", * * "password": "-", @@ -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) { @@ -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", @@ -615,10 +633,10 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d } } //Default - def outputGrid = "geojson" + def outputGrid = "fgb" 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" @@ -891,13 +916,13 @@ def extractProcessingParameters(def processing_parameters) { 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() } } @@ -942,7 +967,7 @@ def extractProcessingParameters(def processing_parameters) { /** - * Save the geoclimate tables into geojson files + * Save the geoclimate tables into files * @param id_zone the id of the zone * @param results a list of tables computed by geoclimate * @param ouputFolder the ouput folder @@ -963,15 +988,15 @@ def saveOutputFiles(def h2gis_datasource, def id_zone, def results, def outputFi } outputFiles.each { if (it == "grid_indicators") { - if (outputGrid == "geojson") { - Geoindicators.WorkflowUtilities.saveToGeojson(results."$it", "${subFolder.getAbsolutePath() + File.separator + it}.geojson", h2gis_datasource, outputSRID, reproject, deleteOutputData) + if (outputGrid == "fgb") { + Geoindicators.WorkflowUtilities.saveInFile(results."$it", "${subFolder.getAbsolutePath() + File.separator + it}.fgb", h2gis_datasource, outputSRID, reproject, deleteOutputData) } else if (outputGrid == "asc") { Geoindicators.WorkflowUtilities.saveToAscGrid(results."$it", subFolder.getAbsolutePath(), it, h2gis_datasource, outputSRID, reproject, deleteOutputData) } } else if (it == "building_height_missing") { Geoindicators.WorkflowUtilities.saveToCSV(results."$it", "${subFolder.getAbsolutePath() + File.separator + it}.csv", h2gis_datasource, deleteOutputData) } else { - Geoindicators.WorkflowUtilities.saveToGeojson(results."$it", "${subFolder.getAbsolutePath() + File.separator + it}.geojson", h2gis_datasource, outputSRID, reproject, deleteOutputData) + Geoindicators.WorkflowUtilities.saveInFile(results."$it", "${subFolder.getAbsolutePath() + File.separator + it}.fgb", h2gis_datasource, outputSRID, reproject, deleteOutputData) } } } 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..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 @@ -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) @@ -300,37 +300,44 @@ 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,inputUrbanAreas) - h2GIS.save(inputBuildings.building,"${file.absolutePath + File.separator}osm_building_${formatedPlaceName}.geojson", true) + extractData.building,extractData.zone,null) + 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") //Hydrography def inputWaterTableName = OSM.InputDataFormatting.formatWaterLayer(h2GIS, extractData.water, extractData.zone_envelope) @@ -338,19 +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 4b10af7d9c..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 @@ -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/" 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") }