0

I have a shapefile of the United states, and I have an m x n array of Cartesian data that represents temperature at each pixel. I am able to load in the shapefile and plot it:

import shapefile as shp
import matplotlib.pyplot as plt

sf = shp.Reader("/path/to/USA.shp")

plt.figure()
for shape in sf.shapeRecords():
    for i in range(len(shape.shape.parts)):
        i_start = shape.shape.parts[i]
        if i==len(shape.shape.parts)-1:
            i_end = len(shape.shape.points)
        else:
            i_end = shape.shape.parts[i+1]
        x = [i[0] for i in shape.shape.points[i_start:i_end]]
        y = [i[1] for i in shape.shape.points[i_start:i_end]]
        plt.plot(x,y, color = 'black')
plt.show()

And I am able to read in my data and plot it:

import pickle
from matplotlib import pyplot as mp
Tfile = '/path/to/file.pkl'
with open(Tfile) as f:
    reshapeT = pickle.load(f)
mp.matshow(reshapeT)

The problem is reshapeT has dimensions of 536 x 592, and is a subdomain of the US. However, I do have information about the top-left corner of the reshapeT grid (lat / long) as well as the spacing between each pixel (0.01)

My question is: How do I overlay the reshapeT data ontop of the shapefile domain?

WX_M
  • 458
  • 4
  • 20

1 Answers1

1

If I understand you correctly you would like to overlay a 536x592 numpy array over a specifc part of a plotted shapefile. I would suggest you use Matplotlib's imwshow() method, with the extent parameter, which allows you to place the image within the plot.

Your way of plotting the shapefile is fine, however, if you have the possibility to use geopandas, it will dramatically simplify things. Plotting the shapefile will reduce to the following lines:

import geopandas as gpd
sf = gpd.read_file("/path/to/USA.shp")
ax1 = sf.plot(edgecolor='black', facecolor='none')

As you have done previously, let's load the array data now:

import pickle
Tfile = '/path/to/file.pkl'
with open(Tfile) as f:
    reshapeT = pickle.load(f)

Now in order to be able to plot numpy array in the correct position, we first need to calculate its extent (the area which it will cover expressed in coordinates). You mentioned that you have information about the top-left corner and the resolution (0.01) - that's all we need. In the following I'm assuming that the lat/lon information about the top-left corner is saved in the the top_left_lat and top_left_lon variables. The extent needs to be passed in a tuple with a value for each of the edges (in the order left, right, bottom, top).

Hence, our extent can be calculated as follows:

extent_mat = (top_left_lon, top_left_lon + reshapeT.shape[1] * 0.01, top_left_lat - reshapeT.shape[0] * 0.01, top_left_lat)

Finally, we plot the matrix onto the same axes object, ax1, on which we already plotted the shape file to the calculated extent:

# Let's turn off autoscale first. This prevents
# the view of the plot to be limited to the image
# dimensions (instead of the entire shapefile). If you prefer
# that behaviour, just remove the following line
ax1.autoscale(False)

# Finally, let's plot!
ax1.imshow(reshapeT, extent=extent_mat)
bexi
  • 1,186
  • 5
  • 9