Image processing using Python, GDAL and Scikit-Image

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).