""" Services implemented: GMS/DEM/calls/LSMS POST: LSMS """ from monitorT.task import Task import os import tempfile from tools.layer import Layer, OutLayer # Fix of python 3.7 upgrade, OTB not working with this version try: import otbApplication except (ModuleNotFoundError, ImportError): pass from GMS.tools.gdal_python import is_overlapping, square_grid_polygon # Fix gdal 2.40 and 3.3 integration problems try: import ogr except ModuleNotFoundError: from osgeo import ogr from MG.tools.postgis_api import PostGis FIX_FACTOR = 0.7 RANGER_MAX = 35 RANGER_MIN = 5 def algorithm_quality(confident): """ Service to get an algorith of quality thought a level of confidence Parameters ---------- confident: int Level of confidence Returns ------- Float level of quality """ return float(confident)*0.5+FIX_FACTOR*0.5 class TaskCustom(Task): """ clase interna itasker """ CONTROL = True def post(task, user=None, aoi: Layer = None, raster_layer: Layer = None, output_layer: OutLayer = None, spatialr: float=None, ranger: float = None, tiled: float = 0): """ Service to get the LSMS Parameters ---------- task user aoi: Layer Layer with entities of interest raster_layer: Layer Raster layer output_layer: OutLayer Output layer spatial_r: float Spatial Reference ranger: float Range tiled: float Default 0 Returns ------- output_layer: OutLayer """ # Generating temporal files tmp_directory = tempfile.mkdtemp(dir='tmp') cut_kml_temp = tempfile.NamedTemporaryFile(suffix='.kml', prefix='cut_', dir=tmp_directory).name # Get feature of AOI layer aoi_ds = ogr.Open(aoi.gdal_layer()) aoi_layer = aoi_ds.GetLayer() aoi_crs = aoi_layer.GetSpatialRef() aoi_feature = aoi_layer.GetFeature(0) aoi_geom = aoi_feature.GetGeometryRef() # Generatin vrt from folder vrt_string = 'gdalbuildvrt {result} {files}' # Vrto of ortos ###################################################### if raster_layer.protocol.is_directory(raster_layer['source']): vrt_files = [] for file in raster_layer.protocol.list_files( raster_layer['source'], raster_layer.driver.format): lidar_parameter = raster_layer.parameters.copy() lidar_parameter['source'] = raster_layer.protocol.join( raster_layer['source'], file) new_layer = Layer(user, lidar_parameter) gdal_url = new_layer.gdal_layer() if is_overlapping(gdal_url, aoi_geom): vrt_files.append(gdal_url) del new_layer vrt_orto = tempfile.NamedTemporaryFile(suffix='.vrt', prefix='lidar_', dir=tmp_directory).name os.system(vrt_string.format(result=vrt_orto, files=' '.join(vrt_files))) else: vrt_orto = raster_layer.protocol.gdal_url(raster_layer['source']) # Tiling analysis area if float(tiled): grid_shape = tempfile.NamedTemporaryFile(suffix=".shp", dir=tmp_directory).name tiles = square_grid_polygon(float(tiled), aoi_geom, gen_grid=grid_shape, offset=5) else: tiles = [aoi_geom] # Running analysis over tiles result = [] postgis_obj = PostGis('public') indexing_table = """CREATE INDEX \"{schema}_{table}\" ON "{schema}".\"{table}\" USING gist({geometry_field})""" for tile_index, tile in enumerate(tiles): orto_png_path = tempfile.NamedTemporaryFile(suffix=".tif", prefix='tile_orto_', dir=tmp_directory).name # First crop by AOI kml_driver = ogr.GetDriverByName('KML') shp_datasource = kml_driver.CreateDataSource(cut_kml_temp) layer = shp_datasource.CreateLayer(cut_kml_temp, geom_type=ogr.wkbPolygon, srs=tile.GetSpatialReference()) feature = ogr.Feature(layer.GetLayerDefn()) feature.SetGeometry(tile) layer.CreateFeature(feature) feature.Destroy() del shp_datasource gdal_warp_sentence = 'gdalwarp -q -multi -cutline "{cut_shape}" ' \ '-dstnodata 0 -overwrite -crop_to_cutline ' \ '-r near "{source}" "{dst}" '\ .format(cut_shape=cut_kml_temp, source=vrt_orto, dst=orto_png_path) os.system(gdal_warp_sentence) shp_geometry = tempfile.NamedTemporaryFile(prefix='OTB_tile_{}_'.format( tile_index), suffix='.shp', dir=tmp_directory).name # (https://www.orfeo-toolbox.org/CookBook/recipes/improc.html#large-scale-mean-shift-lsms-segmentation) app = otbApplication.Registry.CreateApplication("LargeScaleMeanShift") app.SetParameterString("in", orto_png_path) app.SetParameterInt("spatialr", int(spatialr)) app.SetParameterFloat("ranger", int(ranger)) app.SetParameterInt("tilesizex", 800) app.SetParameterInt("tilesizey", 800) app.SetParameterInt("minsize", 10) app.SetParameterString("mode.vector.out", shp_geometry) app.ExecuteAndWriteOutput() command = 'ogr2ogr -f "{}" -lco OVERWRITE=yes -lco LAUNDER=no PG:"host={} port={} dbname={} user={} password={} ' \ 'schemas={}" -t_srs EPSG:{} {}'.format( 'PostgreSQL', postgis_obj.ip, postgis_obj.port, postgis_obj.dbname, postgis_obj.user, postgis_obj.passw, user, aoi_crs.GetAuthorityCode(None), shp_geometry) os.system(command) otb_table = shp_geometry.split(os.sep)[-1].split('.')[0] sql_command = indexing_table.format(schema=user, table=otb_table, geometry_field='wkb_geometry') postgis_obj.send_sql_command(sql_command) result.append(otb_table) ################################################################### # Creating a master table create_master_table = "create table \"{schema}\".\"{master_table}\" " \ "(ogc_fid serial NOT NULL," \ "label INTEGER," \ "\"nbPixels\" double precision," \ "\"meanB0\" double precision," \ "\"meanB1\" double precision," \ "\"meanB2\" double precision," \ "\"varB0\" double precision," \ "\"varB1\" double precision," \ "\"varB2\" double precision," \ "\"wkb_geometry\" geometry);" master_table = postgis_obj.get_unique_name(user, 'master_otb_') postgis_obj.send_sql_command(create_master_table.format( schema=user, master_table=master_table)) merge_table = "INSERT INTO \"{schema}\".\"{master_table}\" " \ "(SELECT (select * from nextval" \ "('\"{schema}\".\"{master_table}_ogc_fid_seq\"')), " \ "label, \"nbPixels\", \"meanB0\", \"meanB1\"," \ "\"meanB2\", \"varB0\", \"varB1\", \"varB2\"," \ " \"wkb_geometry\" FROM \"{schema}\".\"{table_in}\");" for otb_table in result: postgis_obj.send_sql_command(merge_table.format( schema=user, master_table=master_table, table_in=otb_table)) postgis_obj.delete_tables(user, otb_table, ) output_layer.residual_tables.extend(result) output_layer.attach_sources((master_table, )) return output_layer