Masking with rasterio in Python

Question

Introduction

I have a shapefile for my area of interest (which in this case is Tehran) and also a TIF file that contains the categorical information(cloud status). In this case, I want to mask the TIF file and clip my area of interest (Tehran) using the shapefile provided.

This TIF file just contains numbers between 0 and 5. Therefor if I run the below code:

with rasterio.open("SVDNB_npp_d20120301.vcld.tif") as src:
  print(src.read(1))

I will have an output like:

[[0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 [0 0 0 ... 0 0 0]
 ...
 [4 4 4 ... 4 4 4]
 [4 4 4 ... 4 4 4]
 [0 0 0 ... 0 0 0]]

I am using the below code to mask the TIF file(The few first lines are just to change coordinate references system) using the shapefiles provided:

import fiona
import rasterio
import rasterio.mask
from rasterio.windows import Window
from pyproj import Transformer
import geopandas as gpd

# To Change the CRS
df = gpd.read_file('tehran.shp').to_crs(4326)
df.to_file("tehran.shp")

# To mask the data
with fiona.open("tehran.shp", "r") as shapefile:
  shapes = [feature["geometry"] for feature in shapefile]
with rasterio.open("SVDNB_npp_d20120301.vcld.tif") as src:
  out_image, transformed = rasterio.mask.mask(src, shapes, crop=True, filled=True)
  out_meta = src.meta

Issue

The problem is that the masked variable, which in this case is represented by out_image variable, is supposed to contain just numbers between 0 to 5 but, unfortunately, whenever I call this variable name, it outputs is like:

array([[255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       ...,
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255],
       [255, 255, 255, ..., 255, 255, 255]], dtype=uint8)

I am confused and not sure why the numbers now are changed to color-bits.

Answer 1

Your code works fine. If you look at SVDNB_npp_d20120301.vcld.tif (src.meta) you'll see nodata = 255 and you specified filled=True which means fill out_image with nodata values. So when you print your array, you are seeing the nodata values around the masked area of actual values.

If you write out the image and load it in a GIS, you'll see:

# To mask the data
with fiona.open("tehran.shp", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

with rasterio.open("SVDNB_npp_d20120301.vcld.tif") as src:
    out_image, transformed = rasterio.mask.mask(src, shapes, crop=True, filled=True)
    out_profile = src.profile.copy()

out_profile.update({'width': out_image.shape[2],'height': out_image.shape[1], 'transform': transformed})
with rasterio.open("tehran.tif", 'w', **out_profile) as dst:
    dst.write(out_image)

If you just want a 1D array of 0-5, just filter out the nodata

out_image[out_image<=5]]

Answer 2

I am not sure what is the problem, but I recommend you, cut your raster using gdal.warp

from osgeo import gdal

if you use a geometry:

gdal.Warp(destNameOrDestDS = outputpath, # directory output
          srcDSOrSrcDSTab  = intputpath, # directoy intput
          cutlineDSName    = shapefile or geojson, # vector file
          cropToCutline    = True, # Select True
          copyMetadata     = True, # optional
          dstNodata        = 0)    # if you have values of zeros

and if you want to use a extent:

gdal.Warp(destNameOrDestDS = outputpath, 
          srcDSOrSrcDSTab  = intputpath,
          outputBounds     = (minx, miny, maxx, maxy),
          cropToCutline    = True,
          copyMetadata     = True)

Answer 3

The issue can also arise if the pixel size is larger than the shapefile, resulting in nodata value. simply adding 'all_touched=True' in mask function solves it