1

I'm trying to interpolate temperature data observed on an urban area formed by 5 locations. I am using cartopy to interpolate and draw the map, however, when I run the script the temperature interpolation is not shown and I only get the layer of the urban area with the color palette. Can someone help me fix this error? The link of shapefile is

https://www.dropbox.com/s/0u76k3yegvr09sx/LimiteAMG.shp?dl=0

https://www.dropbox.com/s/yxsmm3v2ey3ngsp/LimiteAMG.cpg?dl=0

https://www.dropbox.com/s/yx05n31dfkggbb6/LimiteAMG.dbf?dl=0

https://www.dropbox.com/s/a6nk0xczgjeen2d/LimiteAMG.prj?dl=0

https://www.dropbox.com/s/royw7s51n2f0a6x/LimiteAMG.qpj?dl=0

https://www.dropbox.com/s/7k44dcl1k5891qc/LimiteAMG.shx?dl=0

The Data is:

      Lat      Lon       tmax
  0   20.8208 -103.4434  22.8
  1   20.7019 -103.4728  17.7
  2   20.6833 -103.3500  24.9
  3   20.6280 -103.4261   NaN
  4   20.7205 -103.3172  26.4
  5   20.7355 -103.3782  25.7
  6   20.6593 -103.4136   NaN
  7   20.6740 -103.3842  25.8
  8   20.7585 -103.3904   NaN
  9   20.6230 -103.4265   NaN
 10   20.6209 -103.5004   NaN
 11   20.6758 -103.6439  24.5
 12   20.7084 -103.3901  24.0
 13   20.6353 -103.3994  23.0
 14   20.5994 -103.4133  25.0
 15   20.6302 -103.3422   NaN
 16   20.7400 -103.3122  23.0
 17   20.6061 -103.3475   NaN
 18   20.6400 -103.2900  23.0
 19   20.7248 -103.5305  24.0
 20   20.6238 -103.2401   NaN
 21   20.4753 -103.4451   NaN

the code is:

import cartopy
import cartopy.crs as ccrs
from matplotlib.colors import BoundaryNorm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import cartopy.io.shapereader as shpreader

from metpy.calc import get_wind_components
from metpy.cbook import get_test_data
from metpy.gridding.gridding_functions import interpolate,remove_nan_obervation 
from metpy.plots import add_metpy_logo
from metpy.units import units

to_proj = ccrs.PlateCarree()

data=pd.read_csv('/home/borisvladimir/Documentos/Datos/EMAs/EstacionesZMG/RedZMG.csv',usecols=(1,2,3),names=['Lat','Lon','tmax'],na_values=-99999,header=0)


fname='/home/borisvladimir/Dropbox/Diversos/Shapes/LimiteAMG.shp'
adm1_shapes = list(shpreader.Reader(fname).geometries())


lon = data['Lon'].values
lat = data['Lat'].values
xp, yp, _ = to_proj.transform_points(ccrs.Geodetic(), lon, lat).T


x_masked, y_masked, t = remove_nan_observations(xp, yp, data['tmax'].values)

tempx, tempy, temp = interpolate(x_masked, y_masked, t, interp_type='cressman', minimum_neighbors=3, search_radius=400000, hres=35000)
temp = np.ma.masked_where(np.isnan(temp), temp)

levels = list(range(-20, 20, 1))
cmap = plt.get_cmap('viridis')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

fig = plt.figure(figsize=(15, 10))
view = fig.add_subplot(1, 1, 1, projection=to_proj)

view.add_geometries(adm1_shapes, ccrs.PlateCarree(),edgecolor='black',  facecolor='white', alpha=0.5)


view.set_extent([-103.8, -103, 20.3, 21.099 ], ccrs.PlateCarree())

ZapLon,ZapLat=-103.50,20.80
GuadLon,GuadLat=-103.33,20.68
TonaLon,TonaLat=-103.21,20.62
TlaqLon,TlaqLat=-103.34,20.59
TlajoLon,TlajoLat=-103.44,20.47

plt.text(ZapLon,ZapLat,'Zapopan',transform=ccrs.Geodetic())
plt.text(GuadLon,GuadLat,'Guadalajara',transform=ccrs.Geodetic())
plt.text(TonaLon,TonaLat,'Tonala',transform=ccrs.Geodetic())
plt.text(TlaqLon,TlaqLat,'Tlaquepaque',transform=ccrs.Geodetic())
plt.text(TlajoLon,TlajoLat,'Tlajomulco',transform=ccrs.Geodetic())

mmb = view.pcolormesh(tempx, tempy,   temp,transform=ccrs.PlateCarree(),cmap=cmap, norm=norm)
plt.colorbar(mmb, shrink=.4, pad=0.02, boundaries=levels)
plt.show()

how to solve this? 

enter image description here

  • With only the .shp file, without .dbf ans .shx, we can do nothing – gene Feb 8 '18 at 16:09
  • Sorry. The links loaded now – user1345283 Feb 8 '18 at 18:51
  • Your problem is not relating to the shapefiles since you are able to plot it. Do concentrate on the pcolormesh only. – swatchai May 29 '18 at 2:23
