2

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.

3
  • Does the TIFF file you're using have any georeferencing data in it? If it does, have you made sure it matches the CRS you're using for tehran.shp (epsg:4326)?
    – Felipe D.
    Jan 25 at 23:37
  • As per the help center please do not include chit chat like statements of appreciation within your posts.
    – PolyGeo
    Jan 25 at 23:39
  • @FelipeD., Yes. As you can see, after reading the shapefile into a Geopandas dataframe, I am using to_crs(4326) to match the coordinate references systems. Jan 26 at 7:04

2 Answers 2

6

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)

enter image description here

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

out_image[out_image<=5]]
4
  • This is very helpful. But still, there is a problem for me. I would like to access this array (which contains the number between 0-5) in python itself. How would I be able to get this array from rasterio? Thanks Jan 26 at 7:09
  • 1
    out_image is a numpy array.
    – user2856
    Jan 26 at 8:10
  • Yes. I am aware of that. But still, it contains 255. Do I have to use third-party software like QGIS to get the numbers scaled in the desired format(0-4)? Jan 26 at 8:19
  • There is a way to get the desired scaled numbers. I can re-read the tehran.tif file after writing it to a file. this would solve my problem entirely, but, of course, not efficiently. I would appreciate it if you could edit your answer in order to avoid writing and re-reading. Jan 26 at 9:01
1

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)

1
  • As per the help center please do not include chit chat like greetings in your posts.
    – PolyGeo
    Jan 26 at 0:27

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.