TransWikia.com

Create polygon after excluding extents from raster data having no-data values

Geographic Information Systems Asked on September 5, 2020

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:UsersDesktopPRODUCT1LandsatT1.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:UsersDesktopPRODUCT1LandsatT1_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.

One Answer

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

Answered by mikewatt on September 5, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP