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

Jump to bottom

OSM and BDTopo data extraction options #924

Merged
merged 3 commits into from
Jan 18, 2024
Merged

OSM and BDTopo data extraction options #924

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
987ab52
This PR introduces two new options to process OSM data
Jan 17, 2024
a7e2cad
Fix bbox utility
Jan 17, 2024
c4ec7ca
Add point location option to bdtopo module
Jan 18, 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
17 changes: 16 additions & 1 deletion bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/BDTopoV2Workflow.groovy
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,13 @@ Integer loadDataFromPostGIS(Object input_database_properties, Object code, Objec
//Check if code is a string or a bbox
//The zone is a osm bounding box represented by ymin,xmin , ymax,xmax,
if (code in Collection) {
//def tmp_insee = code.join("_")
if(code.size()==3){
if(code[2]<100){
error("The distance to create a bbox from a point must be greater than 100 meters")
return
}
code = BDTopoUtils.bbox(code[0], code[1],code[2])
}
String inputTableName = """(SELECT
ST_INTERSECTION(st_setsrid(the_geom, $commune_srid), ST_MakeEnvelope(${code[1]},${code[0]},${code[3]},${code[2]}, $commune_srid)) as the_geom, CODE_INSEE from $commune_location where
st_setsrid(the_geom, $commune_srid)
Expand Down Expand Up @@ -306,6 +312,13 @@ def filterLinkedShapeFiles(def location, float distance, LinkedHashMap inputTabl
//Check if code is a string or a bbox
//The zone is a osm bounding box represented by ymin,xmin , ymax,xmax,
if (location in Collection) {
if(location.size()==3){
if(location[2]<100){
error("The distance to create a bbox from a point must be greater than 100 meters")
return
}
location = BDTopoUtils.bbox(location[0], location[1],location[2])
}
debug "Loading in the H2GIS database $outputTableName"
h2gis_datasource.execute("""DROP TABLE IF EXISTS $outputTableName ; CREATE TABLE $outputTableName as SELECT
ST_INTERSECTION(the_geom, ST_MakeEnvelope(${location[1]},${location[0]},${location[3]},${location[2]}, $sourceSRID)) as the_geom, CODE_INSEE from ${inputTables.commune} where the_geom
Expand All @@ -314,6 +327,8 @@ def filterLinkedShapeFiles(def location, float distance, LinkedHashMap inputTabl
debug "Loading in the H2GIS database $outputTableName"
h2gis_datasource.execute("""DROP TABLE IF EXISTS $outputTableName ; CREATE TABLE $outputTableName as SELECT $formatting_geom,
CODE_INSEE FROM ${inputTables.commune} WHERE CODE_INSEE='$location' or lower(nom)='${location.toLowerCase()}'""".toString())
}else{
return
}
def count = h2gis_datasource."$outputTableName".rowCount
if (count > 0) {
Expand Down
15 changes: 15 additions & 0 deletions bdtopo/src/main/groovy/org/orbisgis/geoclimate/bdtopo/BDTopoV3Workflow.groovy
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,14 @@ def filterLinkedShapeFiles(def location, float distance, LinkedHashMap inputTabl
debug "Loading in the H2GIS database $outputTableName"
def communeColumns = h2gis_datasource.getColumnNames(inputTables.commune)
if(communeColumns.contains("INSEE_COM")) {
if(location.size()==3){
if(location[2]<100){
error("The distance to create a bbox from a point must be greater than 100 meters")
return
}
location = BDTopoUtils.bbox(location[0], location[1],location[2])
}

h2gis_datasource.execute("""DROP TABLE IF EXISTS $outputTableName ; CREATE TABLE $outputTableName as SELECT
ST_INTERSECTION(the_geom, ST_MakeEnvelope(${location[1]},${location[0]},${location[3]},${location[2]}, $sourceSRID)) as the_geom, INSEE_COM AS CODE_INSEE from ${inputTables.commune} where the_geom
&& ST_MakeEnvelope(${location[1]},${location[0]},${location[3]},${location[2]}, $sourceSRID) """.toString())
Expand Down Expand Up @@ -261,6 +269,13 @@ Integer loadDataFromPostGIS(Object input_database_properties, Object code, Objec
if (code in Collection) {
def communeColumns = h2gis_datasource.getColumnNames(commune_location)
if(communeColumns.contains("INSEE_COM")) {
if(code.size()==3){
if(code[2]<100){
error("The distance to create a bbox from a point must be greater than 100 meters")
return
}
code = BDTopoUtils.bbox(code[0], code[1],code[2])
}
String inputTableName = """(SELECT
ST_INTERSECTION(st_setsrid(the_geom, $commune_srid), ST_MakeEnvelope(${code[1]},${code[0]},${code[3]},${code[2]}, $commune_srid)) as the_geom, INSEE_COM as CODE_INSEE from $commune_location where
st_setsrid(the_geom, $commune_srid)
Expand Down
18 changes: 12 additions & 6 deletions osm/src/main/groovy/org/orbisgis/geoclimate/osm/WorkflowOSM.groovy
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -232,6 +232,9 @@ Map workflow(def input) {
error "The all parameter must be a boolean value"
return null
}

def osm_date= inputParameters.get("date")

def osm_size_area = inputParameters.get("area")

if (!osm_size_area) {
Expand All @@ -241,7 +244,6 @@ Map workflow(def input) {
error "The area of the bounding box to be extracted from OSM must be greater than 0 km²"
return null
}

def overpass_timeout = inputParameters.get("timeout")
if (!overpass_timeout) {
overpass_timeout = 900
Expand Down Expand Up @@ -360,7 +362,7 @@ Map workflow(def input) {
def logTableZones = postfix("log_zones")
Map osmprocessing = osm_processing(h2gis_datasource, processing_parameters, locations.findAll { it }, file_outputFolder, outputFileTables,
outputDatasource, outputTables, outputSRID, downloadAllOSMData, deleteOutputData, deleteOSMFile, logTableZones, osm_size_area,
overpass_timeout, overpass_maxsize)
overpass_timeout, overpass_maxsize, osm_date)
if (!osmprocessing) {
h2gis_datasource.save(logTableZones,"${file_outputFolder.getAbsolutePath() + File.separator}logzones.geojson", true)
return null
Expand Down Expand Up @@ -401,7 +403,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d
def outputSRID, def downloadAllOSMData,
def deleteOutputData, def deleteOSMFile,
def logTableZones, def bbox_size,
def overpass_timeout, def overpass_maxsize) {
def overpass_timeout, def overpass_maxsize,def overpass_date) {
//Store the zone identifier and the names of the tables
def outputTableNamesResult = [:]
//Create the table to log on the processed zone
Expand All @@ -421,7 +423,7 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d
error "Cannot find any geometry to define the zone to extract the OSM data"
return
}
def srid = h2gis_datasource.getSpatialTable(zone).srid
def srid = h2gis_datasource.getSrid(zone)
def reproject = false
if (outputSRID) {
if (outputSRID != srid) {
Expand All @@ -432,7 +434,11 @@ Map osm_processing(JdbcDataSource h2gis_datasource, def processing_parameters, d
}
//Prepare OSM extraction
//TODO set key values ?
def query = "[timeout:$overpass_timeout][maxsize:$overpass_maxsize]" + OSMTools.Utilities.buildOSMQuery(zones.envelope, null, OSMElement.NODE, OSMElement.WAY, OSMElement.RELATION)
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)

if (downloadAllOSMData) {
//Create a custom OSM query to download all requiered data. It will take more time and resources
Expand All @@ -441,7 +447,7 @@ 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]" + 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.envelope, keysValues, OSMElement.NODE, OSMElement.WAY, OSMElement.RELATION)
}

def extract = OSMTools.Loader.extract(query)
Expand Down
61 changes: 61 additions & 0 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 @@ -663,6 +663,7 @@ class WorflowOSMTest extends WorkflowAbstractTest {
"input" : [
"locations": [location],//["Pont-de-Veyle"],//[nominatim["bbox"]],//["Lorient"],
"area": 2800,
//"date":"2017-12-31T19:20:00Z",
/*"timeout":182,
"maxsize": 536870918,
"endpoint":"https://lz4.overpass-api.de/api"*/],
Expand Down Expand Up @@ -809,6 +810,66 @@ class WorflowOSMTest extends WorkflowAbstractTest {
assertTrue h2gis.firstRow("select count(*) as count from $building where HEIGHT_WALL>0 and HEIGHT_ROOF>0").count > 0
}

@Test
void testEstimateBuildingWithAllInputHeightFromPoint() {
String directory = folder.absolutePath + File.separator + "test_building_height"
File dirFile = new File(directory)
dirFile.delete()
dirFile.mkdir()
def osm_parmeters = [
"geoclimatedb": [
"folder": dirFile.absolutePath,
"name" : "geoclimate_chain_db",
"delete": false
],
"input" : [
"locations": [[43.726898,7.298452,100]]],
"output" : [
"folder": ["path" : directory,
"tables": ["building", "zone"]]],
"parameters" :
[
rsu_indicators: ["indicatorUse" : ["LCZ"]]
]
]
Map process = OSM.WorkflowOSM.workflow(osm_parmeters)
def tableNames = process.values()
def building = tableNames.building[0]
H2GIS h2gis = H2GIS.open("${directory + File.separator}geoclimate_chain_db;AUTO_SERVER=TRUE")
assertTrue h2gis.firstRow("select count(*) as count from $building where HEIGHT_WALL>0 and HEIGHT_ROOF>0").count > 0
}

@Disabled //Because it takes some time to build the OSM query
@Test
void testEstimateBuildingWithAllInputHeightDate() {
String directory = folder.absolutePath + File.separator + "test_building_height"
File dirFile = new File(directory)
dirFile.delete()
dirFile.mkdir()
def osm_parmeters = [
"geoclimatedb": [
"folder": dirFile.absolutePath,
"name" : "geoclimate_chain_db",
"delete": false
],
"input" : [
"locations": [[43.726898,7.298452,43.727677,7.299632]],
"date":"2015-12-31T19:20:00Z"],
"output" : [
"folder": ["path" : directory,
"tables": ["building", "zone"]]],
"parameters" :
["distance" : 0,
rsu_indicators: ["indicatorUse" : ["LCZ"]]
]
]
Map process = OSM.WorkflowOSM.workflow(osm_parmeters)
def tableNames = process.values()
def building = tableNames.building[0]
H2GIS h2gis = H2GIS.open("${directory + File.separator}geoclimate_chain_db;AUTO_SERVER=TRUE")
assertTrue h2gis.firstRow("select count(*) as count from $building where HEIGHT_WALL>0 and HEIGHT_ROOF>0").count > 0
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't be interesting to check two dates and verify that the number of buildings (for example) is different ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be but note that this test is disable because query OSM with date takes a while

}

@Test
void testOneSeaLCZ() {
String directory = folder.absolutePath + File.separator + "test_sea_lcz"
Expand Down
10 changes: 9 additions & 1 deletion osmtools/src/main/groovy/org/orbisgis/geoclimate/osmtools/utils/Utilities.groovy
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -638,7 +638,8 @@ Geometry geometryFromNominatim(def bbox) {
*
* @author Erwan Bocher (CNRS LAB-STICC)
*
* @param bbox 4 values to define a bbox
* @param bbox 4 values to define a bbox or 3 values to define a point (lat/long)
* with a distance around it, expressed in meters
* @return a JTS polygon
*/
Geometry geometryFromValues(def bbox) {
Expand All @@ -651,6 +652,13 @@ Geometry geometryFromValues(def bbox) {
if (bbox.size() == 4) {
return buildGeometry([bbox[1], bbox[0], bbox[3], bbox[2]]);
}
else if (bbox.size()==3){
if(bbox[2]<100){
error("The distance to create a bbox from a point must be greater than 100 meters")
return null
}
return getAreaFromPoint(bbox[1], bbox[0], bbox[2])
}
}


Expand Down
Loading