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 e5e44ea9b3..305c91c41f 100644 --- a/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy +++ b/osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy @@ -1148,31 +1148,37 @@ String formatSeaLandMask(JdbcDataSource datasource, String coastline, String zon SELECT the_geom, EXPLOD_ID FROM st_explode('( SELECT ST_LINEMERGE(st_accum(THE_GEOM)) AS the_geom, NULL FROM $coastLinesIntersects)');""".toString() - datasource.execute """ + //Reduce the zone to limit intersection rounding with JTS + String reducedZone = postfix("reduced_zone") + datasource.execute("""DROP TABLE IF EXISTS $reducedZone; + CREATE TABLE $reducedZone as select st_buffer(the_geom, -0.01,2) as the_geom from $zone;""") + + datasource.execute(""" CREATE TABLE $mergingDataTable AS SELECT THE_GEOM FROM $coastLinesIntersects UNION ALL - SELECT st_tomultiline(st_buffer(the_geom, -0.01,2)) - from $zone + SELECT the_geom + from $reducedZone UNION ALL SELECT st_tomultiline(the_geom) from $water ; CREATE TABLE $sea_land_mask (THE_GEOM GEOMETRY,ID serial, TYPE VARCHAR, ZINDEX INTEGER) AS SELECT THE_GEOM, EXPLOD_ID, 'land', 0 AS ZINDEX FROM st_explode('(SELECT st_polygonize(st_union(ST_NODE(st_accum(the_geom)))) AS the_geom FROM $mergingDataTable)') - as foo where ST_DIMENSION(the_geom) = 2 AND st_area(the_geom) >0; """.toString() + as foo where ST_DIMENSION(the_geom) = 2 AND st_area(the_geom) >0; """) - datasource.execute """ + datasource.execute(""" CREATE SPATIAL INDEX IF NOT EXISTS ${sea_land_mask}_the_geom_idx ON $sea_land_mask (THE_GEOM); CREATE SPATIAL INDEX IF NOT EXISTS ${islands_mark}_the_geom_idx ON $islands_mark (THE_GEOM); CREATE TABLE $coastLinesPoints as SELECT ST_LocateAlong(the_geom, 0.5, -0.01) AS the_geom FROM st_explode('(select ST_GeometryN(ST_ToMultiSegments(st_intersection(a.the_geom, b.the_geom)), 1) as the_geom from $islands_mark as a, - $zone as b WHERE a.the_geom && b.the_geom AND st_intersects(a.the_geom, b.the_geom))'); + $reducedZone as b WHERE a.the_geom && b.the_geom AND st_intersects(a.the_geom, b.the_geom))'); CREATE TABLE $coastLinesIntersectsPoints as SELECT the_geom FROM st_explode('$coastLinesPoints'); - CREATE SPATIAL INDEX IF NOT EXISTS ${coastLinesIntersectsPoints}_the_geom_idx ON $coastLinesIntersectsPoints (THE_GEOM);""".toString() + CREATE SPATIAL INDEX IF NOT EXISTS ${coastLinesIntersectsPoints}_the_geom_idx ON $coastLinesIntersectsPoints (THE_GEOM); + DROP TABLE IF EXISTS $reducedZone;""") //Perform triangulation to tag the areas as sea or water datasource.execute """ 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 27ec7fdfbf..d3bbad4f09 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy @@ -286,16 +286,16 @@ class InputDataFormattingTest { if (!file.exists()) { file.mkdir() } - def h2GIS = H2GIS.open("${file.absolutePath + File.separator}osm_gislayers;AUTO_SERVER=TRUE".toString()) - def zoneToExtract = "Lorient" + + def zoneToExtract = "Vannes" //def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(zoneToExtract) // zoneToExtract = nominatim.bbox //zoneToExtract = [62.2, 28.2, 62.4, 28.4] - //zoneToExtract =[44.795480,12.323227,45.004622,12.627411] + //zoneToExtract =[45.784554,4.861279,45.796015,4.883981] Map extractData = OSM.InputDataLoading.extractAndCreateGISLayers(h2GIS, zoneToExtract) String formatedPlaceName = zoneToExtract.join("-").trim().split("\\s*(,|\\s)\\s*").join("_"); @@ -346,10 +346,11 @@ class InputDataFormattingTest { println("Vegetation formatted") //Hydrography + h2GIS.save(extractData.water,"${file.absolutePath + File.separator}water_osm.fgb", true) def inputWaterTableName = OSM.InputDataFormatting.formatWaterLayer(h2GIS, extractData.water, extractData.zone_envelope) //Impervious - String imperviousTable = OSM.InputDataFormatting.formatImperviousLayer(h2GIS, extractData.impervious, + String imperviousTable = OSM.InputDataFormatting.formatImperviousLayer(h2GIS, extractData.impervious, extractData.zone_envelope) h2GIS.save(imperviousTable,"${file.absolutePath + File.separator}impervious.fgb", true) @@ -362,7 +363,6 @@ class InputDataFormattingTest { //Sea/Land mask def inputSeaLandTableName = OSM.InputDataFormatting.formatSeaLandMask(h2GIS, extractData.coastline, extractData.zone_envelope, inputWaterTableName) - println(inputSeaLandTableName.isEmpty()) h2GIS.save(inputSeaLandTableName, "${file.absolutePath + File.separator}sea_land_mask.fgb", true) println("Sea land mask formatted") 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 464199f1f1..4a720ac424 100644 --- a/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy +++ b/osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy @@ -649,12 +649,15 @@ class WorflowOSMTest extends WorkflowAbstractTest { File dirFile = new File(directory) dirFile.delete() dirFile.mkdir() - def location = "Redon" - //def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location) + def location = "Villeurbanne" + def nominatim = org.orbisgis.geoclimate.osmtools.OSMTools.Utilities.getNominatimData(location) def grid_size = 250 - //location = nominatim.bbox - //location=[46.178404,6.095524,46.222959,6.190109] - //location =[47.693558, -3.5087318, 47.814945, -3.284718] + location = nominatim.bbox + //location=[46.178404,6.095524,46.222959,6.190109] + /* location =[ 48.84017284026897, + 2.3061887733275785, + 48.878115442982086, + 2.36742047202511]*/ def osm_parmeters = [ "description" : "Example of configuration file to run the OSM workflow and store the result in a folder", "geoclimatedb": [ @@ -677,25 +680,25 @@ class WorflowOSMTest extends WorkflowAbstractTest { "rsu_indicators" : [ "indicatorUse": ["LCZ"]//, "TEB"] //, "UTRF"] - ], /*"grid_indicators" : [ + ], "grid_indicators" : [ "x_size" : grid_size, "y_size" : grid_size, - //"rowCol": true, - "indicators": [//"BUILDING_FRACTION","BUILDING_HEIGHT", "BUILDING_POP", - //"BUILDING_TYPE_FRACTION", - //"WATER_FRACTION","VEGETATION_FRACTION", - //"ROAD_FRACTION", "IMPERVIOUS_FRACTION", - "LCZ_PRIMARY", - //"BUILDING_HEIGHT_WEIGHTED", "BUILDING_SURFACE_DENSITY", - //"SEA_LAND_FRACTION", - //"ASPECT_RATIO", - //"SVF", - // "HEIGHT_OF_ROUGHNESS_ELEMENTS", "TERRAIN_ROUGHNESS_CLASS", - "URBAN_SPRAWL_AREAS", - "URBAN_SPRAWL_DISTANCES", - "URBAN_SPRAWL_COOL_DISTANCES"]*//*, - "lcz_lod":1*//* - ], "worldpop_indicators": true*/ + "indicators": [ + "BUILDING_FRACTION", + "BUILDING_HEIGHT", + "BUILDING_POP", + "BUILDING_TYPE_FRACTION", + "WATER_FRACTION", + "VEGETATION_FRACTION", + "ROAD_FRACTION", + "IMPERVIOUS_FRACTION", + "LCZ_FRACTION", + "LCZ_PRIMARY", + "SEA_LAND_FRACTION", + "ASPECT_RATIO", + "SVF" ] + //"lcz_lod":1 + ], "worldpop_indicators": true /* "road_traffic" : true, diff --git a/osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy b/osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy index 28bc44a46f..57611dc3e1 100644 --- a/osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy +++ b/osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy @@ -298,7 +298,6 @@ String extractWaysAsPolygons(JdbcDataSource datasource, String osmTablesPrefix, return outputTableName } } - datasource """ CREATE TABLE $waysPolygonTmp AS SELECT ST_TRANSFORM(ST_SETSRID(ST_MAKEPOLYGON(ST_MAKELINE(the_geom)), 4326), $epsgCode) AS the_geom, id_way @@ -415,13 +414,10 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix error "The tag list cannot be null" return } - + def tablesToDrop =[] def outputTableName = postfix "RELATION_POLYGONS" def osmTableTag = prefix osmTablesPrefix, "relation_tag" - def countTagsQuery = OSMTools.TransformUtils.getCountTagsQuery(osmTableTag, tags) - def columnsSelector = OSMTools.TransformUtils.getColumnSelector(osmTableTag, tags, columnsToKeep) def tagsFilter = OSMTools.TransformUtils.createWhereFilter(tags) - if (!datasource.hasTable(osmTableTag)) { debug "No tags table found to build the table. An empty table will be returned." datasource """ @@ -431,59 +427,56 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix return outputTableName } - if (datasource.firstRow(countTagsQuery.toString()).count <= 0) { + def relationFilteredKeys = postfix "RELATION_FILTERED_KEYS" + tablesToDrop<=2) GROUP BY id_relation; """.toString()) @@ -504,8 +497,8 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix datasource """ DROP TABLE IF EXISTS $relationsPolygonsInner; CREATE TABLE $relationsPolygonsInner AS - SELECT ST_LINEMERGE(ST_ACCUM(the_geom)) the_geom, id_relation - FROM( + SELECT st_linemerge(ST_ACCUM(the_geom)) the_geom, id_relation + FROM ( SELECT ST_TRANSFORM(ST_SETSRID(ST_MAKELINE(the_geom), 4326), $epsgCode) AS the_geom, id_relation, role, id_way FROM( SELECT( @@ -514,8 +507,8 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix SELECT n.id_node, n.the_geom, wn.id_way idway FROM ${osmTablesPrefix}_node n, ${osmTablesPrefix}_way_node wn WHERE n.id_node = wn.id_node ORDER BY wn.node_order) - WHERE idway = w.id_way) the_geom, w.id_way, br.id_relation, br.role - FROM ${osmTablesPrefix}_way w, ${osmTablesPrefix}_way_member br ${inner_condition}) geom_table + WHERE idway = br.id_way) the_geom, br.id_way, br.id_relation, br.role + FROM ${ways_list_relations} as br WHERE br.role='inner') geom_table WHERE st_numgeometries(the_geom)>=2) GROUP BY id_relation; """.toString() @@ -559,6 +552,8 @@ def extractRelationsAsPolygons(JdbcDataSource datasource, String osmTablesPrefix CREATE INDEX ON $relationsMpHoles(id_relation); """.toString() + //select the column to keep + def columnsSelector = OSMTools.TransformUtils.getColumnSelector(osmTableTag, tags, columnsToKeep) String allRelationPolygons = postfix("all_relation_polygons") def query = """ DROP TABLE IF EXISTS $allRelationPolygons; 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 9b230b37ea..e0239daeeb 100644 --- a/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy +++ b/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/TransformTest.groovy @@ -208,7 +208,7 @@ class TransformTest extends AbstractOSMToolsTest { String result = OSMTools.Transform.toPolygons(ds, prefix, epsgCode, tags, columnsToKeep) assertFalse result.isEmpty() def table = ds.getSpatialTable(result) - assertEquals 2, table.rowCount + assertEquals 1, table.rowCount table.each { switch (it.row) { case 1: @@ -216,11 +216,6 @@ class TransformTest extends AbstractOSMToolsTest { assertEquals 2, it.the_geom.getDimension() assertEquals "lake", it.water break - case 2: - assertEquals "r1", it.id - assertEquals 2, it.the_geom.getDimension() - assertEquals "lake", it.water - break } } @@ -369,18 +364,9 @@ class TransformTest extends AbstractOSMToolsTest { String result = OSMTools.Transform.extractRelationsAsPolygons(ds, prefix, epsgCode, tags, columnsToKeep) assertFalse result.isEmpty() + ds.save(result, "/tmp/building.geojson", true) def table = ds.getTable(result) - assertEquals 1, table.rowCount - table.each { - switch (it.row) { - case 1: - assertEquals "r1", it.id - assertEquals 2, it.the_geom.getDimension() - assertEquals "house", it.building - assertEquals "lake", it.water - break - } - } + assertEquals 0, table.rowCount //Test not existing tags result = OSMTools.Transform.extractRelationsAsPolygons(ds, prefix, epsgCode, [toto: "tata"], []) @@ -391,17 +377,7 @@ class TransformTest extends AbstractOSMToolsTest { CREATE TABLE ${prefix}_node_tag (id_node int, tag_key varchar, tag_value varchar);""".toString() result = OSMTools.Transform.extractRelationsAsPolygons(ds, prefix, epsgCode, [], []) table = ds.getTable(result) - assertEquals 1, table.rowCount - table.each { - switch (it.row) { - case 1: - assertEquals "r1", it.id - assertEquals 2, it.the_geom.getDimension() - assertEquals "house", it.building - assertEquals "lake", it.water - break - } - } + assertEquals 0, table.rowCount //Test column to keep absent result = OSMTools.Transform.extractRelationsAsPolygons(ds, prefix, epsgCode, tags, ["landscape"]) diff --git a/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/utils/TransformUtilsTest.groovy b/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/utils/TransformUtilsTest.groovy index b440b5d339..92c802c0e3 100644 --- a/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/utils/TransformUtilsTest.groovy +++ b/osmtools/src/test/groovy/org/orbisgis/geoclimate/osmtools/utils/TransformUtilsTest.groovy @@ -602,7 +602,7 @@ ${osmTablesPrefix}_way_member, ${osmTablesPrefix}_way_not_taken_into_account, ${ result = OSMTools.TransformUtils.toPolygonOrLine(polygonType, h2gis, prefix, epsgCode, tags, columnsToKeep) assertNotNull result table = h2gis.getTable(result) - assertEquals 2, table.rowCount + assertEquals 1, table.rowCount table.each { switch (it.row) { case 1: @@ -611,32 +611,16 @@ ${osmTablesPrefix}_way_member, ${osmTablesPrefix}_way_not_taken_into_account, ${ assertEquals 2, it.the_geom.getDimension() assertEquals "lake", it.water break - case 2: - assertEquals "house", it.building - assertEquals "r1", it.id - assertEquals 2, it.the_geom.getDimension() - assertEquals "lake", it.water - break } } //Test no way tags - h2gis.execute """DROP TABLE ${prefix}_way_tag; - CREATE TABLE ${prefix}_way_tag (id_way int, tag_key varchar, tag_value varchar);""".toString() + h2gis.execute("""DROP TABLE ${prefix}_way_tag; + CREATE TABLE ${prefix}_way_tag (id_way int, tag_key varchar, tag_value varchar);""") result = OSMTools.TransformUtils.toPolygonOrLine(polygonType, h2gis, prefix, epsgCode, tags, columnsToKeep) assertNotNull result table = h2gis.getTable(result) - assertEquals 1, table.rowCount - table.each { - switch (it.row) { - case 1: - assertEquals "house", it.building - assertEquals "r1", it.id - assertEquals 2, it.the_geom.getDimension() - assertEquals "lake", it.water - break - } - } + assertEquals 0, table.rowCount result = OSMTools.TransformUtils.toPolygonOrLine(lineType, h2gis, prefix, epsgCode, tags, columnsToKeep) assertNotNull result table = h2gis.getTable(result)