Tag Archives: ASCII

Plotting PRISM ASCII arrays using Basemap

prism_projected
Plotting NumPy arrays directly using Matplotlib was simple, but what I really wanted though was to show the data in a projected coordinate system.  With a little bit of searching it seemed like the Python basemap module should be able to handle this task pretty easily.  The examples on the basemap website  were really good but the example on Peak5398’s Introduction to Python site also really helped with some of the details.

Just like last time, the first step was reading in the PRISM ASCII raster to a NumPy array.  The only difference is that I read the data into a masked array.

import numpy as np

# Read in PRISM header data
prism_name = 'us_tmax_2011.14.asc'
with open(prism_name, 'r') as prism_f:
    prism_header = prism_f.readlines()[:6]

# Pull out nodata value from 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])

# Readin the PRISM array
prism_array = np.loadtxt(prism_name, dtype=np.float, skiprows=6)

# Build a masked array
prism_ma = np.ma.masked_equal(prism_array, prism_nodata)

# PRISM data is stored as an integer but scaled by 100
prism_ma *= 0.01

# Calculate the extent of the PRISM array
# XMin, XMax, YMin, YMax
prism_extent = [
    prism_xll, prism_xll + prism_cols * prism_cs,
    prism_yll, prism_yll + prism_rows * prism_cs]

I choose to use an Albers Equal Area projection, but it seems really simple to set just about any projection.  I used a masked array so that I could use the maskoceans function to remove the PRISM values over the ocean.  Just setting these cells to NaN didn’t work.

Overall it seems really easy to plot array data in a projected coordinate system.

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, maskoceans

# Build arrays of the lat/lon points for each cell
lons = np.arange(prism_extent[0], prism_extent[1], prism_cs)
lats = np.arange(prism_extent[3], prism_extent[2], -prism_cs)
lons, lats = np.meshgrid(lons, lats)

# Just the contingous United States
m = Basemap(
    width=5000000, height=3400000, projection='aea', resolution='l',
    lat_1=29.5, lat_2=45.5, lat_0=38.5, lon_0=-96.)

# Build the figure and ax objects
fig, ax = plt.subplots()
ax.set_title('PRISM 2011 Tmax')

# Built-in map boundaries
m.drawcoastlines()
m.drawlsmask(land_color='white', ocean_color='0.8', lakes=False)
m.drawcountries()
m.drawstates()

# Set ocean pixels to nodata
prism_ma_land = maskoceans(lons, lats, prism_ma)

# Draw PRISM array directly to map
im = m.pcolormesh(lons, lats, prism_ma_land, cmap=plt.cm.jet, latlon=True)

# Build the colorbar on the right
cbar = m.colorbar(im, 'right', size='5%', pad='5%')
# Build the colorbar on the bottom
# cbar = m.colorbar(im, 'bottom', size='5%', pad='8%')
cbar.set_label('Degrees C')

# Meridians & Parallels
m.drawmeridians(np.arange(0, 360, 10), labels=[False,False,False,True])
m.drawparallels(np.arange(-90, 90, 10), labels=[True,False,False,False])

# Save and show
plt.show()

 

Plotting PRISM ASCII arrays using Matplotlib imshow

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

 

 

prism_geographic_a

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

prism_geographic_b

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

prism_geographic_c