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. – HMadadi 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.

| improve this answer | |
  • Thanks for your help. In my result there is not any raster layer for output in TOC of QGIS. – HMadadi Jul 11 '17 at 10:00
  • This is incomplete and not using raster data but vector points anyways. – bugmenot123 17 hours ago

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.