I have been working to display an image of a USGS 7.5 minute quad sheet.
These are provided at various locations about the Web. Since the
range of colors on these maps is limited, the *.tif files appear to
use an indexed color map wherein each pixel has a value 0 to 255 and
the color is found from a table with 256 entries having triplets of
the RGB in the range of 0-255. I have not been able to sort out
how to get the gdal package to give me the color map from within python, so I dumped it
from the image file using gdalinfo and then cut and pasted to get
the following script:
···
---------------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import ListedColormap
import osgeo.gdal as gdal
from osgeo.gdalconst import *
gd = gdal.Open('o37122d1.tif')
#Setup to compute the colormap from the RGB triplets from the geotif file
div = np.zeros( (256,3), np.float32)
div = 255.0
ctab = np.array([
[ 255,255,255],
[ 0,0,0],
[ 255,255,255],
[ 91,159,230],
[ 230,45,30],
[ 162,96,71],
[ 210,255,177],
[ 197,101,197],
[ 255,240,0],
[ 202,225,245],
...snip.... (deleted many lines:)
[ 250,202,250],
[ 230,230,230],
[ 222,167,146],
[ 255,255,255],
[ 255,255,255] ], dtype=np.uint8)
#Compute colors in range 0.0 to 1.0.
fctab= ctab/div
usgscm = ListedColormap(fctab, name='usgs',N=None)
doq = gd.ReadAsArray()
# doq proves to be a uint8 array: 8802 rows and 7066 columns
# Cut out a subset from the main array--to large to display and
# slow as 'molasses in January' to process:)
suba = doq[0:1101 ,0:1801 ]
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(suba, cmap=usgscm, origin='upper')
plt.show()
---------------------------------------------------------------------------------
This script does give me an image--but in badly wrong colors:( The script does
properly display gray-scaled Digital Ortho-quadrangles using cm.gray as the color
map in imshow. Consequently something is not quite correct with respect to the
definition or the use of the color map. It appears that each map, and there
are about 56,000 of them available on one site, could have its own color map.
Thus my application must be able to compute a color map unique to each of the
topographic maps.
Questions:
1. What am I missing to get imshow to pick out the correct colors from the
color map?
2. Should I be using the name, usgs, given in the ListedColormap instance someplace?
Not clear to me what role that name plays.
3. Is there a way to use ctab directly in the ListedColormap instance? Class Colormap
has a bytes argument, False by default, but I am not yet sure if it has any bearing
on my problem.
I am using Matplotlib 0.98 with the latest numpy, and gdal. I am not using
Basemap, after some study, because it contains far more "horsepower" than my
simple topographic-map viewer application requires. Also, this is my first
foray into "image processing". I am finding the learning curve a bit steep,
as usual, with multiple names for the same thing popping up in descriptions
from various sources--business as usual in the computer business:)
Thanks from a happy matplotlib user in California!
Delbert