Tengo archivos ráster con los datos de precipitación (4x4 km), mensuales, por cerca de 100 años. Tengo un archivo de forma de la cuenca (que incluye a varios de los polígonos como sub-cuencas). Necesito la precipitación media de la cuenca para cada mes. No quiero hacer esto a mano en ArcMap porque probablemente tendrá que ser repetido un montón de veces (para otras cuencas, y cuando se dispone de más datos), así que estoy de codificación.
He escrito el código en R ahora, usando los paquetes raster y rgdal. Funciona muy bien en el sentido de que tiene una función para calcular el promedio ponderado por sub-cuenca; comprueba que la fracción de una celda ráster cae dentro del polígono y la utiliza para calcular un promedio ponderado para el polígono. Debido a que algunas de las sub-cuencas son pequeñas en comparación a la trama, esta es una gran característica.
Y además, porque puedo acceder a el shapefile de información, tales como el área de la sub-cuencas, puedo calcular un promedio ponderado para toda la cuenca.
Sin embargo - esto es lento. También lo dice en la descripción de la función extract(). Calcular el valor ponderado de 4-5 segundos, así que cuando me encontré este script hace 30 años el valor de los datos mensuales, se llevó a 37 minutos (hay un par de cálculos más además de la media ponderada).
Tengo la sensación de que podría ser más rápido, el uso de Python, pero no estoy seguro de cómo acercarse a este. Me hizo empezar con esto antes de que el código R (ver abajo), pero cambió mi atención a R. parecía un poco complicado, ya que ha sido un tiempo que me codificado en Python, y nunca he usado Python para el análisis geoespacial.
También, una prueba inicial parece mostrar que en Python, no hay ponderado promedio sucediendo, y que los pequeños polígonos que no contienen el centro de una celda ráster, se asignaron 0 o None.
Así que mis preguntas son:
- Si es posible en Python, se ponderado promedio de esa manera ser más rápido?
- Cómo acceder a otros datos asociados con el shapefile (tales como el área)?
Abajo está la muy muy básico código que tengo hasta ahora.
from rasterstats import zonal_stats
import gdal
basin="Steinhatchee_Project.shp"
rain="PRISM_ppt_stable_4kmM3_201504_bil.bil"
# NOTE From Python documentation: "There is no right or wrong way to rasterize a vector.
# The default strategy is to include all pixels along the line render path (for lines),
# or cells where the center point is within the polygon (for polygons)
stats=zonal_stats(basin,rain)
# get average for every polygon (there are 19)
data=([f['mean'] for f in stats])
# take out the "None" values...
data2=[x for x in data if x is not None]
# ... so I can calculate the mean for the whole watershed
mean(data2)
data ve de la siguiente manera. Nota: el None...
Out[76]:
[115.64999389648438,
83.07333374023438,
92.33000183105469,
73.35416666666667,
97.36499786376953,
92.03500366210938,
80.30999755859375,
102.8699951171875,
94.3499984741211,
110.99250030517578,
80.18000030517578,
67.6749979654948,
107.825,
82.63999938964844,
88.375,
59.769999186197914,
63.48249816894531,
None,
None]
Mi especificaciones del sistema: Windows 7, Python 2.7.9 (utilizando el Enthought Dosel IDE) de 32 bits, GDAL-1.11.3, rasterio-0.29.0, bien formada-1.5.13
