0

I want to work with TIFF files of aeronautical charts supplied by the FAA here.

The Python module rasterio seems to be the right tool for my purposes and I’ve been through several tutorials and reference docs. They all work for me, but not with the the data I want to use. I cannot seem to get it to display properly in color.

The main difference seems to be that my data is single band (as reported by dataset.count) whereas all the tutorials have three or more bands.

For example, the following code should display the chart:

import rasterio
from rasterio.plot import show
dataset = rasterio.open ('Green_Bay_101/Green Bay SEC 101.tif')
show(dataset)

The resulting plot has the colors all wrong.

Screen shot of image produced by above code

But it should look like this:

Screen shot of image shown with macOS Finder

I’m sure that I’m missing something basic, but I cannot find it. Would someone point it out?

I first posted this on StackOverflow but I think that's the wrong place.

Vince
  • 20,017
  • 15
  • 45
  • 64
Harry Dolan
  • 21
  • 2
  • Welcome to GIS SE. Please avoid cross-posting (https://stackoverflow.com/q/64578207/1446289) for the following reasons: https://meta.stackexchange.com/q/64068/641151. Instead, choose the site that is the best fit and delete the other. – Aaron Oct 29 '20 at 02:03

3 Answers3

2

Yes, there are better ways. The following worked for me:

import rasterio
import matplotlib.pyplot as pyplot
import numpy as np

dataset = rasterio.open('Some_SEC_Chart_from_FAA.tif')
image2 = dataset.read(1)
color_array = np.asarray(
    [dataset.colormap(1)[i] for i in range(256)], dtype=np.uint8)
image3 = color_array[image2]
rasterio.plot.show(np.rollaxis(image3, 2), transform=dataset.transform)

Same idea as yours, it creates separate array image3, but works faster and uses rasterio.plot.show with transform to show geocoded extents.

If you want to use image2 directly and keep low memory, you can define the colormap as you mentioned:

from matplotlib.colors import ListedColormap

rasterio.plot.show(image2, transform=dataset.transform,
                   cmap=ListedColormap(color_array/255),
                   interpolation='none')

I turn off interpolation because the custom colormap from the dataset is not color-continuous, otherwise some weird colors show up at the edges.

Michael
  • 21
  • 3
1

I have run more experiments and determined that this has to do with the difference between how matplotlib and rasterio handle the image. The following code demonsrates the problem. matplotlib does it properly but rasterio doesn't.

import rasterio
import matplotlib.pyplot as pyplot

Obtain from https://aeronav.faa.gov/content/aeronav/sectional_files/Anchorage_106.zip


file = 'Anchorage-Fairbanks_TAC_85/Fairbanks TAC 85.tif'


image1 = pyplot.imread(file)
pyplot.imshow(image1)
pyplot.show()


dataset = rasterio.open(file)
image2 = dataset.read(1)
pyplot.imshow(image2)
pyplot.show()

Here's what matplotlib produces:

matplotlib plot

And here it is from rasterio:

rasterio plot

I suspect it's a colormap issue but I don't know what to do.

Harry Dolan
  • 21
  • 2
  • Now I am certain that it’s a color map issue. I printed out and examined the (huge) image1 and image2 arrays. The first seems to be an array of pixel values (each a group of 4 8-bit values) and the second a array of single 8-bit integers. When I lined them up I discovered that the image2 values are in fact indices into the corresponding pixel values. Then I verified this by looking at dataset.colormap(1) and found the exact same table.

    It seems the I should be able to give the rasterio.imshow() call a parameter to tell it to use this colormap.

     – Harry Dolan Oct 30 '20 at 21:09
1

OK, the following code displays the image properly, but it’s incredibly slow and kludgey. Since I don’t know how to tell imshow to use the colormap from the dataset, I instead create a new image file and map the colors myself.

There’s got to be a better way!

import rasterio
import matplotlib.pyplot as pyplot
import numpy as np

Obtain from https://aeronav.faa.gov/content/aeronav/sectional_files/Anchorage_106.zip


file = 'Anchorage-Fairbanks_TAC_85/Fairbanks TAC 85.tif'


dataset = rasterio.open(file)
image2 = dataset.read(1)


#############################################################################


Start kludge.  Create image3 from image2.


colors = dataset.colormap(1).values()
mycolors = np.zeros((len(colors),4),dtype=np.uint8)
i=0
for color in colors:
    mycolors[i] = np.array(list(color))
    i += 1


shape3 = (image2.shape[0], image2.shape[1], 4)
image3 = np.zeros(shape3, dtype=np.uint8)
for i in range(image2.shape[0]):
    for j in range(image2.shape[1]):
        val = image2[i][j]
        image3[i][j] = mycolors[val]


End kludge


#############################################################################


pyplot.imshow(image3)  # Showing image3 instead of image2


pyplot.show()
Harry Dolan
  • 21
  • 2