3

In a dictionary such as below, I have a series of x,y coordinates as key and distance measured as values in dictionary.

I want to convert this dictionary to a raster layer with PyQGIS, Can I do it?

dictionary example:

mergeDict = {(231111,3.709e+06): 0, (231611,3.7085e+06): 0, (231611,3.7088e+06): [129.4689990991847], (232111,3.7093e+06): [230.72989620660874], (231011,3.7095e+06): 0, (231711,3.7088e+06): [56.121848744267226], (231211,3.7086e+06): 0, (231611,3.709e+06): [126.10543581749376], (231111,3.7093e+06): 0, (231411,3.7096e+06): [88.44457301849167], (232311,3.7096e+06): 0, (232211,3.7087e+06): 0, (231411,3.7089e+06): [147.3714773298335], (231511,3.7093e+06): [84.42704048812357], (231311,3.7094e+06): [91.15471831360921], (231411,3.7087e+06): [126.00603268731915], (231411,3.7094e+06): [96.62363082268323], (231711,3.7091e+06): 0, (231811,3.7089e+06): [55.93565232902348],...}

I tried to create an empty raster then disaggregate x, y coordinates and values by below code but I don't know how to add this data to raster layer.

outDs = drive.Create("/Data/aaSamplePop1.tif", src_cols, src_rows, 1, gdal.GDT_Float32)
outBand = outDs.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

for i in mergeDict.keys():
    for k, j in itertools.combinations(i, 2):
        print k, j
  • Hi, have a look at this: pcjericks.github.io/py-gdalogr-cookbook/… – nickves Jul 10 '17 at 19:23
  • Hi @nickves, thanks for your attention, yes I saw it, my first problem is to extract x, y coordinate and values from dictionary. I can do it by two order, (1) "for i in mergeDict.keys(): for k, j in itertools.combinations(i, 2): print k, j" (2) "for v in mergeDict.values(): print v" while I need extract them in one order and then convert them raster. – nickan Jul 10 '17 at 19:38
3

As I said before related to dictionaries, they don't preserve order as lists and it is necessary to be careful to order data adequately for avoiding errors. However, it's also necessary that points are regularly spaced and form a complete set as in my next example code. In this case, you can get a memory point layer and rasterize it for its values extracted by using grass:v.to.rast.attribute method from processing tools. Complete code is:

import processing

mergeDict = {(355008,4.47339e+06):0, (355008,4.47332e+06):0, 
             (355008,4.47324e+06):[129.4689990991847], (355008,4.47317e+06):[230.72989620660874], 
             (355008,4.4731e+06):0,  (355008,4.47302e+06):[56.121848744267226], 
             (355082,4.47339e+06):0, (355082,4.47332e+06):[230.72989620660874], 
             (355082,4.47324e+06):0, (355082,4.47317e+06):[88.44457301849167], 
             (355082,4.4731e+06):0,  (355082,4.47302e+06):0, 
             (355156,4.47339e+06):[147.3714773298335], (355156,4.47332e+06):[84.42704048812357], 
             (355156,4.47324e+06):[91.15471831360921], (355156,4.47317e+06):[126.00603268731915], 
             (355156,4.4731e+06):[96.62363082268323],  (355156,4.47302e+06):0, 
             (355230,4.47339e+06):[55.93565232902348], (355230,4.47332e+06):0, 
             (355230,4.47324e+06):0, (355230,4.47317e+06):0, 
             (355230,4.4731e+06):[91.15471831360921],  (355230,4.47302e+06):0, 
             (355304,4.47339e+06):0, (355304,4.47332e+06):[230.72989620660874], 
             (355304,4.47324e+06):0, (355304,4.47317e+06):0, 
             (355304,4.4731e+06):0, (355304,4.47302e+06):[126.10543581749376]}

mergeList = [ [ point, mergeDict[point] ] for point in mergeDict ]

points = [ QgsPoint(item[0][0], item[0][1]) for item in mergeList ]
values = [ item[1] for item in mergeList ]