1

With the input data modified and saved as lonlat_tmax.csv.

Lat,Lon,tmax
  0,   20.8208, -103.4434,  22.8
  1,   20.7019, -103.4728,  17.7
  2,   20.6833, -103.3500,  24.9
  3,   20.6280, -103.4261,NaN
  4,   20.7205, -103.3172,  26.4
  5,   20.7355, -103.3782,  25.7
  6,   20.6593, -103.4136,NaN
  7,   20.6740, -103.3842,  25.8
  8,   20.7585, -103.3904,NaN
  9,   20.6230, -103.4265,NaN
 10,   20.6209, -103.5004,NaN
 11,   20.6758, -103.6439,  24.5
 12,   20.7084, -103.3901,  24.0
 13,   20.6353, -103.3994,  23.0
 14,   20.5994, -103.4133,  25.0
 15,   20.6302, -103.3422,NaN
 16,   20.7400, -103.3122,  23.0
 17,   20.6061, -103.3475,NaN
 18,   20.6400, -103.2900,  23.0
 19,   20.7248, -103.5305,  24.0
 20,   20.6238, -103.2401,NaN
 21,   20.4753, -103.4451,NaN

Here is the working code that plot the filled contour and towns locations.

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from matplotlib.mlab import griddata
from matplotlib.colors import BoundaryNorm

# some settings
useproj = ccrs.PlateCarree()

# read data into pandas df
data = pd.read_csv('lonlat_tmax.csv', \
                   header=0, \
                   usecols=(1,2,3), \
                   names=['lat','lon','tmax'])

# flag data, select good/bad data for proper plotting
data['NaN_flag'] = data.tmax.isna()
good_data = data[ data.NaN_flag == False ]
bad_data = data[ data.NaN_flag == True ]

# extents of original data
xmin,xmax = data.lon.min(),data.lon.max()
ymin,ymax = data.lat.min(),data.lat.max()

# tmax limits
zmin, zmax = good_data.tmax.min(), good_data.tmax.max()

# create 2D meshgrid covering target area
numcols, numrows = 100, 100
xi = np.linspace(good_data.lon.min(), good_data.lon.max(), numcols)
yi = np.linspace(good_data.lat.min(), good_data.lat.max(), numrows)
xi, yi = np.meshgrid(xi, yi)

# interpolate for zi of all grid points
# data points, values: x, y, Z
# interpolate points: (xi, yi)
# (output) interpolated values: zi

# use `good_data` dataframe for interpolation
zi = griddata(good_data.lon.values, \
              good_data.lat.values, \
              good_data.tmax, \
              xi, yi, interp='linear')

fig = plt.figure(figsize=(10, 8))
rect = [0.05, 0.05, 0.95, 0.95]  # for map extent
ax = fig.add_axes( rect, projection=useproj )

# prep colormap and color norm
levels = list(range(int(zmin)-1, int(zmax)+1, 1))
cmap = plt.get_cmap('viridis')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

# contours plot
conf = ax.contourf(xi, yi, zi, zorder=4, alpha=0.95, cmap=cmap, norm=norm)
#ax.coastlines(linewidth=0.2, color='black') #nothing plotted in our case

# plot location names
ZapLon,ZapLat=-103.50,20.80
GuadLon,GuadLat=-103.33,20.68
TonaLon,TonaLat=-103.21,20.62
TlaqLon,TlaqLat=-103.34,20.59
TlajoLon,TlajoLat=-103.44,20.47

plt.text(ZapLon,ZapLat,'Zapopan', transform=ccrs.Geodetic(), zorder=9, fontsize=12)
plt.text(GuadLon,GuadLat,'Guadalajara', transform=ccrs.Geodetic(), zorder=9, fontsize=12)
plt.text(TonaLon,TonaLat,'Tonala', transform=ccrs.Geodetic(), zorder=9, fontsize=12)
plt.text(TlaqLon,TlaqLat,'Tlaquepaque', transform=ccrs.Geodetic(), zorder=9, fontsize=12)
plt.text(TlajoLon,TlajoLat,'Tlajomulco', transform=ccrs.Geodetic(), zorder=9, fontsize=12)

#cbar = plt.colorbar(conf, orientation='horizontal', fraction=.04, pad=0.05)
cbar = plt.colorbar(conf, shrink=.5, pad=0.02, boundaries=levels)
cbar.set_label("Max Temp - celcius")

# set limits to cover all data points
pad = .07
ax.set_xlim(xmin-pad, xmax+pad)
ax.set_ylim(ymin-pad, ymax+pad)

plt.show()

The resulting plot.

enter image description here

  • Could you edit and provide a quick "TL;DR" explanation of why this works compared to OP's approach? Otherwise anyone with the same problem will have to compare the code line by line! – Simbamangu Jun 2 '18 at 7:05
  • @Simbamangu That will be TL;DR. I will wait for the OP's response first. Then only when I have spare time. Thanks for the comment anyway. – swatchai Jun 2 '18 at 9:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.