Image processing using Python, GDAL and Scikit-Image

Question

I am struggling with a processing and hopefully I will be able to solve here.

I work with Remote Sensing applied to Forestry, especially working with LiDAR data. The idea is to use Scikit-image for tree top detection. Since I'm new in Python, I considered a great personal triumph to do the following:

1. Import a CHM (with matplotlib);
2. Run a gaussian filter (with scikit-image package);
3. Run a maxima filter (with scikit-image package);
4. Run the peak_local_max (with scikit-image package);
5. Show the CHM with the local maxima (with matplotlib);

Now my problem. When I import with matplot, the image loses its geographic coordinates. So the coordinates I have are just basic image coordinates (i.e. 250,312). What I need is to get the value of the pixel under the local maxima dot in the image (red dots in the image). Here in the forum I saw one guy asking the same thing (Getting pixel value of GDAL raster under OGR point without NumPy?), but he already had the points in a shapefile. In my case the points were computed with scikit-image (It is an array with the coordinates of each tree top). So I do not have the shapefile.

In conclusion, what I want in the end is a txt file with the coordinates of each local maxima in geographic coordinates, for example:

525412 62980123 1150 ...

Answer 1

Firstly, welcome to the site!

Numpy arrays don't have a concept of coordinate systems inbuilt into the array. For a 2D raster they are indexed by column and row.

Note I'm making the assumption that you're reading a raster format that is supported by GDAL.

In Python the best way to import spatial raster data is with the rasterio package. The raw data imported by rasterio is still a numpy array without access to coordinate systems, but rasterio also gives you access to an affine method on the source array which you can use to transform raster columns and rows to projected coordinates. For example:

import rasterio

# The best way to open a raster with rasterio is through the context manager
# so that it closes automatically

with rasterio.open(path_to_raster) as source:

    data = source.read(1) # Read raster band 1 as a numpy array
    affine = source.affine

# ... do some work with scikit-image and get an array of local maxima locations
# e.g.
# maxima = numpy.array([[0, 0], [1, 1], [2, 2]])
# Also note that convention in a numy array for a 2d array is rows (y), columns (x)

for point in maxima: #Loop over each pair of coordinates
    column = point[1]
    row = point[0]
    x, y = affine * (column, row)
    print x, y

# Or you can do it all at once:

columns = maxima[:, 1]
rows = maxima[:, 0]

xs, ys = affine * (columns, rows)

And from there you can write your results our to a text file however you like (I'd suggest taking a look at the inbuilt csv module for example).

Answer 2

From a quick look at matplotlib, i'd say you have to alter the axis scales after the import.

Answer 3

Please try with following piece of code. This can be used to read image data from raster and write processed data to raster(.geotiff file).

from PIL import Image,ImageOps
import numpy as np
from osgeo import gdal
#from osgeo import gdal_array
#from osgeo import osr
#from osgeo.gdalconst import *
#import matplotlib.pylab as plt

#from PIL import Image, ImageOps
#import gdal
#from PIL import Image
gdal.AllRegister()

################## Read Raster #################
inRaster='C:\python\Results\Database\Risat1CRS\CRS_LEVEL2_GEOTIFF\scene_HH\imagery_HH.tif'

inDS=gdal.Open(inRaster,1)
geoTransform = inDS.GetGeoTransform()
band=inDS.GetRasterBand(1)
datatype=band.DataType
proj = inDS.GetProjection()
rows = inDS.RasterYSize
cols=inDS.RasterXSize
data=band.ReadAsArray(0,0,cols,rows)#extraction of data to be processed#
############write raster##########
driver=inDS.GetDriver()
outRaster='C:\\python\\Results\\Database\\Temporary data base\\clipped_26July2017\\demo11.tif'
outDS = driver.Create(outRaster, cols,rows, 1,datatype)
geoTransform = inDS.GetGeoTransform()
outDS.SetGeoTransform(geoTransform)
proj = inDS.GetProjection()
outDS.SetProjection(proj)
outBand = outDS.GetRasterBand(1)
outBand.WriteArray(data1,0,0)
#data is the output array to written in tiff file
outDS=None 
im2=Image.open('C:\\python\\Results\\Database\\Temporary data base\\clipped_26July2017\\demo11.tif');
im2.show()