[GIS] Suitable Raster processing methods to extract field boundaries from images

digitizinggimpgrassqgisraster

I'm trying to extract field boundaries and features from a lot of aerial images of agricultural land – and avoid digitizing these features manually. I've managed to get a crude raster file of the features that I want by manipulating the threshold colour values in GIMP and geocoding the resulting .tif in QGIS.

Farm outline and resulting Raster overlay

The nine broad contour lines are what I want to extract (and clean up a little) but more importantly I want to remove the randomly distributed noise seen in the middle of the fields. After cleaning the image I plan to convert to a vector and clip the existing farm boundary with it.

I am at a loss as to how to process these data?

I have both QGIS 2.8.1 and GRASS 7 available.

Best Answer

I'm not sure how well this will work without trying on your specific raster, but a few pointers. I assume from your image you used Threshold in GIMP to give you a binary raster (inside field/outside field)

Removing noise

You can use GDAL Sieve to remove the 'salt and pepper' noise - those little holes within the fields.

You can do this from QGIS (Raster > Analysis > Sieve). That should fill any 'holes' (smaller than an area you specify) with the value of the surrounding raster area. If you're lucky and guess the minimum area well, it might also clean up the noise in the gaps between fields too.

Using QGIS to get edges

Depending on how successful that is, and how well you choose the sieve settings, you may be able to convert the noise-free raster straight to polygons. There seems to be a lot of noise in the gaps between fields so there's a good chance you end up with some fields joined together.

Using SAGA to get edges

If you get a reasonably clean noise-free image from Sieve

SAGA GIS has a Raster Connectivty Filter (Grid > Filter > Connectivity Analysis). This will give two outputs

a raster, with each contiguous area of pixels given its own ID
a shapefile with the boundary of each contiguous area of pixels

Using GIMP to get edges

GIMP has a number of noise-removal and edge filters (e.g. Sobel, Difference of Guassians etc), you might want to try these too.

Post-processing

You might also need to simplify geometries to get rid of the 'staircasing' (the polygon edges will trace round pixel boundaries)

Related Solutions

[GIS] How to extract percentage of land cover from a raster according to multiple shapefiles

One solution would be to convert each of your polygons to rasters. Then take the portion of cell counts in each landcover class by the total cell count from the clip. Here is an example function:

    PortionClassInPoly <- function(MySingleShape, MainRaster) {

    # Start with one shape from your list
    # Rasterize this shape according to the parameters of your background raster
       mini_rast <- rasterize(MySingleShape, MainRaster, fun='last')  
    # divide by itself - to make new set new raster values to 1 (NA values stays as NA). 
       mini_rast <- mini_rast/mini_rast
    # get the total number of cells for this shape (save for below)
       total_cells <- cellStats(mini_rast, 'sum')

    # multiply by original - to make a mask layer from your original raster
       my_cutout <- MainRaster*mini_rast
    # Count the number of cells for each discrete landcover class
       in_this_poly_unit <- freq(my_cutout)
    # Divide the cell count for each of these classes by the total 
    # number of cells in your current shape
       class_percents <- in_this_poly_unit[,2]/total_cells

    # make the function return a dataframe with
    # % landcover in each class for each shape
       output <- data.frame(myclass=in_this_poly_unit[,1], portion=class_percents)

       return(output)
    }

With your data, copy and paste the above & then try this function in your loop:

    PortionClassInPoly(MySingleShape=list.shape[[k]], MainRaster=r.landcover)

    # you should get something like this:
    # percentage of every land cover for a shapefile 

      myclass   portion
    1       1 0.1726333
    2       2 0.1662689
    3       3 0.1718377
    4       4 0.1789976
    5       5 0.1662689
    6       6 0.1439936

[GIS] How to extract coordinates from a raster in QGIS

Yes, it is possible to extract a shapefile but it is very difficult to handle it because you need a lot a memory. It is preferable to "burn" these coordinates in two raster (one for each coordinates values; x and y).

However, I tried to produce the "shapefile" (in fact an "event theme" from a text file) with next approach. First, I selected an arbitrary area (at Brazil) with your exact cell number (2105x1388 cells) and helped with my two plugins (as it can be observed at the next image).

With the next code ran at the Python Console of QGIS, I got a text file to be loaded in QGIS.

layer = iface.activeLayer()

extent = layer.extent()

xmin = extent.xMinimum()
ymax = extent.yMaximum()

rows = layer.height()
columns = layer.width()

xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

k = 1

xinit = xmin + xsize/2
yinit = ymax - ysize/2

pfile = open('points.txt', 'w')

pfile.write('id x y \n')

print "Processing..."

for i in range(rows):
    for j in range(columns):
        x = xinit + j*xsize 
        y = yinit
        pfile.write(str(k) + " " + str(x) + " " + str(y) + "\n")
        k +=1
    xinit = xmin + xsize/2
    yinit -= ysize

pfile.close()

At the next image, you can observe a point layer (cell centres in lon/lat) created as delimited text file but, it is very difficult to handle. It is possible but inconvenient.

Editing Note:

Next code add the raster values to the text file:

from osgeo import gdal
import os

layer = iface.activeLayer()

provider = layer.dataProvider()

my_path = provider.dataSourceUri()

dataset = gdal.Open(my_path)
band = dataset.GetRasterBand(1)

data = band.ReadAsArray(0, 0, band.XSize, band.YSize)

extent = layer.extent()

xmin = extent.xMinimum()
ymax = extent.yMaximum()

rows = layer.height()
columns = layer.width()

xsize = layer.rasterUnitsPerPixelX()
ysize = layer.rasterUnitsPerPixelY()

k = 1

xinit = xmin + xsize/2
yinit = ymax - ysize/2

pfile = open('points.txt', 'w')

pfile.write('id x y \n')

print "Processing..."

for i in range(rows):
    for j in range(columns):
        x = xinit + j*xsize 
        y = yinit
        pfile.write(str(k) + " " + str(x) + " " + str(y) + " " + str(data[i][j]) + "\n")
        k +=1
    xinit = xmin + xsize/2
    yinit -= ysize

pfile.close()

dataset = None

print "Done!"

Best Answer

Related Solutions

[GIS] How to extract percentage of land cover from a raster according to multiple shapefiles

[GIS] How to extract coordinates from a raster in QGIS

Related Question