I am trying to take a raster (landsat_tm.img) and clip it to the spatial extent of the vector (clip.shp). My problem is that the band 1 values are being copied for all 7 bands.

    #import modules
    import numpy, os, ogr, osr, sys
    import gdal,gdalconst
    from gdalconst import *

    #set working directory

    #set drivers
    vector_driver=ogr.GetDriverByName("ESRI Shapefile")

    #define data sources

    #open data source
    raster_ds= gdal.Open(fn,GA_ReadOnly)

    #check if raster exists
    if raster_ds is None:
        print "Couldn't open."

    #get basic image info
    cols= raster_ds.RasterXSize
    rows= raster_ds.RasterYSize

    #get basic info on image coordinate and pixel resolution
    geotransform= raster_ds.GetGeoTransform()
    originX= geotransform[0]
    originY= geotransform[3]
    pixelWidth= geotransform[1]
    pixelHeight= geotransform[5]

    #open vector data source and read in extent
    vector_ds= ogr.Open(vector_fn)
    vector_layer = vector_ds.GetLayer()
    x_min, x_max, y_min, y_max = vector_layer.GetExtent()

    # Get vector x and y resolution info
    x_res = abs(int((x_max - x_min) / pixelWidth))
    y_res = abs(int((y_max - y_min) / pixelHeight))

    #get x and y offset
    xOffset=int((x_min - originX)/pixelWidth)
    yOffset=int((y_max- originY)/pixelHeight)

    #access bands in data source
    band= raster_ds.GetRasterBand(1)

    #transform image to array

    #create subset image
    subset=raster_driver.Create("subset.img", x_res, y_res, 7, GDT_Float32)

    #add data to blank subset image
    index = 0
    for i in range (bands):
        index = index + 1
        outband= subset.GetRasterBand (i + 1)
        outband.WriteArray(image_array, 0, 0)

    #set geotransform
    transformation=[x_min, geotransform[1], geotransform[2], y_min,geotransform[4], geotransform[5]]

    #set projection
    gdal.SetConfigOption("HFA_USE_RRD", "YES")

    #build pyramids
    subset= None

Update: For loop should look like this for anyone wondering.

    for i in range (raster_ds.RasterCount):
        i += 1
        inband = raster_ds.GetRasterBand(i)
        imgdata= inband.ReadAsArray(xOffset,yOffset,x_res,y_res).astype(numpy.float)
        subset_array= subset_ds.GetRasterBand(i)
        subset_array.WriteArray(imgdata, 0, 0)

This line gives you a hint: ValueError: expected array of dim 2.

This usually means that the data that you're writing doesn't have the right dimensions. You'd better inspect the variable image_array to check that it is a two dimensional array.

I'm not really familiar with how the previous functions work, you I can't help troubleshoot this. But do something like print image_array.shape to view the current shape of that array, and check that you're forming it right.

| improve this answer | |

The error is resulting from iterating through the bands and tryin to write all bands each time:

for i in range (bands):
    outband.WriteArray(image_array, 0, 0)

image_array contains all bands and is a numpy array with 3 dimensions. You need to slice it to the band you want to write to get a 2-dimensional array you can write into the band. i.e. image_array[0, :, :] is a 2D array containing all pixels from the first band.

Most probably you want something like this.

for i in range (bands):
    outband.WriteArray(image_array[i, :, :], 0, 0)
| improve this answer | |
  • yes, you're right. All bands are writing information from the first band. So, even though I got the subset image to output, the values for all bands are the same. When I tried your line of code, I got outband.WriteArray(image_array[i,:,:], 0, 0) IndexError: too many indices for array – CodingGrandma Feb 26 '18 at 19:25

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.