{"id":57,"date":"2018-07-04T02:34:56","date_gmt":"2018-07-04T02:34:56","guid":{"rendered":"http:\/\/wcwp.rf.gd\/2018\/07\/04\/https-landscapearchaeology-org-2018-numpy-loops\/"},"modified":"2018-07-04T02:34:56","modified_gmt":"2018-07-04T02:34:56","slug":"numpy-loops","status":"publish","type":"post","link":"https:\/\/landscapeanalysis.org\/?p=57","title":{"rendered":"Looping through numpy arrays (e.g. moving\/rolling window)"},"content":{"rendered":"

Numpy is the cornerstone of matrix based calculations in QGIS (and elsewhere). It does wonders with raster data (unless it hits the limit of available live memory\u2026).<\/p>\n

A recurrent problem with Numpy is the implementation of various looping routines, such as the sliding window which is frequently used in image filtering and other approaches focused on cell neighbourhood. Below is the illustration of the problem: for each cell the window needs to query a specified neighbourhood (square, circular or other).<\/p>\n

\"\"<\/p>\n

Now, there exists a solution in the form of \u201cstride tricks\u201d within the Numpy library. However it is really cryptic (the SciPy cookbook<\/a> uses the \u201cdevious\u201d epithet); after all the very name \u201ctricks\u201d implies it being \u2026 tricky. Numpy documentation actually discourages its use<\/a>. We can do it in a more understandable manner – which also means more configurable – and without any particular computational overhead.<\/p>\n

Let\u2019s first examine the problem. Without Numpy we would need four nested loops: two for traversing the matrix and two for the analysed window.<\/p>\n

\n
\n
#Level 1: traversing the matrix \n<\/span>for<\/span> y<\/span> in<\/span> range<\/span>(<\/span>rows<\/span>):<\/span>\n    for<\/span> x<\/span> in<\/span> range<\/span>(<\/span>columns<\/span>):<\/span>\n\n        #Level 2: traversing the window (3x3 size) \n<\/span>        for<\/span> win_y<\/span> in<\/span> range<\/span>(<\/span>y<\/span>-<\/span>1<\/span>,<\/span>y<\/span>+<\/span>2<\/span>):<\/span>\n            for<\/span>  win_x<\/span> in<\/span> range<\/span>(<\/span>x<\/span>-<\/span>1<\/span>,<\/span> x<\/span>+<\/span>2<\/span>):<\/span>\n                temp_sum<\/span> +=<\/span> matrix<\/span>[<\/span>win_y<\/span>,<\/span> win_x<\/span>]<\/span>\n\n       #Analysis inside the window is finished\n<\/span>       #Write result to a copy of the input matrix\n<\/span>       output<\/span>[<\/span>y<\/span>,<\/span> x<\/span>]<\/span> =<\/span> calculate_something<\/span>(<\/span>temp_sum<\/span>)<\/span> \n<\/code><\/pre>\n<\/div>\n<\/div>\n
We can immediately see that the level 2 loop can be easily replaced with a Numpy window: win<\/span> =<\/span> matrix<\/span>[<\/span>y<\/span>-<\/span>1<\/span>:<\/span> y<\/span>+<\/span>2<\/span>,<\/span> x<\/span>-<\/span>1<\/span>:<\/span> x<\/span>+<\/span>2<\/span>]<\/span><\/code>. Easy enough \u2026 but this is not the ideal approach. The window has 9 cells, while the matrix would have millions: shouldn\u2019t we focus on the level 1 loop instead?<\/p>\n
To avoid the level 1 loop we need to move the entire matrix, as when overlapping two sheets of paper. We\u2019re now working on offsets between the two – which is the only tricky part. The (0,0) cell of the window now corresponds to all cells in the matrix that are on the upper left relative to analysed cells – which simply means shifting the entire matrix to the upper left. In order to move the matrix in that direction, we need to calculate offsets in relation to the central point. That would be window_index<\/span> -<\/span> window_radius<\/span><\/code> , which for the upper left cell gives (-1, -1) (see below).<\/p>\n
\"\"<\/p>\n
The only problem resides in handling these offsets on the level of the entire matrix. We need to extract matrix views that do not spill over the borders AND that match each other in size and shape. The overlap issue is represented on the illustration below.<\/p>\n
\"\"<\/p>\n
This problem can be avoided with padding: adding extra row(s) and column(s) around the edges. The depth of the padded area should correspond to the maximum possible offset. The solution is now as follows:<\/p>\n
\n
\n
import<\/span> numpy<\/span> as<\/span> np<\/span>\n\n#a mock dataset\n<\/span>data<\/span>=<\/span>np<\/span>.<\/span>random<\/span>.<\/span>rand<\/span>(<\/span>5<\/span>,<\/span>5<\/span>)<\/span>\n\nrows<\/span>,<\/span> columns<\/span> =<\/span> data<\/span>.<\/span>shape<\/span>\ntemp_sum<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>((<\/span>rows<\/span>,<\/span> columns<\/span>))<\/span>\n\n# create a padded copy\n<\/span>pad<\/span> =<\/span> 1<\/span> \nmatrix<\/span> =<\/span> np<\/span>.<\/span>pad<\/span>(<\/span>data<\/span>,<\/span> pad<\/span>,<\/span> 'edge'<\/span>)<\/span>\n\n# Level 2: traversing the window (3x3 size)\n# You could (should) use numpy.ndenumerate() function, as well\n<\/span>for<\/span> y<\/span> in<\/span> range<\/span>(<\/span>3<\/span>):<\/span>\n    for<\/span> x<\/span> in<\/span> range<\/span>(<\/span>3<\/span>):<\/span>\n\n        # Level 1: handling the matrix \n<\/span>        # (rows, columns = data.shape !)\n<\/span>        temp_sum<\/span> +=<\/span> matrix<\/span>[<\/span>y<\/span> :<\/span> rows<\/span> +<\/span> y<\/span>,<\/span>\n                          x<\/span> :<\/span> columns<\/span> +<\/span> x<\/span>]<\/span>\n\n#Analysis inside the window is finished\n<\/span>output<\/span> =<\/span> calculate_something<\/span>(<\/span>temp_sum<\/span>)<\/span> \n<\/code><\/pre>\n<\/div>\n<\/div>\n
Note that the pad width should be half the window size: pad<\/span> =<\/span> window_size<\/span> \/\/<\/span> 2<\/span><\/code> where \/\/ is integer division. Window frame size should always be in odd numbers, otherwise it won\u2019t be placed symetrically.<\/p>\n
Refinement: unpadded solution<\/h3>\n
The above solution should be good enough for most of typical uses (for example, a small moving average filter). But, it\u2019s kind of ugly to use pads. We are producing copies of the original data, which can be particularly problematic for large offsets. Furthermore, padding may influence the result in edge zones in unexpected ways.<\/p>\n
Returning to the illustration above, we can see that the cropped window view cannot match the main view (the one pointing to the original and\/or temporary data). The latter also needs to be cropped, but in the opposite direction. This is logical: the top row cannot have a neighbour to the north anyway.<\/p>\n
Therefore, we need to adjust both views at the same time. Note that these views are opposite to each other, which means that we can play with swapping: when the window is placed to the north-west (-1, -1), the main view is to the south-east (1,1), but when the window is on the south-east (1,1), the main view will be on the north-west (-1, -1) (see the illustration above).<\/p>\n
To finish, let\u2019s throw in an additional parameter for steps – the range to be skipped between analysed cells (see numpy docs<\/a> for steps). Let\u2019s now wrap all this to a function which will return matching views:<\/p>\n
\n
\n
"""\nFunction returning two matching numpy views for moving window routines.\n- offset_y and offset_x refer to the shift in relation to the analysed (central) cell \n- size_y and size_x refer to the size of the data matrix (not of the window!)\n- view_in is the shifted view and view_out is the position of central cells\n\nAll comments are welcome at landscapeanalysis.org\/2018\/numpy-loops\n"""<\/span>\nimport<\/span> numpy<\/span> as<\/span> np<\/span>\n\ndef<\/span> view<\/span> (<\/span>offset_y<\/span>,<\/span> offset_x<\/span>,<\/span> size_y<\/span>,<\/span> size_x<\/span>,<\/span> step<\/span>=<\/span>1<\/span>):<\/span>\n \n    x<\/span> =<\/span> abs<\/span>(<\/span>offset_x<\/span>)<\/span>\n    y<\/span> =<\/span> abs<\/span>(<\/span>offset_y<\/span>)<\/span>\n \n    x_in<\/span> =<\/span> slice<\/span>(<\/span>x<\/span> ,<\/span> size_x<\/span>,<\/span> step<\/span>)<\/span> \n    x_out<\/span> =<\/span> slice<\/span>(<\/span>0<\/span>,<\/span> size_x<\/span> -<\/span> x<\/span>,<\/span> step<\/span>)<\/span>\n \n    y_in<\/span> =<\/span> slice<\/span>(<\/span>y<\/span>,<\/span> size_y<\/span>,<\/span> step<\/span>)<\/span>\n    y_out<\/span> =<\/span> slice<\/span>(<\/span>0<\/span>,<\/span> size_y<\/span> -<\/span> y<\/span>,<\/span> step<\/span>)<\/span>\n \n    # the swapping trick    \n<\/span>    if<\/span> offset_x<\/span> <<\/span> 0<\/span>:<\/span> x_in<\/span>,<\/span> x_out<\/span> =<\/span> x_out<\/span>,<\/span> x_in<\/span>                                 \n    if<\/span> offset_y<\/span> <<\/span> 0<\/span>:<\/span> y_in<\/span>,<\/span> y_out<\/span> =<\/span> y_out<\/span>,<\/span> y_in<\/span>\n \n    # return window view (in) and main view (out)\n<\/span>    return<\/span> np<\/span>.<\/span>s_<\/span>[<\/span>y_in<\/span>,<\/span> x_in<\/span>],<\/span> np<\/span>.<\/span>s_<\/span>[<\/span>y_out<\/span>,<\/span> x_out<\/span>]<\/span>\n \n#a mock dataset\n<\/span>data<\/span>=<\/span>np<\/span>.<\/span>random<\/span>.<\/span>rand<\/span>(<\/span>5<\/span>,<\/span>5<\/span>)<\/span>\nwindow<\/span> =<\/span> 3<\/span>\n\nradius<\/span> =<\/span> int<\/span>(<\/span>window<\/span>\/<\/span>2<\/span>)<\/span>\n \nrows<\/span>,<\/span> columns<\/span> =<\/span> data<\/span>.<\/span>shape<\/span>\ntemp_sum<\/span> =<\/span> np<\/span>.<\/span>zeros<\/span>((<\/span>rows<\/span>,<\/span> columns<\/span>))<\/span>\n \n# Our window loop  \n<\/span>for<\/span> y<\/span> in<\/span> range<\/span>(<\/span>window<\/span>):<\/span>\n\n    #we need offsets from centre !\n<\/span>    y_off<\/span> =<\/span> y<\/span> -<\/span> radius<\/span>\n\n    for<\/span> x<\/span> in<\/span> range<\/span>(<\/span>window<\/span>):<\/span>\n        x_off<\/span> =<\/span> x<\/span> -<\/span> radius<\/span>\n \n        view_in<\/span>,<\/span> view_out<\/span> =<\/span> view<\/span>(<\/span>y_off<\/span>,<\/span> x_off<\/span>,<\/span> rows<\/span>,<\/span> columns<\/span>)<\/span>\n \n        temp_sum<\/span>[<\/span>view_out<\/span>]<\/span> +=<\/span> data<\/span>[<\/span>view_in<\/span>]<\/span>\n\n#Analysis inside the window is finished\n<\/span>output<\/span> =<\/span> calculate_something<\/span>(<\/span>temp_sum<\/span>)<\/span> \n<\/code><\/pre>\n<\/div>\n<\/div>\n
Perhaps we could do better, but the advantage of this approach is that we can adapt it to any imaginable shape of the moving window – because we conserve the loop over window cells. We are also working on views, as opposed to hard copies of the data, which is what Numpy likes a lot.<\/p>\n
Related posts<\/h2>\n
You can find an implementation of sliding window logic in a script for GIS raster analysis (Topographic position index)<\/a>.<\/p>\n
Any suggestions are welcome!<\/p>\n","protected":false},"excerpt":{"rendered":"
Numpy is the cornerstone of matrix based calculations in QGIS (and elsewhere). It does wonders with raster data (unless it hits the limit of available live memory\u2026). A recurrent problem with Numpy is the implementation of various looping routines, such as the sliding window which is frequently used in image filtering and other approaches focused […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[8],"class_list":["post-57","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-code"],"_links":{"self":[{"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=\/wp\/v2\/posts\/57","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=57"}],"version-history":[{"count":0,"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=\/wp\/v2\/posts\/57\/revisions"}],"wp:attachment":[{"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=57"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=57"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/landscapeanalysis.org\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=57"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}