1

I was trying to mask the building given its TIFF file and GeoJSON file creating a new raster TIFF file. I don't know what is the problem, the result is all black, here is the code below:

from osgeo import gdal, ogr

def create_building_mask(rasterSrc, vectorSrc, npDistFileName='', 
                        noDataValue=0, burn_values=1):

    source_ds = ogr.Open(vectorSrc)
    source_layer = source_ds.GetLayer()


    srcRas_ds = gdal.Open(rasterSrc)
    cols = srcRas_ds.RasterXSize
    rows = srcRas_ds.RasterYSize


    memdrv = gdal.GetDriverByName('GTiff') 
    dst_ds = memdrv.Create(npDistFileName, cols, rows, 1, gdal.GDT_Byte)
    dst_ds.SetGeoTransform(srcRas_ds.GetGeoTransform())
    dst_ds.SetProjection(srcRas_ds.GetProjection())
    band = dst_ds.GetRasterBand(1)
    band.SetNoDataValue(noDataValue)    

    gdal.RasterizeLayer(dst_ds, [1], source_layer, burn_values=[burn_values])

    dst_ds = 0

    return 

create_building_mask("img.tif", "img.geojson", npDistFileName='mask.tif', noDataValue=0, burn_values=1)

The result of mask.tif is all black, instead of masked img. enter image description here

Improve this question
3
  • You could try gdal.RasterizeLayer(dst_ds, [1], source_layer, burn_values=[burn_values], options='-init 0') which should initialize the raster to 0 before burning in 1 for your buildings assuming that your source_layer has the same CRS as your template srcRas_ds.
    – Michael Stimson
    Commented Jan 9, 2020 at 1:43
  • You have created tiff that is of type Byte. If you use pixel value 1 for the masked area then it is pretty close to black for viewers which are showing the whole range 0-255 on screen. You can for example create mask image as 1 bit tiff, use pixel value 255 for a white mask, or do LUT stretch with your image viewer.
    – user30184
    Commented Jan 9, 2020 at 8:45
  • @user30184 I changed the' burn_bvalues' to 255 as you suggested, it worked. Thank you, your explanation made me clearly understand.
    – user8005429
    Commented Jan 9, 2020 at 16:06

1 Answer 1

Reset to default
2

You may be interested in rioxarray to do this type of operation. It uses rasterio which is an alternative python wrapper for GDAL.


import rioxarray
import json


# load in the geojson file
with open("img.geojson") as igj:
    data = json.load(igj)
# if GDAL 3+
crs = data["crs"]["properties"]["name"]
# crs = "EPSG:4326" # if GDAL 2
geoms = [feat["geometry"] for feat in data["features"]]

# create empty mask raster based on the input raster
rds = rioxarray.open_rasterio("img.tif").isel(band=0)
rds.values[:] = 1
rds.rio.write_nodata(0, inplace=True)

# clip the raster to the mask
clipped = rds.rio.clip(geoms, crs, drop=False)

enter image description here

# or clip it inverted if that is what you wanted
clipped = rds.rio.clip(geoms, crs, drop=False, invert=True)

enter image description here

# write output to file
clipped.rio.to_raster("mask.tif", dtype="uint8")

Alternatively, if you just want to mask out the buildings in the original tif file:

rds = rioxarray.open_rasterio("img.tif")
clipped = rds.rio.clip(geoms, crs, drop=False, invert=True)

enter image description here

Improve this answer
3
  • Thank you for your reply. It is really helpful. I used your code, and I ploted the output using " input_img = plt.imread('mask.tif') plt.imshow('input_img')", it showed masked result as expected, however if I open the "mask.tif" file itself, it is all black. And I noticed that the original "img.tif" has 'interleave': 'pixel' and the saved output file which is "mask.tif" has 'interleave': 'band', is that the reason why I open the 'mask.tif' file itself is all black? Thank you very much.
    – user8005429
    Commented Jan 9, 2020 at 4:27
  • You could pass in interleave="pixel" into the to_raster function. What are you using to read in the file that shows it as all black?
    – snowman2
    Commented Jan 9, 2020 at 5:04
  • I am using Preview to open the "mask.tif" file and it is all black, but I converted it to .png file. Problem solved. Many thanks .
    – user8005429
    Commented Jan 9, 2020 at 6:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.