Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix issue 985 #986

Merged
merged 5 commits into from
Jul 3, 2024
Merged

Fix issue 985 #986

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
3845d38
Fix issue 985
ebocher Jul 3, 2024
6414968
Merge branch 'master' into issue_985
ebocher Jul 3, 2024
42b0c1c
Fix test
ebocher Jul 3, 2024
d1a2aec
Merge remote-tracking branch 'ebocher/issue_985' into issue_985
ebocher Jul 3, 2024
70f502e
Fix test
ebocher Jul 3, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 13 additions & 7 deletions osm/src/main/groovy/org/orbisgis/geoclimate/osm/InputDataFormatting.groovy
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading