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.

Alex
  • 41
  • 4
  • Do you have any code to start from? In python you will want QgsVectorLayer.Extent https://qgis.org/pyqgis/master/core/QgsVectorLayer.html#qgis.core.QgsVectorLayer.extent but first you might want to call updateExtents() https://qgis.org/pyqgis/master/core/QgsVectorLayer.html#qgis.core.QgsVectorLayer.updateExtents to ensure that the extents match the points in the file. – Michael Stimson Feb 07 '19 at 22:43
  • I've edited the question to be more specific – Alex Feb 08 '19 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 '19 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 '19 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 https://pcjericks.github.io/py-gdalogr-cookbook/projection.html – Michael Stimson Feb 11 '19 at 23:33
  • I have the shapefile and the raster in the same coordinate system (ETRS89 / UTM zone 31N; EPSG:25831) but when I run the code I got this error log:

    ERROR 3: Free disk space available is 330339983360 bytes, whereas 195259767827628 are at least necessary. You can disable this check by defining the CHECK_DISK_FREE_SPACE configuration option to FALSE.

    Any idea?

    – Alex Feb 12 '19 at 11:14

2 Answers2

2

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

ndawson
  • 27,620
  • 3
  • 61
  • 85
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.

Tim C.
  • 1,506
  • 4
  • 22
  • 32