I was trying to display an old PRISM ASCII raster file, but didn’t want to have to load ArcGIS and then import the raster from the ASCII file. The new PRISM arrays are in a BIL format, which is much more convenient for use in ArcGIS, but requires using something like GDAL to read into NumPy arrays.
Reading the ASCII raster format directly into a NumPy array is really simple though. The format of the ASCII raster is defined as:
NCOLS xxx NROWS xxx XLLCENTER xxx | XLLCORNER xxx YLLCENTER xxx | YLLCORNER xxx CELLSIZE xxx NODATA_VALUE xxx row 1 row 2 ... row n
The only other issue with the old PRISM ASCII rasters is that they are stored as integers but scaled by 100.
import numpy as np
# Read in PRISM header data
with open('us_tmax_2011.14.asc', 'r') as prism_f:
prism_header = prism_f.readlines()[:6]
# Read the PRISM ASCII raster header
prism_header = [item.strip().split()[-1] for item in prism_header]
prism_cols = int(prism_header[0])
prism_rows = int(prism_header[1])
prism_xll = float(prism_header[2])
prism_yll = float(prism_header[3])
prism_cs = float(prism_header[4])
prism_nodata = float(prism_header[5])
# Read in the PRISM array
prism_array = np.loadtxt(prism_path, dtype=np.float, skiprows=6)
# Set the nodata values to nan
prism_array[prism_array == prism_nodata] = np.nan
# PRISM data is stored as an integer but scaled by 100
prism_array *= 0.01
The next step was using the imshow function in matplotlib to display the array. I thought a colorbar and grid lines would look nice also.
import matplotlib.pyplot as plt
# Plot PRISM array again
fig, ax = plt.subplots()
ax.set_title('PRISM 2011 Precipitation')
# Get the img object in order to pass it to the colorbar function
img_plot = ax.imshow(prism_array, cmap='jet')
# Place a colorbar next to the map
cbar = fig.colorbar(img_plot)
ax.grid(True)
plt.show()
The figure clearly shows the data, but I thought it would be nice if the x and y ticks where in degrees instead of the row and column indices. It would also be nice to have the colorbar match up with the y-axis of the plot and have a label.
I tried setting the ticks manually for awhile before I realized you could just pass an extent object to imshow. It is pretty easy to manually set the colorbar parameters as well, but I’m sure there is an automated way to do this.
# Calculate the extent of the PRISM array
# Extent is Xmin, Xmax, Ymin, Ymax
prism_extent = [
prism_xll, prism_xll + prism_cols * prism_cs,
prism_yll, prism_yll + prism_rows * prism_cs]
img_plot = ax.imshow(prism_array, extent=prism_extent)
# Place a colorbar next to the map
cbar = plt.colorbar(img_plot, orientation='vertical', shrink=0.5, aspect=14)
cbar.set_label('Degrees C')
If you want the colorbar to be on the bottom, that is a pretty easy fix as well.
cbar = plt.colorbar(img_plot, orientation='horizontal')
cbar.set_label('Degrees C')
PyMorton 