new_values = []

for item in values:
    if type(item) == list:
        new_values.append(item[0])
    else:
        new_values.append(item)

epsg = 32612

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer&field=value:real""&index=yes"

mem_layer = QgsVectorLayer(uri,
                          'points',
                          'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i, new_values[i]])
    feat.setGeometry(QgsGeometry.fromPoint(points[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

xmin, ymin, xmax, ymax = mem_layer.extent().toRectF().getCoords()

TYPE = 1
HSPACING = 73.9887
VSPACING = 73.9887

extent = str(xmin - HSPACING/2)+ ',' + str(xmax + HSPACING/2)+ ',' +str(ymin - VSPACING/2)+ ',' +str(ymax + VSPACING/2) 

path = processing.runalg('grass:v.to.rast.attribute', 
                         mem_layer, #input <ParameterVector>
                         0,         #use <ParameterSelection> 0 for attribute
                         'value',   #column <parameters from input>
                         extent,    #GRASS_REGION_PARAMETER <ParameterExtent>
                         HSPACING,  #GRASS_REGION_CELLSIZE_PARAMETER <ParameterNumber>
                         0,         #GRASS_SNAP_TOLERANCE_PARAMETER <ParameterNumber>
                         0,         #GRASS_MIN_AREA_PARAMETER <ParameterNumber>
                         None)      #output <OutputRaster>

raster = QgsRasterLayer(path['output'],
                        'raster_points')

QgsMapLayerRegistry.instance().addMapLayer(raster)

After running it at the Python Console I got:

enter image description here

where it can be also observed points that were rasterized in this complete process. I hope that it helps.

  • Thanks for your help. In my result there is not any raster layer for output in TOC of QGIS. – nickan Jul 11 '17 at 10:00
2

This is not a complete PyQGIS solution but it might help.

Start by writing the keys and values of the dictionary to a csv file. Since some of the values are lists they need to be "unlisted":

import csv

mergeDict = {(231111,3.709e+06): 0, (231611,3.7085e+06): 0, (231611,3.7088e+06): [129.4689990991847], (232111,3.7093e+06): [230.72989620660874], (231011,3.7095e+06): 0, (231711,3.7088e+06): [56.121848744267226], (231211,3.7086e+06): 0, (231611,3.709e+06): [126.10543581749376], (231111,3.7093e+06): 0, (231411,3.7096e+06): [88.44457301849167], (232311,3.7096e+06): 0, (232211,3.7087e+06): 0, (231411,3.7089e+06): [147.3714773298335], (231511,3.7093e+06): [84.42704048812357], (231311,3.7094e+06): [91.15471831360921], (231411,3.7087e+06): [126.00603268731915], (231411,3.7094e+06): [96.62363082268323], (231711,3.7091e+06): 0, (231811,3.7089e+06): [55.93565232902348]}

coordcsv=r'C:\folder\toRaster.csv'

with open(coordcsv,'w') as csvfile:
    coordwriter = csv.writer(csvfile, delimiter=',')
    coordwriter.writerow(['X','Y','VALUE'])
    for key, value in mergeDict.iteritems():
        x,y=key
        if type(value)==list:
            newvalue=value[0]
        else:
            newvalue=value
        coordwriter.writerow((x,y,newvalue))

This will give you a csv looking like this:

X,Y,VALUE
231611,3708800.0,129.4689990991847
231111,3709300.0,0
231611,3708500.0,0
231411,3709400.0,96.62363082268323

Convert this to a raster using "Raster > Conversion > Rasterize" menu with a "Comma Separated Value (.csv)"

  • Hi @BERA, Thanks for your answer, Your solution is very nice to solve my first problem in extract x, y, value from a dictionary, but my goal is create a raster from a dictionary directly in pyqgis. – nickan Jul 11 '17 at 8:23
  • 1
    Nice. Maybe someone else can show us how to convert the csv to raster using pyqgis – BERA Jul 11 '17 at 8:24

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.