My objective is to classify a Landsat 8 image using supervised classification. I have a training point shapefile data and a preprocessed Imagery. I want to use Random Forest classification method.

I want to rasterize all the vectors (actually, there is only one point shapefile) in the given directory into a single image using GDAL. The labeled_pixels are those whose land use land cover extracted (rasterized) from training vector (shapefile) data. The labeled_pixels have 5 classes [1,2,3,4,5]. However when trying to execute the following code from this source, I get

Traceback (most recent call last):
  File "C:/Users/Kaleab/Desktop/Machine_Learning/classify.py", line 73, in <module>
    labeled_pixels = vectors_to_raster(shapefiles, rows, cols, geo_transform, proj)
  File "C:/Users/Kaleab/Desktop/Machine_Learning/classify.py", line 35, in vectors_to_raster
    labeled_pixels += band.ReadAsArray()
ValueError: operands could not be broadcast together with shapes (1,2) (1678,2218)


import numpy as np
import os

from osgeo import gdal,ogr
from sklearn import metrics
from sklearn.ensemble import RandomForestClassifier

# A list of "random" colors (for a nicer output)
COLORS = ["#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941"]

def create_mask_from_vector(vector_data_path, cols, rows, geo_transform,
                            projection, target_value=1):
    """Rasterize the given vector (wrapper for gdal.RasterizeLayer)."""
    data_source = ogr.Open('./Springfield_UrbanLC/training/training_data.shp')
    if data_source is None:
        print ("Could not open data_source")
    layer = data_source.GetLayerByIndex(0)
    driver = gdal.GetDriverByName('MEM')  # In memory dataset
    target_ds = driver.Create('', cols, rows, 1, gdal.GDT_UInt16)
    gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[target_value])
    return target_ds

    def vectors_to_raster(file_paths, rows, cols, geo_transform, projection):
        """Rasterize all the vectors in the given directory into a single image."""
        labeled_pixels = np.nonzero((rows, cols))
        print ( " Rows", rows,"Columns", cols)
        for i, path in enumerate(file_paths):
            label = i+1
            ds = create_mask_from_vector(path, cols, rows, geo_transform,
                                         projection, target_value=label)
            band = ds.GetRasterBand(1)
            labeled_pixels += band.ReadAsArray()
            ds = None
        return labeled_pixels

def write_geotiff(fname, data, geo_transform, projection):
    """Create a GeoTIFF file with the given data."""
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    dataset = driver.Create(fname, cols, rows, 1, gdal.GDT_Byte)
    band = dataset.GetRasterBand(1)
    dataset = None  # Close the file

raster_data_path = './Springfield_UrbanLC/Landsat8_GeoTiff_Image/Sangamon_Image_Clip.tif'
output_fname = 'RFclassification.tiff'
train_data_path = './Springfield_UrbanLC/training'
validation_data_path = './Springfield_UrbanLC/test'

raster_dataset = gdal.Open(raster_data_path, gdal.GA_ReadOnly)
geo_transform = raster_dataset.GetGeoTransform()
proj = raster_dataset.GetProjectionRef()
bands_data = []
for b in range(1, raster_dataset.RasterCount+1):
    band = raster_dataset.GetRasterBand(b)

bands_data = np.dstack(bands_data)
rows, cols, n_bands = bands_data.shape

files = [f for f in os.listdir(train_data_path) if f.endswith('.shp')]
classes = [f.split('.')[0] for f in files]
shapefiles = [os.path.join(train_data_path, f)
              for f in files if f.endswith('.shp')]

labeled_pixels = vectors_to_raster(shapefiles, rows, cols, geo_transform, proj)
is_train = np.nonzero(labeled_pixels)
training_labels = labeled_pixels[is_train]
training_samples = bands_data[is_train]

classifier = RandomForestClassifier(n_jobs=1, n_estimators=10)
classifier.fit(training_samples, training_labels)

n_samples = rows*cols
flat_pixels = bands_data.reshape((n_samples, n_bands))
result = classifier.predict(flat_pixels)

classification = result.reshape((rows, cols))
write_geotiff(output_fname, classification, geo_transform, proj)

shapefiles = [os.path.join(validation_data_path, "%s.shp"%c) for c in classes]
verification_pixels = vectors_to_raster(shapefiles, rows, cols, geo_transform, proj)
for_verification = np.nonzero(verification_pixels)
verification_labels = verification_pixels[for_verification]
predicted_labels = classification[for_verification]

print("Confussion matrix:\n%s" %
      metrics.confusion_matrix(verification_labels, predicted_labels))
target_names = ['Class %s' % s for s in classes]
print("Classification report:\n%s" %
      metrics.classification_report(verification_labels, predicted_labels,
print("Classification accuracy: %f" %
      metrics.accuracy_score(verification_labels, predicted_labels))

I have limited knowledge about NumPy, however my online search indicates a broadcasting rule issue in NumPy, where the shapes of the training shapefile are not the same shape (not compatible) as those of each band in the mask raster. Can you please explain the issue and suggest a practical solution?

  • Can you edit your question and post the full exception, not just the ValueError message. – user2856 Jul 20 '17 at 11:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.