""" Services implemented: GMS/DEM/calls/conatenate POST: to concatenate raster """ from monitorT.task import Task import os import tempfile from tools.layer import Layer, OutLayer from GMS.tools.gdal_python import is_overlapping # Fix gdal 2.40 and 3.3 integration problems try: import gdal import ogr except ModuleNotFoundError: from osgeo import gdal, ogr class TaskCustom(Task): """ clase interna itasker """ CONTROL = True def post(task, user=None, aoi: Layer = None, orig_layer: Layer = None, layer_to_concat: Layer = None, output_layer: OutLayer = None): """ Servicio que concatena raster Parameters ---------- task user aoi: Layer Capa de entidades de interes orig_layer: Layer Capa de origen layer_to_concat: Layer Capa con la que concatenar output_layer: OutLayer Capa de salida Returns ------- output_layer: OutLayer devuelve el objeto OutLayer """ # Generating temporal files tmp_directory = tempfile.mkdtemp(dir='tmp') orig_crop = tempfile.NamedTemporaryFile(suffix='.tif', prefix='orto_', dir=tmp_directory).name concat_crop = tempfile.NamedTemporaryFile(suffix='.tif', prefix='concat_crop_', dir=tmp_directory).name concat_crop_changed = tempfile.NamedTemporaryFile( suffix='.tif', prefix='concat_crop_changed_', dir=tmp_directory).name result = tempfile.NamedTemporaryFile(suffix='.tif', prefix='result_crop_changed_', dir=tmp_directory).name # Get feature of AOI layer aoi_ds = ogr.Open(aoi.gdal_layer()) aoi_layer = aoi_ds.GetLayer() aoi_feature = aoi_layer.GetFeature(0) aoi_geom = aoi_feature.GetGeometryRef() # Generating vrt from folder vrt_string = 'gdalbuildvrt {result} {files}' # Vrto of ortos # if orig_layer.protocol.is_directory(orig_layer['source']): vrt_files = [] for file in orig_layer.protocol.list_files(orig_layer['source'], orig_layer.driver.format): lidar_parameter = orig_layer.parameters.copy() lidar_parameter['source'] = orig_layer.protocol.join( orig_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 = orig_layer.protocol.gdal_url(orig_layer['source']) # Generating vrt from folder vrt_string = 'gdalbuildvrt {result} {files}' # Vrto of ortos if layer_to_concat.protocol.is_directory(layer_to_concat['source']): vrt_files = [] for file in layer_to_concat.protocol.list_files( layer_to_concat['source'], layer_to_concat.driver.format): lidar_parameter = layer_to_concat.parameters.copy() lidar_parameter['source'] = layer_to_concat.protocol.join( layer_to_concat['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 to_concat_vrt = tempfile.NamedTemporaryFile(suffix='.vrt', prefix='to_concat_vrt_', dir=tmp_directory).name os.system(vrt_string.format(result=to_concat_vrt, files=' '.join(vrt_files))) else: to_concat_vrt = layer_to_concat.protocol.gdal_url( layer_to_concat['source']) # First crop orig by AOI gdal_warp_sentence = 'gdalwarp -q -multi -cutline "{cut_shape}" ' \ '-dstnodata 0 ' \ '-overwrite -crop_to_cutline ' \ '-r near "{source}" "{dst}" '\ .format(cut_shape=aoi.gdal_layer(), source=vrt_orto, dst=orig_crop) os.system(gdal_warp_sentence) # First crop orig by AOI gdal_warp_sentence = 'gdalwarp -q -multi -cutline "{cut_shape}" ' \ '-dstnodata 0' \ ' -overwrite -crop_to_cutline ' \ '-r near "{source}" "{dst}" '\ .format(cut_shape=aoi.gdal_layer(), source=to_concat_vrt, dst=concat_crop) os.system(gdal_warp_sentence) # Getting size x and size y, resizing to concat layer orig_ds = gdal.Open(orig_crop) gt = orig_ds.GetGeoTransform() pixel_size_x = gt[1] pixel_size_y = -gt[5] # ## Transforming to concat raster gdal_translate_command = "gdal_translate -tr {pixel_x} {pixel_y} " \ "{source} {output}" \ .format(pixel_x=pixel_size_x, pixel_y=pixel_size_y, source=concat_crop, output=concat_crop_changed) # ## Merging gdal_merge = "gdal_merge.py -ot Float32 -separate -o {result} {inputs}" \ .format(result=result, inputs=" ".join([orig_crop, concat_crop_changed])) os.system(gdal_translate_command) os.system(gdal_merge) output_layer.attach_sources((result, )) output_layer.residual_sources.append(tmp_directory) return output_layer