2

I want to create an azimuth raster layer base of single point shapefile, for this I write below code that create a memory point layer that each point in memory layer is center of each pixel and have a field by name of "Azimuth", in qgis when run this program it's ok, but when I want to load created raster layer "Azimuth_lyr.tif", display this message:

enter image description here

What is error(s) in this code?

my code:

def azimuth_raster(self):
    # Chosen Point layer with GDAL
    pnt_layer = self.dlg.ui.comboBox_sourcelayer.currentIndex()
    pntId = self.dlg.ui.comboBox_sourcelayer.itemData(pnt_layer)
    pntLayer = QgsMapLayerRegistry.instance().mapLayer(pntId)
    feat = pntLayer.getFeatures().next()
    pntxy = feat.geometry().asPoint()

    def pixel2coord(x, y):
        """Returns global coordinates from pixel x, y coords"""
        xp = (g_pixelWidth * x) + g_originX + (g_pixelWidth/2)
        yp = (g_pixelHeight * y) + g_originY + (g_pixelHeight /2)
        return(xp, yp)

    # Open tif file
    ds = QgsRasterLayer("C:\Temp_NoiseOnePoint\prox_Distance.tif","prox_Distance")
    global g_pixelWidth, g_pixelHeight
    g_pixelWidth = ds.rasterUnitsPerPixelX()
    g_pixelHeight = ds.rasterUnitsPerPixelY()

    # extent of the layer
    ext = ds.extent()
    global g_originX ,g_originY
    g_originX ,g_originY = (ext.xMinimum(),ext.yMinimum())
    cols = ds.width()
    rows = ds.height()

    #Create temporary vector layer and add to map
    vl = QgsVectorLayer("Point?crs=", "temporary_points", "memory")
    vl.dataProvider().addAttributes([ QgsField("X_coord", QVariant.Int),
                                      QgsField("Y_coord", QVariant.Int),
                                      QgsField("Azimuth", QVariant.Int) ])
    fet = QgsFeature()
    for row in range(0,rows):
        for col in range(0,cols):
            rspnt = pixel2coord(col,row)
            fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(rspnt[0], rspnt[1])))
            fet.initAttributes(3)
            fet.setAttribute(0, rspnt[0] )
            fet.setAttribute(1, rspnt[1])
            az = pntxy.azimuth(QgsPoint(rspnt[0], rspnt[1]))
            if az < 0.0:
                azim = 360 + az
            else:
                azim = az
            fet.setAttribute(2, azim)
            vl.dataProvider().addFeatures( [fet] )
            vl.updateExtents()
            vl.updateFields()

    QgsMapLayerRegistry.instance().addMapLayers([vl])

    # Convert azimuth point layer to azimuth raster layer

    mem_layer_provider = vl.dataProvider()

    srs = osr.SpatialReference()
    srs.ImportFromWkt(vl.crs().toWkt())
    out_rst = gdal.GetDriverByName('GTiff').Create("C:\Temp_NoiseOnePoint\Azimuth_lyr.tif",
                                                   cols, rows, 1, gdal.GDT_Byte)
    out_rst.SetGeoTransform(g_x_min, g_pixelWidth, 0, g_y_max, 0, -g_pixelHeight)
    out_rst.SetProjection(srs.ExportToWkt())
    band_rst = out_rst.GetRasterBand(1)
    band_rst.SetNoDataValue(g_NoData_value)
    band_rst.FlushCache()

    gdal.RasterizeLayer(out_rst, [1], mem_layer_provider, options = ["ATTRIBUTE=Azimuth"])

    # close bands
    band_rst = None

    # close rasters
    out_rst = None
  • Might not solve your problem but when it comes to setting paths, I always use double slashes for the final directory to where the file is stored (e.g. ds = QgsRasterLayer("C:/Temp_NoiseOnePoint//prox_Distance.tif","prox_Distance") and "C:/Temp_NoiseOnePoint//Azimuth_lyr.tif") – Joseph Jan 27 '16 at 10:32
  • Hi Joseph, thanks, but not solved. – nickan Jan 27 '16 at 11:17
2

The problem probably could be here:

.
.
.
    out_rst.SetGeoTransform(g_x_min, g_pixelWidth, 0, g_y_max, 0, -g_pixelHeight)
.
.
.

because 'SetGeoTransform' method needs a list as parameter.

Try out this:

.
.
.
    out_rst.SetGeoTransform( [ g_x_min, g_pixelWidth, 0, g_y_max, 0, -g_pixelHeight ] )
.
.
.

I hope that it helps

Editing Note:

I tried to solve this 'gdal.RasterizeLayer' issue unsuccessfully. So, I used another approach: the gdal_rasterize command. The code was:

import os

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

extent = layers[0].extent() #vector layer

xmin = extent.xMinimum()
ymax = extent.yMaximum()
xmax = extent.xMaximum()
ymin = extent.yMinimum()

xsize = 73.9887232281
ysize = 73.9887232281

vlayer = layers[0]  #vector layer
vprovider = vlayer.dataProvider()
uri = vprovider.dataSourceUri()

uri = uri.split("|")

(root, filename) = os.path.split(uri[0])

filename = filename[:-4]

#extent of raster
extent = "-te " + str(xmin) + " " + str(ymin) + " " + str(xmax) + " " + str(ymax) + " "

#resolution of raster
resolution = " -tr " + str(xsize) + " " +  str(ysize) + " "

attribute_name = " -a id "

input_vector_line = " -l " + filename + " " + uri[0]
output_raster = " /home/zeito/pyqgis_data/random_points.tif"

cmd = "gdal_rasterize -at " + \
                     extent + \
                 resolution + \
          input_vector_line + \
             attribute_name + \
              output_raster

print cmd

os.system(cmd)

I tested the script with this shapefile (it will be rasterized by id field):

enter image description here

After running the script the result was successfully. It could be corroborated at the next image where I 'zoom in' on point with id = 23.

enter image description here

  • Thanks, It's good suggestion and work but in raster layer "Azimuth_lyr.tif" created all pixels have zero value no "Azimuth" values from memory layer!! why? – nickan Jan 27 '16 at 13:22
  • 1
    I think that the problem is now here: gdal.RasterizeLayer(out_rst, [1], mem_layer_provider, options = ["ATTRIBUTE=Azimuth"]) – xunilk Jan 27 '16 at 17:25
  • That's right, problem is in "mem_layer_provider", I don't know how get memory layer, I think "mem_layer_provider = vl.dataProvider() " is not correct method to get memory layer. – nickan Jan 27 '16 at 17:36
  • Your memory layer is vl. The mem_layer_provider is an object of QgsVectorDataProvider. – xunilk Jan 27 '16 at 17:54
  • 1
    I used another approach, the 'gdal_rasterize' command, and it worked. – xunilk Jan 28 '16 at 11:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.