2

I have a 0.5 X 0.5 metre resolution raster. I would like to extract the maximum value points from the raster for every 2 metre square cell. And the extracted points should be in the centre of the source pixel.

enter image description here

In the shown example, the vertex point 18 should be extracted which is the maximum value within the 2m x 2m. And the position of the vertex point 18 should not be altered. I tried this using resampling tool and then extracting the vertices, but this puts the extracted value in the middle of the 2m x 2m cell.

  • Make vector polygon grid 2x2 m and use grass statistic plugin. Check min, max options. – nagib Nov 13 '19 at 15:19
2

This can be done with 'ReadAsArray' method of Python GDAL module by using a 4x4 "window" and storing maximum points for each window in a memory layer. For testing my approach, I created a raster of resolution 0.5x0.5, with 20x20 (rows,columns) and values between 1 and 1000 for avoiding repeated values.

My code is as follows; where print lines were used for corroborating that code works as expected (so, them could be eliminated).

from osgeo import gdal, osr
import numpy as np

driver = gdal.GetDriverByName('GTiff')
filename = "/home/zeito/pyqgis_data/test_0_5_0_5_clip.tif"
dataset = gdal.Open(filename)
band = dataset.GetRasterBand(1)

transform = dataset.GetGeoTransform()

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]

cols = dataset.RasterXSize
rows = dataset.RasterYSize

nb = 4

idxs = []

for i in range(0, rows, nb):
    for j in range(0, cols, nb):
        data = band.ReadAsArray(j, i, nb, nb)
        print (data)
        max = np.max(data)
        print ("max: ", max)
        print ("indices maximum original matrix")

        try:
            for m, item in enumerate(data):
                for n, element in enumerate(item):
                    if element == max:
                        idxs.append([m + i, n + j])
                        print(m + i, n + j)

        except TypeError:
            print ("there is no matrix")

x0 = xOrigin + pixelWidth/2
y0 = yOrigin - pixelHeight/2

points = []

for item in idxs:
    x = x0 + item[1]*pixelWidth
    y = y0 - item[0]* pixelHeight
    points.append(QgsPointXY(x,y))

wkt = dataset.GetProjection()

crs = QgsCoordinateReferenceSystem()
crs.createFromWkt(wkt)

epsg = crs.postgisSrid()

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

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

prov = mem_layer.dataProvider()

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

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

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

After running above code, I got following result printed in Python Console only for corroboration purpose.

[[557 231 752 950]
 [937   3 282 817]
 [171 548 846 977]
 [929 788 445  64]]
max:  977
indices maximum original matrix
2 3
[[240 503   7 417]
 [ 84 691 877   7]
 [762 163 243 556]
 [644 822 484 860]]
max:  877
indices maximum original matrix
1 6
[[ 38 124 342 813]
 [435 548 524  62]
 [286 101 439 690]
 [598  64   6 589]]
max:  813
indices maximum original matrix
0 11
[[773 594 866 548]
 [427 921 883 621]
 [730 262 270 757]
 [464 158 197 969]]
max:  969
indices maximum original matrix
3 15
[[786 997 865 821]
 [166 118 365 254]
 [483 398 212 419]
 [116 349 555 835]]
max:  997
indices maximum original matrix
0 17
[[640 608 836 124]
 [583 300 605 725]
 [888 927 848 377]
 [866 858 101 396]]
max:  927
indices maximum original matrix
6 1
[[161 301 214 915]
 [760 707 844 739]
 [325 645 883 899]
 [120 307 909 711]]
max:  915
indices maximum original matrix
4 7
[[376 232 369 914]
 [454 300 472 453]
 [970 350 651 959]
 [920 102 706 305]]
max:  970
indices maximum original matrix
6 8
[[ 83  45 620 414]
 [553 493 856 315]
 [572 532 872 384]
 [162 176  91 839]]
max:  872
indices maximum original matrix
6 14
[[222 100   6 767]
 [626 286 227 360]
 [ 58 138 392 251]
 [295 870 719 885]]
max:  885
indices maximum original matrix
7 19
[[168 751 925 455]
 [950  27 809 331]
 [290 800 581 509]
 [ 92 400 144 934]]
max:  950
indices maximum original matrix
9 0
[[891 305 267 809]
 [ 19 457 455 554]
 [303 727 537 865]
 [ 66 253  40 740]]
max:  891
indices maximum original matrix
8 4
[[ 67 831 329 703]
 [804 990 563 263]
 [384 761 801 629]
 [849 441 620 479]]
max:  990
indices maximum original matrix
9 9
[[980 912 581 658]
 [375 879 202 317]
 [ 11 691 138 782]
 [548 538 145 217]]
max:  980
indices maximum original matrix
8 12
[[512 486 145 897]
 [ 79 549   7 618]
 [401 378 826 577]
 [944  35 716 875]]
max:  944
indices maximum original matrix
11 16
[[334 815 132  13]
 [103 767 803 610]
 [532 493 457 809]
 [420  11 311 198]]
max:  815
indices maximum original matrix
12 1
[[201 471 685 237]
 [670   8 707 373]
 [  8 580 690 447]
 [630 869   1   7]]
max:  869
indices maximum original matrix
15 5
[[108 490 358 911]
 [669 285 108 640]
 [148 843 546 293]
 [701 563  79 870]]
max:  911
indices maximum original matrix
12 11
[[339 768 761 774]
 [423  30 326 145]
 [ 22  32 270 764]
 [832 512 672 950]]
max:  950
indices maximum original matrix
15 15
[[416 194 733 435]
 [937 868 344 236]
 [343 852 713 597]
 [163 536 378 610]]
max:  937
indices maximum original matrix
13 16
[[ 64 780 655 329]
 [710 583  48 390]
 [256  73 133  96]
 [446 139 587 703]]
max:  780
indices maximum original matrix
16 1
[[522 609 373 742]
 [117 780 443   2]
 [151 755 598  50]
 [891 365 398 589]]
max:  891
indices maximum original matrix
19 4
[[986 237 501 787]
 [237 586 376 773]
 [187 597 884 789]
 [944 770 441 730]]
max:  986
indices maximum original matrix
16 8
[[915 868  32 780]
 [156 528 133  86]
 [707 945 425 766]
 [190   5 142 960]]
max:  960
indices maximum original matrix
19 15
[[942 404 576 803]
 [282 122 496 181]
 [892 781 924  54]
 [ 30 149 168 742]]
max:  942
indices maximum original matrix
16 16

In following image can be observed points situated in the center of each pixel where maximum values are (also corroborated with Value Tool QGIS plugin). It was added a red grid (4x4; 2m x 2m) to facilitate the visualization.

enter image description here

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.