j3r3m1
commented
7 months ago Which area were you expecting to obtain ? You probably know that you can check what is the OSM area corresponding to your location using Nominatim. Here is what you get for your string location: https://nominatim.openstreetmap.org/ui/search.html?q=Senlis%2C+Oise%2C+Hauts-de-France%2C+France+m%C3%A9tropolitaine%2C+60300%2C+France
Note that the first geometry will be taken. When I use the GeoClimate workflow with the corresponding string I get the same area.
balglave
Hi, I want to run the following groovy code with line
groovy geoclimate_with_urban_sprawl.groovy ["Senlis, Oise, Hauts-de-France, France métropolitaine, 60300, France"]But I get a totally different area. Could you tell me how to solve the problem.
Groovy code
@GrabResolver(name='orbisgis', root='https://oss.sonatype.org/content/repositories/snapshots/') @Grab(group='org.orbisgis.geoclimate', module='geoclimate', version='0.0.2-SNAPSHOT')
/**
import org.orbisgis.geoclimate.osm.OSM import org.orbisgis.data.H2GIS
String location = args[0]
//Directory to store all the results String outputDirectory ="/home/balglave/Desktop/Research/geoclimate/res/"
if(!outputDirectory){ println "The output directory to store the result cannot be null or empty" }
File dirFile = new File(outputDirectory) if(!dirFile.exists()){ println "Create the output directory because it doesn't exist" dirFile.mkdir() }
//Compute the BBOX of the ICU layer
def local_database_name = "geoclimate_sprawl"
def h2gis_db_parameters = [ "folder": outputDirectory, "name" : "${local_database_name};AUTO_SERVER=TRUE".toString(), "delete": false ]
if(!location){ println("The location value cannot be null") }
//Create before the local H2GIS database that will be used to compute the BBOX of the area H2GIS h2GIS = H2GIS.open(h2gis_db_parameters.folder + File.separator + h2gis_db_parameters.name) if(h2GIS==null){ println("Cannot create the local H2GIS database") return }
def outputTableNames = runGeoclimate(h2gis_db_parameters, location, outputDirectory )
if(!outputTableNames){ println("Cannot compute the GeoClimate indicators") }
def rsu_lcz = outputTableNames.rsu_lcz
//Extract sprawl area String urbanSprawlH2GIS = computeUrbanSprawl(h2GIS, rsu_lcz )
if(!urbanSprawlH2GIS){ println("Cannot compute the urban sprawl") return }
//Save the data in a geojson file h2GIS.save(urbanSprawlH2GIS, outputDirectory+File.separator+"osm_"+location+ File.separator+"urban_sprawl.geojson", true)
/**
@return */ def runGeoclimate(def h2gis_db_parameters, def zone, def outputDirectory) {
if (!zone) { println "The name of the zone to process cannot be null or empty" }
/*================================================================================
Exemple with OSM input areas configuration */
def input = ["locations": [zone], "delete": true, "area": 2800]
/*================================================================================
/*================================================================================
WORKFLOW PARAMETERS */ def workflow_parameters = [ "description" : "Run the Geoclimate chain and export result to a folder", "geoclimatedb": h2gis_db_parameters, "input" : input, "output" : output, "parameters" : [ "rsu_indicators" : [ "indicatorUse": ["LCZ"] ], "worldpop_indicators": true ] ]
Map process = OSM.workflow(workflow_parameters) return process.values().collectEntries() }
/**
@return */ String computeUrbanSprawl(H2GIS h2GIS, String lczTableName ){
def mergingLcz = "merging_lcz"
h2GIS.execute(""" DROP TABLE IF EXISTS $mergingLcz; CREATE TABLE $mergingLcz as select EXPLOD_ID as id_rsu, the_geom from ST_EXPLODE('( select st_removeholes(st_union(st_accum(the_geom))) as the_geom FROM $lczTableName where LCZ_PRIMARY IN (1,2,3,4,5,6,7,8,9,10,105))')
""".toString())
//Create a grid to generalize the urban areas h2GIS.execute(""" DROP TABLE IF EXISTS grid; CREATE TABLE grid as select * from ST_MAKEGRID('$mergingLcz', 100,100); """.toString())
//Compute the percentage by grid h2GIS.execute(""" DROP TABLE IF EXISTS grid_lcz, grid_lcz_area, grid_raster; CREATE TABLE grid_lcz as select st_area(st_intersection(a.the_geom, b.the_geom)) as area, a.id from grid as a, $mergingLcz as b where a.the_geom && b.the_geom and st_intersects(a.the_geom, b.the_geom); CREATE TABLE grid_lcz_area as select sum(area) as sum_lcz, id from grid_lcz group by id; CREATE INDEX ON grid_lcz_area (ID); CREATE INDEX ON grid (ID); CREATE TABLE grid_raster as select a.id, a.the_geom from grid as a left join grid_lcz_area as b on a.id=b.id where sum_lcz/(100*100) > 0.4 """.toString())
DROP TABLE IF EXISTS $urban_sprawl; CREATE TABLE $urban_sprawl as select EXPLOD_ID as id_rsu, the_geom from ST_EXPLODE('( select st_removeholes(st_union(st_accum(the_geom))) as the_geom FROM grid_raster)') ; DROP TABLE IF EXISTS grid_lcz, grid_raster, grid, grid_lcz_area, mergingLcz; """.toString())
}