2

I want to get extent of raster dataset after excluding the extents containing no-data values in the raster data which I then want to export as a shapefile. Till now I have used the following code:

import rasterio

with rasterio.open(r'C:\Users\Desktop\PRODUCT1\Landsat\T1.tif') as src: bounds = src.bounds

from shapely.geometry import box geom = box(*bounds)

import geopandas as gpd df = gpd.GeoDataFrame({"id":1,"geometry":[geom]}) df.crs = src.crs df.to_file(r"C:\Users\Desktop\PRODUCT1\Landsat\T1_boundary.shp")

After running the above code I am getting raster extents containing no-data values.

Can someone please help me out in getting the raster extents excluding the extents having no-data values.

RRSC NGP
  • 658
  • 1
  • 8
  • 19

1 Answers1

3

If you're looking for a polygon of the actual valid data, you can use rasterio.features.shapes() to polygonize a boolean mask of your data as the source, e.g. (data !== src.nodata).astype(np.uint8).

Example:

import numpy as np

import rasterio from rasterio import features from shapely.geometry import shape

path = '/vsicurl/http://download.osgeo.org/geotiff/samples/usgs/m30dem.tif'

with rasterio.open(path) as f: image = f.read(1) # use f.nodata if possible; it's not defined on this particular image nodata = -32768 # create a binary image, 0 where there's nodata, 1 where it's valid is_valid = (image != nodata).astype(np.uint8) # vectorize the binary image, supplying the transform so it returns maps coords for coords, value in features.shapes(is_valid, transform=f.transform): # ignore polygons corresponding to nodata if value != 0: # convert geojson to shapely geometry geom = shape(coords) print(geom)

If you're truly looking for the extent of valid data, create that same array, feed it into np.where() to get array indices/image coordinates where valid data exists, take the min/max of the axes, and use the src.transform.xy() method to convert back to map coords

mikewatt
  • 5,083
  • 9
  • 21
  • I tried your suggestion. Inside "with rasterio.open('....')", I had used following commands: " image = src.read(1) " & then "image = (image != src.nodata).astype(np.uint8)". Post this when I am using "plt.imshow(image)" then the plot contains the whole image including "No Data" region. – RRSC NGP Aug 22 '20 at 14:49
  • @RRSCNGP you should be looking at a binary mask which you can then polygonize. I've added a full example – mikewatt Aug 24 '20 at 17:41
  • Your suggestion worked for me. Thanks for the help. – RRSC NGP Aug 26 '20 at 06:32