3

I have a BigTiff image and an associate point shapefile with features. I would like to know how can I create a bounding box from that shapefile to create/clip images around that points to train a model.

EDIT: I'm trying to do it automatically in a Jupyter Notebook and I've followed this threat but all exported images are completely black.

from osgeo import ogr, gdal

InputImage = 'XXX.tif'
Shapefile = 'XXX.shp'
RasterFormat = 'GTiff'
PixelRes = 0.5
VectorFormat = 'ESRI Shapefile'

# Open datasets
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()

VectorDriver = ogr.GetDriverByName(VectorFormat)
VectorDataset = VectorDriver.Open(Shapefile, 0) # 0=Read-only, 1=Read-Write
layer = VectorDataset.GetLayer()
FeatureCount = layer.GetFeatureCount()
print("Feature Count:",FeatureCount)

# Iterate through the shapefile features
Count = 0
for feature in layer:
    Count += 1
    print("Processing feature "+str(Count)+" of "+str(FeatureCount)+"...")

    geom = feature.GetGeometryRef() 
    minX, maxX, minY, maxY = geom.GetEnvelope()
    xmin, ymin = minX - 50, minY - 50
    xmax, ymax = maxX + 50, maxY + 50

    # Create raster
    OutTileName = str(Count)+'.SomeTileName.tif'
    OutTile = gdal.Warp(OutTileName, Raster, format=RasterFormat, outputBounds=[xmin, xmax, ymin, ymax], xRes=PixelRes, yRes=PixelRes, dstSRS=Projection, resampleAlg=gdal.GRA_NearestNeighbour, options=['COMPRESS=DEFLATE'])
    OutTile = None # Close dataset

# Close datasets
Raster = None
VectorDataset.Destroy()
print("Done.")

Moreover, I feel this is not the most simple or efficient way to do the bounding box.

  • Do you have any code to start from? In python you will want QgsVectorLayer.Extent qgis.org/pyqgis/master/core/… but first you might want to call updateExtents() qgis.org/pyqgis/master/core/… to ensure that the extents match the points in the file. – Michael Stimson Feb 7 at 22:43
  • I've edited the question to be more specific – Alex Feb 8 at 10:56
  • Please confirm that your points and your raster have the same coordinate reference system.. Although unlikely there might also be a problem with your compression, try with default (none) compression and see if that helps. – Michael Stimson Feb 10 at 23:38
  • You are right, I was working with different coordinate reference systems. Is there any way to correct without starting again? – Alex Feb 11 at 19:27
  • 1
    After geom = feature.GetGeometryRef() add a line geom.TransformTo(Projection) to transform the geometry to the same spatial reference obtained from the raster. For this to work the shapefile and raster must have a coordinate reference system set, the geometry will obtain the CRS from the layer it is derived from so it's not necessary to set it in the code. See example pcjericks.github.io/py-gdalogr-cookbook/projection.html – Michael Stimson Feb 11 at 23:33
2

From the Processing Toolbox run the "Extract layer extent" algorithm.

0

I advise you to use rasterio library for raster's manipulations. Here you have the command to easily extract the bounding box of your raster.

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.