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.

Improve this question
6
  • 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. Commented Feb 7, 2019 at 22:43
  • I've edited the question to be more specific Commented Feb 8, 2019 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. Commented Feb 10, 2019 at 23:38
  • You are right, I was working with different coordinate reference systems. Is there any way to correct without starting again? Commented Feb 11, 2019 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 Commented Feb 11, 2019 at 23:33

2 Answers 2

Reset to default
2

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

Improve this answer
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.

Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.