6

I have two rasters, one has a 1 km resolution. The other has a spatial resolution of 30 m. This latter pixelwise is a classification. I want to overlay the 1 km raster on the classification raster and extract for each 1 km pixel the distribution of classes corresponding to the 30 m classification.

enter image description here

The 30 m classification includes 7 classes so I would want this process to produce 7 new layers at 1 km. One for each class. Each pixel of 1 km will now have associated a corresponding % of the presence of each class present in the 30 m classification. If one made a table where each row is a pixel then one would have a table like the following:

enter image description here

One additional difficulty is that the rasters are quite big. I would really like to achieve this using python.

EDIT:

I have added subsets of my two datasets a this link (at 30m and 1000m resolutions respectively), in case anyone wants to propose a proof of concept; would be amazing.

  • Not a full solution but just an idea: Maybe you could generate a fishnet grid of 1km x 1km resolution, and then use rectangles generated as inputs into the zonal statistics (or zonal statistics as table) tool. you could then use those statistics to calculate the % that you need, and then convert the rectangles to a raster. – BruceDoh Aug 14 '15 at 15:46
  • Can you post your data for me to work on? I think I have a solution, I want to test it on your data first though. I am thinking this would pretty simple in R but I can probably work it up in arcpy also... anyway, I look forward to trying to resolve this! – c0ba1t Aug 17 '15 at 16:31
  • How can I post data? I think it's quite big. Is it allowed to link to dropbox or something? I don't have arcgis, would it be really hard to implement it in python as is? Sorry for being so picky. – JEquihua Aug 17 '15 at 16:41
  • 1
    I apologize, I was under the impression you were using arc, no that is no worries, no apology necessary... as far as what is allowed -- see this post meta.stackexchange.com/questions/15821/… a dropbox link would be fine – c0ba1t Aug 17 '15 at 16:48
  • 1
    What do you mean by quite big ? If you haven't seen them, there is some useful slides from C. Garrard about handling raster in python with gdal. I guess these slides suggest you can iterate over each pixel of the 1*1km layer, then use this "window' for defining a block of pixels to read on the 30*30m layer and compute the % of each value in this block (just a guess). Anyway if you find the slides helpful, change the '4' to '5' and/or to '6' in the url for more slides about raster and python. – mgc Aug 17 '15 at 18:22
6
+50

Step 1

Make bit rasters for each of the unique classes. This can be a 1-band rasters for each class, or a single raster with a band for each class (e.g. GeoTIFF). If using GTiff, you can use the creation option NBITS=1 to conserve space. You may also want to consider twobit rasters to store three-valued logic where the third (e.g. 2) is NODATA, which would require NBITS=2.

For each band / class, the pixels will be 0 where that class does not exist, or 1 where that class exists.

import numpy as np
from osgeo import gdal
gdal.UseExceptions()

# Get data from raster with classifications
ds = gdal.Open('cropped_30m.tif')
band = ds.GetRasterBand(1)
class_ar = band.ReadAsArray()
gt = ds.GetGeoTransform()
pj = ds.GetProjection()
ds = band = None  # close

# Define the raster values for each class, to relate to each band
class_ids = (np.arange(31) + 1).tolist()

# Make a new bit rasters
drv = gdal.GetDriverByName('GTiff')
ds = drv.Create('bit_raster.tif', class_ar.shape[1], class_ar.shape[0],
                len(class_ids), gdal.GDT_Byte,
                ['NBITS=1', 'COMPRESS=LZW', 'INTERLEAVE=BAND'])
ds.SetGeoTransform(gt)
ds.SetProjection(pj)
for bidx in range(ds.RasterCount):
    band = ds.GetRasterBand(bidx + 1)
    # create boolean
    selection = (class_ar == class_ids[bidx])
    band.WriteArray(selection.astype('B'))
ds = band = None  # save, close

E.g. for Red = Band 10, Green = Band 13 and Blue = Band 27, it can be visualised:

30m

Black indicates it is another class / band.

Step 2

Use gdalwarp to get the average of 0 and 1 values from a coarser sample grid. This method requires GDAL >= 1.10.0. A really simply way to do this is to specify a target resolution on the command line, e.g. 1 km:

$ gdalwarp -ot Float32 -tr 1000 1000 -r average bit_raster.tif average_1km.tif

The result will have as many bands as classes, and each band will describe the fraction of that class for each pixel, between 0 and 1. If you prefer percent, then multiply it by 100 (e.g. see gdal_calc.py).

There is also a Python interface to gdalwarp (e.g. this answer), which can use a template raster for shape and where NODATA should be used. For example:

# Open raster from step 1
src_ds = gdal.Open('bit_raster.tif')

# Open a template or copy array, for dimensions and NODATA mask
cpy_ds = gdal.Open('cropped_1000m.tif')
band = cpy_ds.GetRasterBand(1)
cpy_mask = (band.ReadAsArray() == band.GetNoDataValue())

# Result raster, with same resolution and position as the copy raster
dst_ds = drv.Create('average2_1km.tif', cpy_ds.RasterXSize, cpy_ds.RasterYSize,
                    len(class_ids), gdal.GDT_Float32, ['INTERLEAVE=BAND'])
dst_ds.SetGeoTransform(cpy_ds.GetGeoTransform())
dst_ds.SetProjection(cpy_ds.GetProjection())

# Do the same as gdalwarp -r average; this might take a while to finish
gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Average)

# Convert all fractions to percent, and apply the same NODATA mask from the copy raster
NODATA = -9999
for bidx in range(dst_ds.RasterCount):
    band = dst_ds.GetRasterBand(bidx + 1)
    ar = band.ReadAsArray() * 100.0
    ar[cpy_mask] = NODATA
    band.WriteArray(ar)
    band.SetNoDataValue(NODATA)

# Save and close all rasters
src_ds = cpy_ds = dst_ds = band = None

E.g. the same three bands (classes) can be visualised as RGB:

1km

Or pull it into QGIS, and use the Identify features tool to inspect the fraction for each band/classification at a single location:

QGIS

These number add up to 100%, which is a good check.

  • It's running, ill let you know if it turned out fine, thank you. By the way gdal and numpy dtypes don't get along so well so I had to modify the code where you write the bands to disk: selection = (class_ar == class_ids[bidx]).astype("u1") band.WriteArray(selection) – JEquihua Aug 22 '15 at 19:31
  • It's going very very slow but getting there. One question, how does this process handle the projection, extent, etc of the output raster? What I'm worried about is how am I going to align the output raster to my original 1km raster? (30m and 1km are in have same projection). – JEquihua Aug 22 '15 at 19:43
  • @JEquihua updated to align to the 1km raster, and copy the NODATA locations, multiply the values to percent – Mike T Aug 23 '15 at 3:44

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.