2

I have a polygon shapefile (the state of Illinois) and a CSV file with (lat, lon, zvalue). I want to plot a smooth contour plot representing those zvalues. Following is my code: 

import glob
import fiona
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.mlab import griddata


# Read in the tabulated data
tabfname = glob.glob("Outputs\\*.csv")[0]
df = pd.read_table(tabfname, sep=",")
print(df.head())
lat, lon, z = list(df.y), list(df.x), list(df["Theil Sen Slope"])
z0, z1, z2 = np.min(z)+0.03, np.mean(z), np.max(z)-0.01


# Read some metadata of the shapefile
shp = glob.glob("GIS\\*.shp")[0]
with fiona.drivers():
    with fiona.open(shp) as src:
        bnds = src.bounds        
extent = [values for values in bnds]
lono = np.mean([extent[0], extent[2]])     
lato = np.mean([extent[1], extent[3]])
llcrnrlon = extent[0]-0.5
llcrnrlat = extent[1]-0.5
urcrnrlon = extent[2]+0.5
urcrnrlat = extent[3]+0.5


# Create a Basemap
fig = plt.figure()
ax = fig.add_subplot(111)

m = Basemap(llcrnrlon=llcrnrlon, llcrnrlat=llcrnrlat,
            urcrnrlon=urcrnrlon, urcrnrlat=urcrnrlat,
             resolution='i', projection='tmerc' , lat_0 = lato, lon_0 = lono)
# Read in and display the shapefile
m.readshapefile(shp.split(".")[0], 'shf', zorder=2, drawbounds=True)


# Compute the number of bins to aggregate data
nx = 100
ny = 100

# Create a mesh and interpolate data
xi = np.linspace(llcrnrlon, urcrnrlon, nx)
yi = np.linspace(llcrnrlat, urcrnrlat, ny)
xgrid, ygrid = np.meshgrid(xi, yi)
xs, ys = m(xgrid, ygrid) 
zs = griddata(lon, lat, z, xgrid, ygrid, interp='nn')

# Plot the contour map
conf = m.contourf(xs, ys, zs, 30, zorder=1, cmap='jet')
cbar = m.colorbar(conf, location='bottom', pad="5%", ticks=(z0, z1, z2))

# Scatter plot of the points that make up the contour
for x, y in zip(lon, lat):
    X, Y = m(x,y)
    m.scatter(X, Y, zorder=4, color='black', s=1)

plt.show()
fig.savefig("Myplot.png", format="png")

And this is the output I got(The scattered black dots are there to show the spatial distribution of the points from which the interpolation was generated. I used Nearest Neighbor interpolation method here.): enter image description here

I basically referred to the examples given in the following two links to plot this:

Now this image has 3 problems:

  1. The interpolated contour does not expand within the whole of the shapefile
  2. The part of the contour plot protruding out of the shapefile boundary is not masked off
  3. The contour is not smooth.

What I want is to overcome these three deficiencies of my plot and generate a smooth and nice looking plot similar to the ones shown below (Source: https://doi.org/10.1175/JCLI3557.1 ): enter image description here Source:https://doi.org/10.1175/JCLI3557.1

How do I achieve that?

SereneWizard
  • 166
  • 3
  • 14
  • 1
    About removing the overlapping parts of your contour plot, have a look if [this](https://stackoverflow.com/a/48062299/2454357) helps you. – Thomas Kühn Feb 02 '18 at 07:16
  • 2
    [This question](https://stackoverflow.com/questions/42426095/matplotlib-contour-contourf-of-concave-non-gridded-data) may be of interest here. – ImportanceOfBeingErnest Feb 02 '18 at 09:02
  • @ImportanceOfBeingErnest Thank you for the link here. But my problem is still not solved. As you can see, the interpolation does not go outside of the outer bounds of the scatter points. Because of that, I cannot get the contour to cover the whole of the shapefile. Do you suggest anything to solve this issue? – SereneWizard Feb 02 '18 at 17:44
  • 1
    Those are really two completely different things. Cropping/clipping may be rather easy, here you may even clip the contour plot by the shape path itself. But to extent the contour to the edges of the shape, there cannot be a general solution, because *you* need to know what data you want to put there if it isn't present anywhere in the code. – ImportanceOfBeingErnest Feb 02 '18 at 17:48
  • By the way, do you have any idea what did these guys do in those two US maps I shared? Looks like it’s not interpolation and contouring. And I, for my life, cannot figure out what process they followed. – SereneWizard Feb 02 '18 at 18:02
  • To me it looks exactly like a contour plot, just with much more data that is much more dense and extents up to the country borders. – ImportanceOfBeingErnest Feb 02 '18 at 18:10
  • @ImportanceOfBeingErnest If you look at the paper from which I borrowed the images, in Figure 1 we see the number and distribution of the points they used. And these points do not necessarily extend into the boundary of the shapefile. But still they got that end product. Do you have any suggestion about what path would be the best for me to follow, with the resources I have, to achieve the desired output? – SereneWizard Feb 02 '18 at 18:19
  • The paper tells you that they used ESRI ARCGIS for visualization of the previously smoothened data via a head-banging algorithm. – ImportanceOfBeingErnest Feb 02 '18 at 18:24
  • @ImportanceOfBeingErnest That is not the path I want to walk into. I guess, I am ok with rough interpolation. What do you suggest I should oo to extend the interpolation to the boundary of the shapefile? I am ok with it going beyond the bounds of the shapefile either, since now I have solved the masking issue. – SereneWizard Feb 02 '18 at 18:30
  • Did you find the solution to your proble @SereneWizard? I would be interested... Thanks! – leroygr Feb 07 '19 at 11:18

0 Answers0