1

I have two raster images, mask and data. The first one is a single-band tif with integer values that refer to classes. The second one, data, is a 3-band (RGB) tif. In the raster containing the classes, I assigned the value 100 to pixels that could not be identified.

Now I tried to set the R = G = B = 0 in data for every pixel where mask has the value 100.

To this end, I wrote the following code:

(For easier understanding of variable names: a = array, b = band, m = mask)

import numpy as np
from osgeo import gdal
mask = "C:/Users/Manuel/Nextcloud/dat/omk/03_2021/ortho_B3_4_mask.tif"
data = "C:/Users/Manuel/Nextcloud/dat/omk/03_2021/ortho_B3_4_CROP.tif"
NoDataVal = 100

from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
def blacken(data, mask, NoDataVal, outFile = "/vsimem/Blackened.tif"):
    dat = gdal.Open(data, GA_ReadOnly)
    msk = gdal.Open(mask, GA_ReadOnly)
    bm = msk.GetRasterBand(1)
    am = BandReadAsArray(bm)
    null = np.where(am == NoDataVal)
    am = None
    bm = None
    msk = None
    drv = gdal.GetDriverByName("GTiff")
    new = drv.Create(outFile, dat.RasterXSize, dat.RasterYSize,
                     dat.RasterCount, dat.GetRasterBand(1).DataType)
    CopyDatasetInfo(dat, new)
    for b in range(dat.RasterCount):
        band = dat.GetRasterBand(b+1)
        a = BandReadAsArray(band)
        a[null] = 0
        new.GetRasterBand(b+1).WriteArray(a)
        a = None
    new.FlushCache()
    dat = None
    return new

However, the values that are set to 0 are not the ones I asked for. It looks like the removed pattern is upside-down. Maybe this is due to the images location in the southern hemisphere? Or some issue in my code?

I attached an imageenter image description here that compares the 3 layers: The original data, the mask layer where the value 100 is the large gray matrix surrounding the patches, and the output image, where the black area is supposed to resemble the gray matrix of the mask.

Obviously, the patches don't fit while everything else seems to line up. Looking at the larger pattern, the black area appears to be vertically mirrored. The images were previously cropped to the same size and they are in the same reference system. I thought I treated them in the same way, so I wonder why the blackened pixels don't match the pixels where mask = 100. Perhaps there is an annoyingly simple mistake somewhere...

2
  • 1
    Test if CopyDatasetInfo(dat, new) was successful with copying the geotransform,
    – user30184
    May 20 '21 at 19:11
  • 1
    Thanks, the word "geotransform" was a good hint. I now used GetGeoTransform() on every variable I could find and eventually found the problem. May 21 '21 at 8:34
1

Your Copy loop does not cover all the bands,

range((dat.RasterCount-1)) will give you bands from [0, 1] for RGB

instead use range(dat.RasterCount)

1
  • Thanks, I already wondered why the colours are more greenish in the output (I added the -1 when trying to adjust the way GDAL and Python count differently (GDAL labels the first layer with 1). However, the black areas are still in the same place after I fixed the missing band... May 20 '21 at 16:15
1

I think this was the error:

In [30]: msk.GetGeoTransform()
Out[30]: 
(637791.2537999959,
 0.011203483858671787,
 0.0,
 7364288.8728,
 0.0,
 0.011203483858671787)

In [31]: new.GetGeoTransform()
Out[31]: 
(637791.2537999959,
 0.011203483858671787,
 0.0,
 7364406.1047,
 0.0,
 -0.011203354357799447)

I probably could have found this problem. But I didn't (partly because I'm not yet familiar with what to look for and how to do it in Python).

The mask was generated from a Shapefile via gdal.RasterizeLayer while the data served as a template. When I created the mask, I had not heared about the CopyDatasetInfo function and consequently, the mask was written "backwards". Fixing this, the output is as expected. Sorry for posting the part of the code that actually was not the problem.

Anyways, thanks for your time.

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.