[GIS] Graph points in image using matplotlib

matplotlibpythonrasterioshapely

I am trying to plot points (latitude, longitude) inside a raster image in Python. First, I read the image using rasterio and plot it using the same package, after transform the image from uint16 to Byte (uint6): 

#Transform to uint8
def scale8bit(image):
    scale = float(256) / (image.max() - image.min())
    return np.clip(np.round(np.multiply(image, scale)), 0, 255).astype(np.uint8)
#Read image and plot visible bands    
src = rasterio.open(path)
array_8b = img_as_ubyte(img_as_float(src.read([3, 2, 1])))
a = scale8bit(src.read([3, 2, 1]))

Now, I have some WKB points which coincide –as far as I can tell from the data– with the bounding box of the image. Hence, I am trying to run the following in order to visualize the images and the points (read using shapely):

from shapely import wkb
from rasterio import plot.show
for i, geom in zip(selected_images_paths[0:10], subset_images["geom"][0:10]):

  #Open file using rasterio and change to int8
  src = rasterio.open(i)
  src_8bit = scale8bit(src.read([3, 2, 1]))

  #Graph!
  fig, ax = plt.subplots(figsize=(30, 30))
  plot.show(src_8bit)

  #Open geom with shapely
  geom_wkb = wkb.loads(geom, hex = True)
  x,y = geom_wkb.xy
  ax.scatter(x, y, marker = 'o', s=50, c='r')

The units of the graphs are clearly on color intensity and not in long, lat, which are the units of the x,y tuple. I've tried different options, but I can figure out how to overlap both data sets.

Is there something wanting, or should I take another approach?

EDIT
I have been following the code described in this question (Plot shapefile on top of raster using plot and imshow from matplotlib). Nonetheless, the result is sub-optimal.

When I use ax = plt.gca() to set the axes in the subplot, the fig definition is completely ignored for some images and the resulting figure is small. I think the error must be on the ax or fig definition, but I've tried other suggestions in the documentation without success.

Best Answer

There might be multiple problems here.

On the larger problem of sharing axes or making rasterio.plot.show not also invoke pyplot.show (which renders the current figure to screen). If you look at the rasterio docs on plotting, you will see that rasterio.plot.show takes a keyword argument called ax which allows for specifying an existing axis. You should pass in the ax variable you create with plt.subplots.

I think you are also encountering a problem with coordinates. rasterio.plot.show is a wrapper around matplotlib.pyplot.imshow. This function plots an image in the image space, not the geographic space. That means the image does not have georeferencing as might be expected. The transform argument to plot.show just handles formatting the labels for the axes; it does not shift the location of the image in the plotting space. Therefore attempting to plot the corners of the image as longitude-latitude pairs will actually result in plotting the points as pixel coordinates. You will need to either transform the lon-lats into pixel coordinates using the inverse geotransform from the image or place the image into its proper geolocation with pyplot.imshow(img, extent=[left, right, bottom, top]) (see an example question that is similar to yours here).

Related Question