Skip to content

Commit

Permalink
Merge pull request #986 from ebocher/issue_985
Browse files Browse the repository at this point in the history 
Fix issue 985
  • Loading branch information
@ebocher
ebocher authored Jul 3, 2024
2 parents 54b00f2 + 70f502e commit a323350
Show file tree
Hide file tree
Showing 6 changed files with 97 additions and 133 deletions.
20 changes: 13 additions & 7 deletions osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 """
Expand Down
10 changes: 5 additions & 5 deletions osm/src/test/groovy/org/orbisgis/geoclimate/osm/InputDataFormattingTest.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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("_");
Expand Down Expand Up @@ -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)

Expand All @@ -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")
Expand Down
47 changes: 25 additions & 22 deletions osm/src/test/groovy/org/orbisgis/geoclimate/osm/WorflowOSMTest.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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": [
Expand All @@ -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,
Expand Down
97 changes: 46 additions & 51 deletions osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/Transform.groovy
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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 """
Expand All @@ -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<<relationFilteredKeys
if (tagsFilter) {
datasource.execute("""
DROP TABLE IF EXISTS $relationFilteredKeys;
CREATE TABLE $relationFilteredKeys AS
SELECT DISTINCT a.id_relation
FROM ${osmTablesPrefix}_relation_tag as a
right JOIN
(SELECT ID_RELATION FROM ${osmTablesPrefix}_relation_tag
WHERE tag_key='type' AND tag_value='multipolygon') as b
ON a.ID_RELATION = b.id_relation
WHERE $tagsFilter;
""")
}else {
datasource.execute("""
DROP TABLE IF EXISTS $relationFilteredKeys;
CREATE TABLE $relationFilteredKeys AS
SELECT DISTINCT id_relation
FROM ${osmTablesPrefix}_relation_tag WHERE tag_key='type' AND tag_value='multipolygon';
""")
}

if (datasource.getRowCount(relationFilteredKeys) <= 0) {
debug "No keys or values found in the relations. An empty table will be returned."
datasource """
DROP TABLE IF EXISTS $outputTableName;
CREATE TABLE $outputTableName (the_geom GEOMETRY(GEOMETRY,$epsgCode));
""".toString()
datasource.dropTable(tablesToDrop)
return outputTableName
}
debug "Build outer polygons"
def relationsPolygonsOuter = postfix "RELATIONS_POLYGONS_OUTER"
def relationFilteredKeys = postfix "RELATION_FILTERED_KEYS"
def outer_condition
def inner_condition
if (!tagsFilter) {
outer_condition = "WHERE w.id_way = br.id_way AND br.role='outer'"
inner_condition = "WHERE w.id_way = br.id_way AND br.role='inner'"
} else {
datasource """
DROP TABLE IF EXISTS $relationFilteredKeys;
CREATE TABLE $relationFilteredKeys AS
SELECT DISTINCT id_relation
FROM ${osmTablesPrefix}_relation_tag as a
WHERE $tagsFilter;
CREATE INDEX ON $relationFilteredKeys(id_relation);
""".toString()

if (columnsToKeep) {
if (datasource.getRowCount(relationFilteredKeys) < 1) {
debug "Any columns to keep. Cannot create any geometry polygons. An empty table will be returned."
datasource """
DROP TABLE IF EXISTS $outputTableName;
CREATE TABLE $outputTableName (the_geom GEOMETRY(GEOMETRY,$epsgCode));
""".toString()
return outputTableName
}
}
//Extract the corresponding way by relation
def ways_list_relations = postfix "WAY_LIST_RELATIONS"
tablesToDrop<<ways_list_relations

outer_condition = """, $relationFilteredKeys g
WHERE br.id_relation=g.id_relation
AND w.id_way = br.id_way
AND br.role='outer'
"""
inner_condition = """, $relationFilteredKeys g
WHERE br.id_relation=g.id_relation
AND w.id_way = br.id_way
AND br.role='inner'
"""
}
datasource.execute("""
CREATE INDEX ON $relationFilteredKeys(id_relation);
DROP TABLE IF EXISTS $ways_list_relations;
CREATE TABLE $ways_list_relations as SELECT br.id_way, br.id_relation, br.role, br.way_order
FROM ${osmTablesPrefix}_way_member br, $relationFilteredKeys as r WHERE br.id_relation=r.id_relation""")

debug "Build outer polygons"
def relationsPolygonsOuter = postfix "RELATIONS_POLYGONS_OUTER"
datasource.execute("""
CREATE INDEX ON $relationFilteredKeys(id_relation);
DROP TABLE IF EXISTS $relationsPolygonsOuter;
CREATE TABLE $relationsPolygonsOuter AS
SELECT ST_LINEMERGE(ST_ACCUM(the_geom)) as the_geom, id_relation
SELECT st_linemerge(ST_ACCUM(the_geom)) as 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(
Expand All @@ -493,8 +486,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 $outer_condition) geom_table
WHERE idway = br.id_way) the_geom, br.id_way, br.id_relation, br.role , br.way_order
FROM ${ways_list_relations} as br WHERE br.role='outer' ) geom_table
WHERE st_numgeometries(the_geom)>=2)
GROUP BY id_relation;
""".toString())
Expand All @@ -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(
Expand All @@ -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()
Expand Down Expand Up @@ -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;
Expand Down
Loading

0 comments on commit a323350

Please sign in to comment.