import sys import numpy as np import laspy as lp import geopandas as gpd import pdal import os import json import subprocess import datetime import traceback import time import fnmatch import csv from shapely.geometry import Polygon #entrada: path a carpeta con lasz desde el directorio de ejecucion # path a carpeta con .shp desde el directorio de ejecucion start_time = time.time() #entrada='./DATA4/2020-21_AND_PNOA_LIDAR/bv/fase2/' #shape_vial='./shapes/concat_urbana_andalucia.shp' #shape_embalses='./shapes/orillas_embalses_julio.shp' #shape_edifs='./shapes/JULIO_ALL_EDIF.shp' #shape_industrial='./shapes/BTN0513S_INS_IND_2.5.shp' #salida='./DATA4/2020-21_AND_PNOA_LIDAR/bv/fase2_procesados/' def cortarshape(entrada,salida,archivo,shape,tipo): lidar=lp.read(entrada+archivo) laspyxmax = round(lidar.header.x_max) #print('\nAltitud maxima x:', laspyxmax) laspyymax = round(lidar.header.y_max) #print('\nAlitud maxima y:', laspyymax) laspyxmin = round(lidar.header.x_min) #print('\nAltitud minima x:', laspyxmin) laspyymin = round(lidar.header.y_min) #print('\nAltitud minima y:', laspyymin) cuadrado_seleccionado = Polygon([(laspyxmin,laspyymin), (laspyxmin,laspyymax), (laspyxmax,laspyymax), (laspyxmax,laspyymin), (laspyxmin,laspyymin)]) shape_2x2 = gpd.read_file(shape,cuadrado_seleccionado) try: shape_2x2 = shape_2x2.dissolve(by=None) except: #shape2 = shape2.explode(column=None) polygons = shape_2x2['geometry'] polygons=polygons.make_valid() shape_2x2 = shape_2x2.set_geometry(polygons.buffer(0.01,1)) shape_2x2 = shape_2x2.dissolve(by=None) shape_2x2 = shape_2x2.assign(ruido=7) shape_2x2 = shape_2x2.assign(suelo=2) shape_2x2 = shape_2x2.assign(edificio=6) shape_2x2 = shape_2x2.assign(industrial=21) pathshape=salida+'2x2_'+tipo+'_'+archivo.rsplit('.',1)[0]+'.shp' if shape_2x2.empty: return None shape_2x2.to_file(pathshape,driver="ESRI Shapefile",crs="EPSG:25830") return pathshape def pdal_overlay(origen,archivo,destino,shape_in): nombrearchivo=archivo.rsplit('.',1)[0] formato=archivo.rsplit('.',1)[1] filename=destino+nombrearchivo+'.'+formato inn=origen+archivo jsonOverlay=[] in1 = { "filename":origen+archivo, "tag":"t1" } hag = { "type":"filters.hag_nn", "tag":"t1h" } coplanar = { "type":"filters.planefit", "knn":"8" } overlay = { "type":"filters.overlay", "dimension":"Classification", "datasource":shape_in[0], "where":"(Classification==4 || Classification==5 || Classification==6) && HeightAboveGround<=4 && HeightAboveGround>=0.6", "column":"ruido", } overlay2 = { "type":"filters.overlay", "dimension":"Classification", "datasource":shape_in[0], "where":"Classification==3", "column":"suelo" } overlay3 = { "type":"filters.overlay", "dimension":"Classification", "datasource":shape_in[1], "where":"Classification==9", "column":"suelo", } overlay4 = { "type":"filters.overlay", "dimension":"Classification", "datasource":shape_in[2], "where":"(Classification==3 || Classification==4 || Classification==5) && HeightAboveGround>=4 && PlaneFit>0.6", "column":"ruido", } overlay4_2 = { "type":"filters.overlay", "dimension":"Classification", "datasource":shape_in[2], "where":"(Classification==3 || Classification==4 || Classification==5) && HeightAboveGround>=4 && PlaneFit<=0.6", "column":"edificio", } overlay5 = { "type":"filters.overlay", "dimension":"Classification", "datasource":shape_in[3], "where":"Classification==6", "column":"industrial", } writers = { "type": "writers.las", "filename": filename } jsonOverlay.append(in1) jsonOverlay.append(hag) jsonOverlay.append(coplanar) if shape_in[0] is not None: jsonOverlay.append(overlay) jsonOverlay.append(overlay2) if shape_in[1] is not None: jsonOverlay.append(overlay3) if shape_in[2] is not None: jsonOverlay.append(overlay4) jsonOverlay.append(overlay4_2) if shape_in[3] is not None: jsonOverlay.append(overlay5) jsonOverlay.append(writers) peticion = json.dumps(jsonOverlay) pipeline = pdal.Pipeline(peticion) count = pipeline.execute() arrays = pipeline.arrays metadata = pipeline.metadata log = pipeline.log def cambiar_clasif(entrada,salida,vial,embalses,edifs,industrial,logs,id): try: file = open(logs+"log_filtrado_clasif_"+id+".txt","w") if(entrada[len(entrada)-1]!='/'): log='Path entrada mal especificado. Ha de acabar en "/"' file.write(log) file.close() print(log) return if(salida[len(salida)-1]!='/'): log='Path salida mal especificado. Ha de acabar en "/"' file.write(log) file.close() print(log) return nohaylases=True for archivo in os.listdir(entrada): try: if fnmatch.fnmatch(archivo,'*.las') or fnmatch.fnmatch(archivo,'*.laz'): log='\n#####################'+entrada+archivo+'#####################\n' file.write(log) print(log) nohaylases = False pathshape2x2vial=cortarshape(entrada,salida,archivo,vial,'vial') pathshape2x2embalses=cortarshape(entrada,salida,archivo,embalses,'embalses') pathshape2x2edifs=cortarshape(entrada,salida,archivo,edifs,'edifs') pathshape2x2industrial=cortarshape(entrada,salida,archivo,industrial,'industrial') shapess=[pathshape2x2vial,pathshape2x2embalses,pathshape2x2edifs,pathshape2x2industrial] pdal_overlay(entrada,archivo,salida,shapess) except Exception: traceback.print_exc() file2=open(logs+"log_error_filtrado_clasif"+id+".txt","a") file2.write(archivo+"\n") file.write(traceback.format_exc()) file2.write(traceback.format_exc()) continue if nohaylases: log= '\n #################### No hay .las/.laz en '+entrada+' #################### \n' file.write(log) print(log) log="--- "+str(time.time() - start_time)+" seconds ---" file.write(log) print(log) file.close() return log="\n--- "+str(time.time() - start_time)+" seconds ---\n" file.write(log) print(log) file.close() except Exception: file = open(logs+"log_error_filtrado_clasif"+id+".txt","a") traceback.print_exc() file.write(traceback.format_exc()) log="\n--- "+str(time.time() - start_time)+" seconds ---\n" file.write(log) print(log) file.close() cambiar_clasif('./in/','./out/',sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],'./logs/',sys.argv[5])