Ráster de atenuación de ruido obtenido mediante PyQGIS para una fuente lineal

En el post anterior se consideró la obtención de un ráster de atenuación de ruido obtenido mediante PyQGIS para una fuente puntual. En éste se tomará en cuenta una fuente lineal de ruido; tal como corresponde a una autopista con mucho tráfico donde los valores pueden ser tan altos como de 90 dBA (Id1) a una distancia d1 de sólo 1 m.

Para una fuente lineal la Intensidad tiene una reducción de 3 dBA al doblar la distancia, por tanto, la fórmula empírica que da cuenta de las Intensidades como función de las distancias se escribe finalmente como:

Idn = Id1 - (3./log10(2)).(log10(dn)-log10(d1))

Para estimar todas las distancias pertinentes en el ráster (dn), se considera aquella que existe entre el centro de la celda (p1) y el punto que resulta de la intersección del segmento más corto que intersecta a la recta y que parte de p1. Estás distancias dn, conjuntamente con los valores de 90 dBA (Id1) a una distancia d1 de 1 m (para una autopista), se emplean para estimar las Intensidades y producir un ráster de 200 filas por 200 columnas empleando los parámetros de Geo transformación de un ráster de referencia. El código completo es el siguiente:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80

from osgeo import gdal, osr from math import log10 import numpy as np registry = QgsMapLayerRegistry.instance() linear_noise_source = registry.mapLayersByName('linear_noise_source') feat_line = linear_noise_source[0].getFeatures().next() utah_lake_part = registry.mapLayersByName('utah_lake_part') rprovider = utah_lake_part[0].dataProvider() raster_filename = rprovider.dataSourceUri() rows = utah_lake_part[0].height() cols = utah_lake_part[0].width() xsize = utah_lake_part[0].rasterUnitsPerPixelX() ysize = utah_lake_part[0].rasterUnitsPerPixelY() dataset = gdal.Open(raster_filename) band = dataset.GetRasterBand(1) transform = dataset.GetGeoTransform() #highway noise parameters d1 = 1 #distance (meters) Id1 = 90 #Intensity (decibels) values = [ [] for i in range(rows) ] x = transform[0] + xsize/2 y = transform[3] - ysize/2 points = [] #for verifiying purposes points_lines = [] distances = [] for i in range(rows): for j in range(cols): pt_geom = feat_line.geometry() pt_line = pt_geom.closestSegmentWithContext(QgsPoint(x,y))[1] dn = pt_geom.distance(QgsGeometry.fromPoint(QgsPoint(x,y))) distances.append(dn) points.append((x,y)) points_lines.append(pt_line) Idn = Id1 - (3./log10(2))*(log10(dn)-log10(d1)) values[i].append(Idn) x += xsize x = transform[0] + xsize/2 y -= ysize data = np.array(values) # Create gtif file driver = gdal.GetDriverByName("GTiff") output_file = "/home/zeito/pyqgis_data/noise_linear_raster.tif" dst_ds = driver.Create(output_file, cols, rows, 1, gdal.GDT_Float32) #writting output raster dst_ds.GetRasterBand(1).WriteArray( data ) #setting extension of output raster # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution dst_ds.SetGeoTransform(transform) wkt = dataset.GetProjection() # setting spatial reference of output raster srs = osr.SpatialReference() srs.ImportFromWkt(wkt) dst_ds.SetProjection( srs.ExportToWkt() ) #Close output raster dataset dataset = None dst_ds = None

Al ejecutar el código anterior en la Python Console de QGIS se obtiene lo siguiente:

Valores de 60 dBA (ruido tolerable) se alcanzan en la zona difusa que corresponde al color azul.