Monday 30 September 2019

pyqgis - convert vector (point shape) to raster (tif) using python gdal lib in qgis


I'm trying to convert a vectorlayer (point- shapefile) to a raster-dataset (tiff) using python gdal-library in qgis.


On Graphical UserInterface I simple choose the menue "raster->conversion->rasterize". To get what I want, I fill in the file names in tool dialog box and choose "resolution in map units per pixel" horizonthal: 0.033 and vertical: 0.0164. That works fine, but when I try to script I orientated on documentation and adapted the code from source: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#convert-an-ogr-file-to-a-raster like this:


# 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])

It results in error message:



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'

It doesn't accept float values for x_res and y_res,
but how can I enter the proper geometry for my dataset?


Additional Info:



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



It's telling you it doesn't like your values as int, they round to 0 and you can't have a 0 row, 0 column raster. That's because you're supplying the cell size and not the rows and columns. When you create a dataset it needs to know how many rows and columns it is, you then tell it how big the cells are in the geotransform:


# 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])


I don't know what the spatial reference of your input is so I don't know if your X and Y cell sizes are appropriate. Yes, you can have different cell sizes for X and Y. If your cell size should be 1 then change the values to 1 or just make cols = int(x_max - x_min) and rows = int(y_max - y_min) as anything divided by 1 is unchanged... but the create raster still needs an int value.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...