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 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 tifcortado(entrada,salida,logs,id,temporal): try: if(temporal==''): temporal='./' 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 try: 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 '+entrada+'* '+temporal+'temporalDTM'+id+'/',shell = True) archivos=os.listdir(temporal+'temporalDTM'+id) for archivo in archivos: 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 ./temporalDTM"+id+'/'+archivo file.write(log) print(log) for vecino in archivos: 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 "+'./temporalDTM'+id+'/'+vecino+' movido a '+'./temporalDTM2'+id+'/'+'\n' file.write(log) print(log) 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(temporal+'temporalDTM2'+id+'/',temporal+'temporalDTM2'+id+'/',formato,rastertiff) except Exception: subprocess.run('mv '+temporal+'temporalDTM2'+id+'/* '+temporal+'temporalDTM'+id+'/',shell = True) log="\n error durante dtm 3x3 raster "+archivo+" sigo \n" print(log) file.write(log) traceback.print_exc() file.write(traceback.format_exc()) continue input_raster_original = gdal.Open(temporal+'temporalDTM2'+id+'/output3x3'+rastertiff+'_DTM.tif') 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 '+temporal+'temporalDTM2'+id+'/* '+temporal+'temporalDTM'+id+'/',shell = True) except Exception: file.write(traceback.format_exc()) traceback.print_exc() continue try: subprocess.run('mv '+temporal+'temporalDTM'+id+'/*.tif '+salida, shell = True) subprocess.run('rm '+salida+'output*',shell=True) subprocess.run('rm -r '+temporal+'temporalDTM'+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('./in/','./out/','./logs/',sys.argv[1],'./temporal/')