convertir vector (punto forma) a raster (tif) usando python gdal lib en qgis

Question

Estoy intentando convertir una capa vectorial (shapefile de puntos) en un conjunto de datos raster (tiff) utilizando la biblioteca gdal de python en qgis.

En la interfaz gráfica de usuario simplemente elijo el menú "raster->conversión->rasterizar". Para obtener lo que quiero, relleno los nombres de los archivos en el cuadro de diálogo de la herramienta y elijo "resolución en unidades de mapa por píxel" Horizontal: 0,033 y vertical: 0,0164. Eso funciona bien, pero cuando intento hacer el script me oriento en la documentación y adapto el código de la fuente: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#convert-an-ogr-file-to-a-raster así:

# 1. Define pixel_size and NoData value of new raster
NoData_value = -9999
x_res = 0.03333378
y_res = 0.01666641
pixel_size = 1

# 2. Filenames for in- and output
_in = r"C:/Users/.../hoppla.shp"
_out = r"C:/Users/.../hoppla.tif"

# 3. Open Shapefile
source_ds = ogr.Open(_in)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 4. Create Target - TIFF
_raster = gdal.GetDriverByName('GTiff').Create(_out, x_res, y_res, 1, gdal.GDT_Byte)
_raster.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
_band = _raster.GetRasterBand(1)
_band.SetNoDataValue(NoData_value)

# 5. Rasterize
gdal.RasterizeLayer(_raster, [1], source_layer, burn_values=[0])

El resultado es un mensaje de error:

Traceback (most recent call last):
  File "", line 1, in 
  File "C:/Users/.../test.py", line 73, in 
    \_raster = gdal.GetDriverByName('GTiff').Create(\_out, x\_res, y\_res, 1, gdal.GDT\_Byte)
  File "C:\\PROGRA~1\\QGISBR~1\\apps\\Python27\\lib\\site-packages\\osgeo\\gdal.py", line 388, in Create
    return \_gdal.Driver\_Create(self, \*args, \*\*kwargs)
TypeError: in method 'Driver\_Create', **argument 3 of type 'int'**

No acepta valores flotantes para x_res y y_res,
pero ¿cómo puedo introducir la geometría adecuada para mi conjunto de datos?

Información adicional:

QGIS 2.6.1 Brighton on Win7
crs: WGS84 - epsg4326
input shape horizonthal point distance: 0.03333378 units
input shape vertical point distance: 0.01666641 units

Answer 1

Te está diciendo que no le gustan tus valores como int, que redondean a 0 y que no puedes tener una trama de 0 filas y 0 columnas. Esto se debe a que está proporcionando el tamaño de la celda y no las filas y columnas. Cuando creas un conjunto de datos necesita saber cuántas filas y columnas tiene, entonces le dices el tamaño de las celdas en la geotransformación:

# 1. Define pixel_size and NoData value of new raster
NoData_value = -9999
x_res = 0.03333378 # assuming these are the cell sizes
y_res = 0.01666641 # change as appropriate
pixel_size = 1

# 2. Filenames for in- and output
_in = r"C:/Users/.../hoppla.shp"
_out = r"C:/Users/.../hoppla.tif"

# 3. Open Shapefile
source_ds = ogr.Open(_in)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 4. Create Target - TIFF
cols = int( (x_max - x_min) / x_res )
rows = int( (y_max - y_min) / y_res )

_raster = gdal.GetDriverByName('GTiff').Create(_out, cols, rows, 1, gdal.GDT_Byte)
_raster.SetGeoTransform((x_min, x_res, 0, y_max, 0, -y_res))
_band = _raster.GetRasterBand(1)
_band.SetNoDataValue(NoData_value)

# 5. Rasterize why is the burn value 0... isn't that the same as the background?
gdal.RasterizeLayer(_raster, [1], source_layer, burn_values=[0])

No sé cuál es la referencia espacial de tu entrada, así que no sé si los tamaños de tus celdas X e Y son apropiados. Sí, puede tener diferentes tamaños de celda para X e Y. Si su tamaño de celda debe ser 1 entonces cambie los valores a 1 o simplemente haga cols = int(x_max - x_min) y rows = int(y_max - y_min) ya que cualquier cosa dividida por 1 no cambia... pero la trama de creación sigue necesitando un valor int.