domingo, 18 de junio de 2017

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

La atenuación del ruido, tal como ya se ha señalado con anterioridad, depende de si la fuente emisora es puntual o lineal. En ausencia de obstáculos, se refiere que cada vez que se dobla la distancia con relación a una fuente puntual el sonido se atenua en 6 dBA (decibelios), mientras que en el segundo caso, cada vez que se dobla la distancia con relación a una fuente lineal (por ejemplo, una autopista con mucho tráfico automotor) la atenuación del sonido será sólo de 3 dBA. Generalmente, los reportes de Tabla para la intensidad del sónido (Id1) producida por una determinada fuente puntual se hacen con base a una distancia fija que denominaremos d1. Por tanto, para una taladradora de construcción, es común que se reporte que produce una intensidad de 110 dBA (Id1) a una distancia d1 de 10 m.


En la referencia anterior, se derivó una fórmula empírica para la atenuación del ruido de una fuente puntual que correspondió a:

dn = 2(Id1 – Idn)/6 . d1

donde dn corresponde a la distancia donde el sonido se atenua o potencia dependiendo de si es mayor o menor que la distancia de referencia d1. Además, se determinó que para una taladradora de construcción el dn sería de 3.225 m para una intensidad de 60 dBA (intensidad para la cual los ruidos comienzan a considerarse molestos) y de 130 dBA a una distancia de sólo 1 m (equivalente al sonido producido por un avión jet despegando y apenas 10 dBA por debajo del umbral de dolor que es 140 dbA). Sin embargo, para los ráster de atenuación de ruido nos interesa la Intensidad como función de la distancia y por tanto hay que logaritmar la ecuación anterior para obtener:

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

Para facilitar la obtención del ráster de atenuación, se utiliza otro cualquiera de la zona de referencia para obtener los parámetros de Geo transformación mediante GDAL. Esto sirve para determinar la distancia entre la fuente puntual y el centro de cada celda del ráster para el número de filas y columnas que deseemos. Las distancias se usan para determinar las intensidades correspondientes mediante la fórmula anterior y se almacenan en una matriz de valores que se convertirán en ráster mediante los métodos de GDAL. 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
from osgeo import gdal, osr
from math import log10
import numpy as np

layer = iface.activeLayer()
provider = layer.dataProvider()
filename = provider.dataSourceUri()

rows = layer.height()
cols = layer.width()
xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

transform = dataset.GetGeoTransform()

#position of point noise source
pt = QgsPoint(386517.303009, 4437489.61456)

#construction drill noise parameters
d1 = 10 #distance (meters)
Id1 = 110 #Intensity (decibels) 

values = [ [] for i in range(rows) ]

x = transform[0] + xsize/2
y = transform[3] - ysize/2

points = [] #for verifiying purposes

for i in range(rows):
    for j in range(cols):
        dn = pt.distance(QgsPoint(x,y))
        points.append((x,y))
        Idn = Id1 - (6./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_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

Después de ejecutado en la Python Console de QGIS se obtuvo lo siguiente:


Más allá de la zona difusa de color azul, los valores son menores de 60 dBA. No obstante, un mapa de visibilidad y un DEM de la zona de referencia permitirían una mayor precisión en los valores de Intensidad (aunque no se tomarían en cuenta los fenómenos de interferencia, constructiva o destructiva, que son susceptibles de sufrir las ondas de sonido al chocar con los accidentes geográficos).

No hay comentarios: