I have Landsat 8 preprocessed image I want to classify using random forest(RF) classification in python. The size of the image is 3,721,804 pixels with 7 bands. I have 250 training data shapefiles which were rasterized and yielded y (labels) and trainingData. The training labels(y) have five classes [1,2,3,4,5] with (250,) dimension. The trainingData has a dimension of (250,7) i.e. 250 sampling rasters each with 7 band.
I am encountering memoryerror when trying to predict the classes of the unknown (class)pixels because the process consumes so much RAM memory attributable to the byte-size of noSamples (rows x columns) in the input Landsat Image. I have read that Random Forest does not do well on sparse matrix (+1000 columns). But others have produced a classification using RF method.
How can I solve this problem?
The Code:
import numpy as np import os import pandas as pd from osgeo import gdal,gdal_array,ogr from sklearn import metrics from sklearn.ensemble import RandomForestClassifier from matplotlib import pyplot as plt # Get the current working directory print(os.getcwd()) # Open the dataset from the file dataset = ogr.Open('./training/training_data.shp') # Make sure the dataset exists -- it would be None if we couldn't open it if not dataset: print('Error: could not open dataset') #Let's get the driver from this file driver = dataset.GetDriver() print('Dataset driver is: {n}\n'.format(n=driver.name)) #How many layers are contained in this Shapefile? layer_count = dataset.GetLayerCount() print('The shapefile has {n} layer(s)\n'.format(n=layer_count)) #What is the name of the 1 layer? layer = dataset.GetLayerByIndex(0) print('The layer is named: {n}\n'.format(n=layer.GetName())) # Get the spatial reference spatial_ref = layer.GetSpatialRef() # Export this spatial reference to something we can read... like the Proj4 proj4 = spatial_ref.ExportToProj4() print('Layer projection is: {proj4}\n'.format(proj4=proj4)) ### How many features are in the layer? feature_count = layer.GetFeatureCount() print('Layer has {n} features\n'.format(n=feature_count)) ### How many fields are in the shapefile, and what are their names? # First we need to capture the layer definition defn = layer.GetLayerDefn() # How many fields field_count = defn.GetFieldCount() print('Layer has {n} fields'.format(n=field_count)) # What are their names? print('Their names are: ') for i in range(field_count): field_defn = defn.GetFieldDefn(i) print('\t{name} - {datatype}'.format(name=field_defn.GetName(), datatype=field_defn.GetTypeName())) #In order to pair up our vector data with our input image raster pixels, #we will need a way of co-aligning the datasets in space. #First we will open our input raster image, to understand how we will want to rasterize our training vector raster_ds = gdal.Open('./Landsat8_GeoTiff_Image/Sangamon_Image_Clip.tif', gdal.GA_ReadOnly) # Fetch number of rows and columns ncol = raster_ds.RasterXSize nrow = raster_ds.RasterYSize # Fetch projection and extent proj = raster_ds.GetProjectionRef() ext = raster_ds.GetGeoTransform() raster_ds = None # Create the raster dataset memory_driver = gdal.GetDriverByName('GTiff') out_raster_ds = memory_driver.Create('./training/training_data.gtif', ncol, nrow, 1, gdal.GDT_Byte) # Set the ROI image(rasterized training vector)'s projection and extent to our input raster's projection and extent out_raster_ds.SetProjection(proj) out_raster_ds.SetGeoTransform(ext) # Fill our output band with the 0 blank, no class label, value b = out_raster_ds.GetRasterBand(1) b.Fill(0) # Rasterize the shapefile layer to our new dataset status = gdal.RasterizeLayer(out_raster_ds, # output to our new dataset [1], # output to our new dataset's first band layer, # rasterize this layer None, None, # don't worry about transformations since we're in same projection [0], # burn value 0 ['ALL_TOUCHED=TRUE', # rasterize all pixels touched by polygons 'ATTRIBUTE=itree_New_'] # put raster values according to the 'itree_New_' field values ) # Close dataset out_raster_ds = None # Check rasterized layer roi_ds = gdal.Open('./training/training_data.gtif', gdal.GA_ReadOnly) roi = roi_ds.GetRasterBand(1).ReadAsArray() if status != 0: print("I don't think it worked...") else: print("Success") # How many pixels are in each class? classes = np.unique(roi) # Iterate over all class labels in the ROI image, printing out some information for c in classes: print('Class {c} contains {n} pixels'.format(c=c, n=(roi == c).sum())) # Tell GDAL to throw Python exceptions, and register all drivers gdal.UseExceptions() gdal.AllRegister() #Read in our image and ROI image img_ds = gdal.Open('./Landsat8_GeoTiff_Image/Sangamon_Image_Clip.tif', gdal.GA_ReadOnly) roi_ds = gdal.Open('./training/training_data.gtif', gdal.GA_ReadOnly) img = np.zeros((img_ds.RasterYSize, img_ds.RasterXSize, img_ds.RasterCount), gdal_array.GDALTypeCodeToNumericTypeCode(img_ds.GetRasterBand(1).DataType)) for b in range(img.shape[2]): img[:, :, b] = img_ds.GetRasterBand(b + 1).ReadAsArray() roi = roi_ds.GetRasterBand(1).ReadAsArray().astype(np.uint8) # In order to plot a multiband image , we need to normalize # the array of each band and also, you should convert the data to float32 or uInt8 for matplotlib. def normalize(img): ''' Function to normalize an input array to 0-1 ''' img_min = img.min() img_max = img.max() return (img - img_min) / (img_max - img_min) normalize(img) index=np.array([4,3,2]) colors=img[:,:,index].astype(np.float64) # Display them plt.subplot(121) plt.imshow(colors,cmap=plt.cm.Spectral) plt.colorbar() #plt.imshow(img[:, :, 4], cmap=plt.cm.Greys_r) plt.title('Normalized Multiband Image') plt.subplot(122) plt.imshow(roi,cmap=plt.cm.Greys_r) plt.title('ROI Training Data') plt.show() def createGeotiff(outRaster, data, geo_transform, projection): # Create a GeoTIFF file with the given data driver = gdal.GetDriverByName('GTiff') rows, cols = data.shape rasterDS = driver.Create(outRaster, cols, rows, 1, gdal.GDT_Byte) rasterDS.SetGeoTransform(geo_transform) rasterDS.SetProjection(projection) band = img_ds.GetRasterBand(1) band.WriteArray(data) dataset = None '''Now, We have the ROI Image(rasterized training point shapefiles) and the input multiband image, let's pair Y with X. The image we want to classify is our X feature inputs, and the ROI with the land cover labels is our Y labeled data, we need to pair them up in NumPy arrays so we may feed them to Random Forest:''' # Find how many non-zero entries we have -- i.e. how many training data samples? n_samples = (roi > 0).sum() print('We have {n} samples'.format(n=n_samples)) # What are our classification labels? labels = np.unique(roi[roi > 0]) print('The training data include {n} classes: {classes}'.format(n=labels.size, classes=labels)) # We will need a "X" matrix containing our features, and a "y" array containing our labels # These will have n_samples rows # In other languages we would need to allocate these and them loop to fill them, but NumPy can be faster X = img[roi > 0, :] # include 8th band, which is Fmask, for now y = roi[roi > 0] # The matrix is number of samples(250) x number of bands(7) print('Our X matrix is sized: {sz}'.format(sz=X.shape)) print('Our y array is sized: {sz}'.format(sz=y.shape)) '''Training the Random Forest Now that we have our X matrix of feature inputs (the spectral bands) and our y array (the labels), we can train our model.''' # Initialize our model with 500 trees rf = RandomForestClassifier(n_estimators=1000, oob_score=True) # Fit our model to training data rf = rf.fit(X, y) print('Our OOB("Out-Of-Bag") prediction of accuracy is: {oob}%'.format(oob=rf.oob_score_ * 100)) '''To help us get an idea of which spectral bands were important, we can look at the feature importance scores:''' bands = [1, 2, 3, 4, 5, 7, 6] for b, imp in zip(bands, rf.feature_importances_): print('Band {b} importance: {imp}'.format(b=b, imp=imp)) # Extract band's data and transform into a numpy array bandsData = [] for b in range(1, img_ds.RasterCount+1): band = img_ds.GetRasterBand(b) bandsData.append(band.ReadAsArray()) bandsData = np.dstack(bandsData) rows, cols, noBands = bandsData.shape # Prepare training data (set of pixels used for training) and labels # The nonzero pixels are training data print("Training Data dimension",X.shape) print("Labels data dimension",y.shape) # Train a Random Forest classifier classifier = RandomForestClassifier(n_jobs=1, n_estimators=10) classifier.fit(X, y) # Predict class label of unknown pixels noSamples = rows*cols print ("The number of pixel Samples from the input image is:", noSamples,"with", noBands, "bands") flat_pixels = bandsData.reshape((noSamples, noBands)) result = classifier.predict(flat_pixels) classification = result.reshape((rows, cols)) # Create a GeoTIFF file with the given data createGeotiff(outRaster, classification, geo_transform, projection) img = Image.open('randomForest.tiff') img.save('randomForest.png','png') 
