import geopandas as gpd #pip install geopandas import fnmatch from shapely.geometry import Polygon import shapefile # Hay que instalar este paquete -> pip install pyshp / conda install -c conda-forge pyshp import os import json import pdal # conda install -c conda-forge pdal python-pdal gdal / pip install pdal / sudo apt-get install python3-pdal /apt-get install pdal import traceback import laspy as lp # conda install -c conda-forge laspy | pip install laspy | pip install laspy[lazrs,laszip] import numpy as np import subprocess import datetime import time import sys import re import csv from osgeo import gdal #conda install -c conda-forge gdal directorio_entrada = sys.argv[1] directorio_salida = sys.argv[2] directorio_logs = sys.argv[3] id_proceso = sys.argv[4] directorio_temporal = sys.argv[5] start_time = time.time() def get_bounds(entrada,raster): #point_cloud = lp.read(entrada+raster) #print('\n segun laspy: \n') #xmax = round(point_cloud.header.x_max) #print('\nAltitud maxima x:', laspyxmax) #ymax = round(point_cloud.header.y_max) #print('\nAlitud maxima y:', laspyymax) #xmin = round(point_cloud.header.x_min) #print('\nAltitud minima x:', laspyxmin) #ymin = round(point_cloud.header.y_min) #print('\nAltitud minima y:', laspyymin, '\n') coordenadas = re.search('[0-9][0-9][0-9]-[0-9][0-9][0-9][0-9]',raster).group(0) xmin = int(coordenadas.split('-')[0]) * 1000 xmax = xmin + 2000 ymax = int(coordenadas.split('-')[1]) * 1000 ymin = ymax - 2000 return (xmin,ymin,xmax,ymax) def pdal_merge_and_DTM3x3(directoriorigen,directoriodestino,formato,rastertiff): jsonDTM=[directoriorigen+'*.'+formato] archivo=rastertiff+"."+formato coordenadas = re.search('[0-9][0-9][0-9]-[0-9][0-9][0-9][0-9]',archivo).group(0) xmin = int(coordenadas.split('-')[0]) * 1000 xmax = xmin + 2000 ymax = int(coordenadas.split('-')[1]) * 1000 ymin = ymax - 2000 xmin=str(xmin-10) ymin=str(ymin-10) destino=directoriodestino+'output3x3'+rastertiff+'_DTM.tif' merge = { "type": "filters.merge" } expression = { "type":"filters.expression", "expression":"Classification == 2" } delaunay = { "type": "filters.delaunay" } faceraster = { "type": "filters.faceraster", "resolution": 1.0, "width": 2020, "height": 2020, "origin_x": xmin, "origin_y": ymin } writers = { "filename": destino, "type": "writers.raster" } jsonDTM.append(merge) jsonDTM.append(expression) jsonDTM.append(delaunay) jsonDTM.append(faceraster) jsonDTM.append(writers) peticion = json.dumps(jsonDTM) pipeline = pdal.Pipeline(peticion) count = pipeline.execute() arrays = pipeline.arrays #print('arrays '+ str(arrays)) metadata = pipeline.metadata #print('metadata '+ str(metadata)) log = pipeline.log def proceso_debug(entrada): path_to_join = entrada listado = os.listdir(entrada) archivos = [] for element in listado: ruta_completa = os.path.join(path_to_join,element) print(f'-- El element es {ruta_completa}') if os.path.isdir(ruta_completa): print("Hemos encontrado una carpeta") sub_archivos = os.listdir(ruta_completa) archivos.extend([os.path.join(ruta_completa,sub_archivo) for sub_archivo in sub_archivos]) # archivos.extend(sub_archivos) if os.path.isfile(ruta_completa): print("Hemos encontrado un archivo") archivos.append(ruta_completa) print("Listado final de archivos") for archivo in archivos: path_archivo,archivo = archivo.rsplit("/",1) print(f'-- El archivo es {archivo}') print(f'++ El path del archivo es {path_archivo}') def tifcortado(entrada,salida,logs,id,temporal): try: if(temporal==''): temporal='./' try: os.mkdirs(temporal) except: pass try: os.mkdirs(logs) except: pass filogs=logs+"log_DTM_"+id+".txt" file = open(filogs,"a") if(salida[len(salida)-1]!='/'): log='Path entrada mal especificado. Ha de acabar en "/"' file.write(log) file.close() print(log) return carpeta_temporal_DTM = temporal+'temporalDTM'+id carpeta_temporal_DTM2 = temporal+'temporalDTM2'+id try: os.mkdir(carpeta_temporal_DTM) os.mkdir(carpeta_temporal_DTM2) # os.mkdir(temporal+'temporalDTM'+id) # os.mkdir(temporal+'temporalDTM2'+id) except: log='\nYa están creadas las carpetas\n' print(log) file.write(log) subprocess.run('cp -r '+entrada+'* '+carpeta_temporal_DTM+'/',shell = True) # subprocess.run('cp -r '+entrada+'* '+temporal+'temporalDTM'+id+'/',shell = True) path_to_join = carpeta_temporal_DTM listado = os.listdir(path_to_join) archivos = [] for element in listado: ruta_completa = os.path.join(path_to_join,element) # print(f'-- El element es {ruta_completa}') if os.path.isdir(ruta_completa): # print("Hemos encontrado una carpeta") sub_archivos = os.listdir(ruta_completa) archivos.extend([os.path.join(ruta_completa,sub_archivo) for sub_archivo in sub_archivos]) # archivos.extend(sub_archivos) if os.path.isfile(ruta_completa): # print("Hemos encontrado un archivo") archivos.append(ruta_completa) # archivos=os.listdir(temporal+'temporalDTM'+id) for archivo in archivos: path_archivo,archivo = archivo.rsplit("/",1) # print(f"-- El path {path_archivo}") # print(f"++ El archivo {archivo}") # print(carpeta_temporal_DTM) if path_archivo == carpeta_temporal_DTM: # print("El archivo se encuentra dentro de la propia carpeta temporal DTM") carpeta_temporal_DTM_this_file = carpeta_temporal_DTM carpeta_temporal_DTM2_this_file = carpeta_temporal_DTM2 else: # print("El archivo NO se encuentra dentro de la propia carpeta temporal DTM") carpeta_temporal_DTM_this_file = f"{carpeta_temporal_DTM}/{path_archivo.split(carpeta_temporal_DTM+'/',1)[1]}" carpeta_temporal_DTM2_this_file = f"{carpeta_temporal_DTM2}/{path_archivo.split(carpeta_temporal_DTM+'/',1)[1]}" try: os.mkdirs(carpeta_temporal_DTM2_this_file) except: pass # print(carpeta_temporal_DTM_this_file) # print(carpeta_temporal_DTM2_this_file) # sys.exit("Fuerzo el terminar para probar una cosa") # print('Despues del continue, esto no deberia de aparecer') try: coordenadas = re.search('[0-9][0-9][0-9]-[0-9][0-9][0-9][0-9]',archivo).group(0) xmin = int(coordenadas.split('-')[0]) * 1000 xmax = xmin + 2000 ymax = int(coordenadas.split('-')[1]) * 1000 ymin = ymax - 2000 log="3x3 para archivo "+path_archivo+'/'+archivo # log="3x3 para archivo ./temporalDTM"+id+'/'+archivo file.write(log) print(log) for vecino in archivos: path_vecino,vecino = vecino.rsplit("/",1) try: coordenadasvecino = re.search('[0-9][0-9][0-9]-[0-9][0-9][0-9][0-9]',vecino).group(0) xminvecino = int(coordenadasvecino.split('-')[0]) * 1000 xmaxvecino = xminvecino + 2000 ymaxvecino = int(coordenadasvecino.split('-')[1]) * 1000 yminvecino = ymaxvecino - 2000 if abs(xmin-xminvecino)<=2000 and abs(ymax-ymaxvecino)<=2000: log="Vecino "+path_vecino+'/'+vecino+' movido a '+carpeta_temporal_DTM2_this_file+'/'+'\n' # log="Vecino "+'./temporalDTM'+id+'/'+vecino+' movido a '+'./temporalDTM2'+id+'/'+'\n' file.write(log) print(log) subprocess.run('mv '+path_vecino+'/'+vecino+' '+carpeta_temporal_DTM2_this_file+'/',shell = True) # subprocess.run('mv '+temporal+'temporalDTM'+id+'/'+vecino+' '+temporal+'temporalDTM2'+id+'/',shell = True) except Exception: log="Error moviendo vecino/s" print(log) file.write(log) traceback.print_exc() file.write(traceback.format_exc()) continue try: formato=archivo.rsplit('.',1)[1] rastertiff=archivo.rsplit('.',1)[0] pdal_merge_and_DTM3x3(carpeta_temporal_DTM2_this_file+'/',carpeta_temporal_DTM2_this_file+'/',formato,rastertiff) # pdal_merge_and_DTM3x3(temporal+'temporalDTM2'+id+'/',temporal+'temporalDTM2'+id+'/',formato,rastertiff) except Exception: subprocess.run('mv '+carpeta_temporal_DTM2_this_file+'/* '+carpeta_temporal_DTM_this_file+'/',shell = True) # subprocess.run('mv '+temporal+'temporalDTM2'+id+'/* '+temporal+'temporalDTM'+id+'/',shell = True) log="\n error durante dtm 3x3 raster "+path_archivo+"/"+archivo+" sigo \n" print(log) file.write(log) traceback.print_exc() file.write(traceback.format_exc()) continue input_raster_original = gdal.Open(carpeta_temporal_DTM2_this_file+'/output3x3'+rastertiff+'_DTM.tif') # input_raster_original = gdal.Open(temporal+'temporalDTM2'+id+'/output3x3'+rastertiff+'_DTM.tif') log='\n'+carpeta_temporal_DTM2_this_file+'/output3x3'+rastertiff+'_DTM.tif\n' # log='\n./temporalDTM2'+id+'/output3x3'+rastertiff+'_DTM.tif\n' file.write(log) bounds=get_bounds(entrada,archivo) print('\n rastertiff \n',rastertiff,'\n') ds=gdal.Warp(salida+rastertiff+'.tif',input_raster_original,format='GTiff',cropToCutline=False,outputBounds=bounds,srcSRS="EPSG:25830",dstSRS="EPSG:25830",resampleAlg="near") log="\n"+salida+rastertiff+'.tif\n' file.write(log) subprocess.run('mv '+carpeta_temporal_DTM2_this_file+'/* '+carpeta_temporal_DTM_this_file+'/',shell = True) # subprocess.run('mv '+temporal+'temporalDTM2'+id+'/* '+temporal+'temporalDTM'+id+'/',shell = True) except Exception: file.write(traceback.format_exc()) traceback.print_exc() continue try: subprocess.run('mv '+carpeta_temporal_DTM+'/* '+salida, shell = True) subprocess.run('rm '+salida+'output*',shell=True) subprocess.run('rm -r '+carpeta_temporal_DTM2,shell=True) subprocess.run('rm -r '+carpeta_temporal_DTM2,shell=True) # subprocess.run('mv '+temporal+'temporalDTM'+id+'/*.tif '+salida, shell = True) # subprocess.run('rm '+salida+'output*',shell=True) # subprocess.run('rm -r '+temporal+'temporalDTM2'+id,shell=True) # subprocess.run('rm -r '+temporal+'temporalDTM2'+id,shell=True) except Exception: log=('\nproblema reorganizando directorios\n') print(log) file.write(log) traceback.print_exc() file.write(traceback.format_exc()) return except Exception: file = open(logs+"log_error_DTM_"+id+".txt","a") traceback.print_exc() file.write(traceback.format_exc()) log="\n--- "+str(time.time() - start_time)+" seconds ---\n" file.write(log) file.close() tifcortado(directorio_entrada, directorio_salida, directorio_logs, id_proceso, directorio_temporal) #proceso_debug(directorio_entrada)