I am using GDAL to process a lot of rasters. I then need to combine all the arrays to create a single output raster that represents all the unique combinations of values. I.E the same thing as r.cross in GRASS or Combine_sa in arcgis.

I'd like to program a similar thing using the GDAL library. What is the algorithm? How would I handle rasters that are too large to fit in RAM? I.e., when chunking them into smaller arrays, how do I ensure that a combination in the nth chunk is given the same value if it appeared in a different chunk?

I'm struggling to come up with anything that will perform well.

  • When you say "a lot of rasters", what kind of number are we looking at? Also, how many unique numbers are in each raster? Doesn't have to be precise, just a rough estimate. – Mikkel Lydholm Rasmussen Apr 23 '15 at 16:18
  • @MikkelLydholmRasmussen There will be between 12-20 input rasters depending on configuration. For each raster, before combining a number of additional values (for each cell) are calculated, meaning that for a single cell, there will be about 30-40 values to combine. Unique numbers in each raster (to begin with) is not that many - usually 6. But with 30 - 40 values for a single cell, that is a lot of possible combinations. I have started programming a possible solution based on a btree to look up values, which I believe r.cross uses. – James Apr 24 '15 at 6:21
  • I think that you may run into problems with there simply being too many possible combinations for a raster to handle. Assuming 20 rasters, each with 6 unique values, that ends up as 6^20 = 3.6561584e+15 combinations, which is to say that you need to operate with 64-bit integer rasters, which are uncommon. Furthermore, the reverse look-up is going to be a pain to work through. – Mikkel Lydholm Rasmussen Apr 24 '15 at 6:54

As stated by @Mikkel you will be in trouble because of the large number of layers (r.cross is limited to 10 layers on purpose). Here is the equation if you are able to reduce the number of layers or if you don't have another solution:

\sum_i {(n^i) * c_i }

where i is the index of the image (starting at 0), n is the number of classes, and c_i is the class value for image i (from 0 to n-1). The result will not depend on the chunking.

As a remark, you will need something like the BIGTIFF format to handle this. You should also use other gdal creation option like COMPRESS and TILE. Finally I suggest you to try OTB for this king of process.

| improve this answer | |
  • @radoxju, is your equation meant for gdal's gdal_calc? – user1186 May 25 '16 at 5:01
  • 1
    @user1186 This is only the algorithm, not something you could use directly. Note that I updated it with a better version. in gdal_calc it would look like this --calc="A+nB+Cn^2+Dn^3+En^4+...+J*n^9" – radouxju May 25 '16 at 13:30

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.