Come posso ottenere il contorno del mio profilo sovrapposto su una mappa di base

Question

Questa è una domanda che ho chiesto diversi mesi fa e sto ancora lottando per arrivare a una soluzione. Il mio codice mi dà una mappa di base e un contorno di contorno affiancati (ma la stampa su file fornisce solo la trama del contorno), ma li voglio sovrapposti. La soluzione migliore sarebbe quella qui https://gist.github.com/oblakeobjet/7546272 ma questo non mostra come introdurre i dati, ed è difficile quando si sta imparando da zero online. Senza stancare persone molto gentili, spero che la soluzione sia semplice come cambiare una linea di codice e che qualcuno possa aiutare. Il mio codice

#!/usr/bin/python 
# vim: set fileencoding=UTF8 

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

#fig = plt.figure(figsize=(10,8)) #when uncommented draws map with colorbar but no contours 

#prepare a basemap 

m = Basemap(projection = 'merc',llcrnrlon = 21, llcrnrlat = -18, urcrnrlon = 34, urcrnrlat = -8, resolution='h') 

# draw country outlines. 

m.drawcountries(linewidth=0.5, linestyle='solid', color='k', antialiased=1, ax=None, zorder=None) 
m.drawmapboundary(fill_color = 'white') 
m.fillcontinents(color='coral',lake_color='blue') 
parallels = np.arange(-18, -8, 2.) 
m.drawparallels(np.arange(-18, -8, 2.), color = 'black', linewidth = 0.5) 
m.drawparallels(parallels,labels=[True,False,False,False]) 
meridians = np.arange(22,34, 2.) 
m.drawmeridians(np.arange(21,36, 2.), color = '0.25', linewidth = 0.5) 
m.drawmeridians(meridians,labels=[False,False,False,True]) 

fig = plt.figure(figsize=(10,8))  # At this position or commented draws teo figures side by side 

#-- Read the data. 

data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True) 

#-- Now gridding data. First making a regular grid to interpolate onto 

numcols, numrows = 300, 300 
xi = np.linspace(data.Lon.min(), data.Lon.max(), numcols) 
yi = np.linspace(data.Lat.min(), data.Lat.max(), numrows) 
xi, yi = np.meshgrid(xi, yi) 

#-- Interpolating at the points in xi, yi 

x, y, z = data.Lon.values, data.Lat.values, data.Z.values 
zi = griddata(x, y, z, xi, yi) 

#-- Display and write the results 

m = plt.contourf(xi, yi, zi) 
plt.scatter(data.Lon, data.Lat, c=data.Z, s=100, 
     vmin=zi.min(), vmax=zi.max()) 
fig.colorbar(m) 
plt.savefig("rainfall.jpg", format="jpg") 

Le trame che ricevo questo aspetto e

e miei dati

32.6 -13.6 41 
27.1 -16.9 43 
32.7 -10.2 46 
24.2 -13.6 33 
28.5 -14.4 43 
28.1 -12.6 33 
27.9 -15.8 46 
24.8 -14.8 44 
31.1 -10.2 35 
25.9 -13.5 24 
29.1 -9.8 10 
25.8 -17.8 39 
33.2 -12.3 44 
28.3 -15.4 46 
27.6 -16.1 47 
28.9 -11.1 31 
31.3 -8.9 39 
31.9 -13.3 45 
23.1 -15.3 31 
31.4 -11.9 39 
27.1 -15.0 42 
24.4 -11.8 15 
28.6 -13.0 39 
31.3 -14.3 44 
23.3 -16.1 39 
30.2 -13.2 38 
24.3 -17.5 32 
26.4 -12.2 23 
23.1 -13.5 27 

Answer 1

Sei quasi arrivato, ma la mappa di base può essere capricciosa e devi gestire l'ordine z di grafici/dettagli mappa. Inoltre, devi trasformare le tue coordinate lon/lat in coordinate di proiezione mappa prima di tracciarle usando la mappa di base.

Ecco una soluzione completa, che fornisce il seguente output. Ho modificato alcuni colori e linee di fuga per rendere tutto più leggibile, YMMV. Ho anche ridimensionato le dimensioni dei punti di dispersione per il valore "medio" normalizzato (data['Z']) - puoi semplicemente rimuoverlo e sostituirlo ad es. 50 se preferisci una dimensione costante (sembrerà il marcatore più grande).

Si prega di specificare anche le unità di pioggia, e la durata della misura che ha generato i mezzi, se possibile:

import numpy as np 
import pandas as pd 
from matplotlib.mlab import griddata 
from mpl_toolkits.basemap import Basemap 
import matplotlib.pyplot as plt 
from matplotlib.colors import Normalize 
%matplotlib inline 

# set up plot 
plt.clf() 
fig = plt.figure(figsize=(10, 8)) 
ax = fig.add_subplot(111, axisbg='w', frame_on=False) 

# grab data 
data = pd.read_csv('../../data/meansr.txt', delim_whitespace=True) 
norm = Normalize() 

# define map extent 
lllon = 21 
lllat = -18 
urlon = 34 
urlat = -8 

# Set up Basemap instance 
m = Basemap(
    projection = 'merc', 
    llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat, 
    resolution='h') 

# transform lon/lat coordinates to map projection 
data['projected_lon'], data['projected_lat'] = m(*(data.Lon.values, data.Lat.values)) 

# grid data 
numcols, numrows = 1000, 1000 
xi = np.linspace(data['projected_lon'].min(), data['projected_lon'].max(), numcols) 
yi = np.linspace(data['projected_lat'].min(), data['projected_lat'].max(), numrows) 
xi, yi = np.meshgrid(xi, yi) 

# interpolate 
x, y, z = data['projected_lon'].values, data['projected_lat'].values, data.Z.values 
zi = griddata(x, y, z, xi, yi) 

# draw map details 
m.drawmapboundary(fill_color = 'white') 
m.fillcontinents(color='#C0C0C0', lake_color='#7093DB') 
m.drawcountries(
    linewidth=.75, linestyle='solid', color='#000073', 
    antialiased=True, 
    ax=ax, zorder=3) 

m.drawparallels(
    np.arange(lllat, urlat, 2.), 
    color = 'black', linewidth = 0.5, 
    labels=[True, False, False, False]) 
m.drawmeridians(
    np.arange(lllon, urlon, 2.), 
    color = '0.25', linewidth = 0.5, 
    labels=[False, False, False, True]) 

# contour plot 
con = m.contourf(xi, yi, zi, zorder=4, alpha=0.6, cmap='RdPu') 
# scatter plot 
m.scatter(
    data['projected_lon'], 
    data['projected_lat'], 
    color='#545454', 
    edgecolor='#ffffff', 
    alpha=.75, 
    s=50 * norm(data['Z']), 
    cmap='RdPu', 
    ax=ax, 
    vmin=zi.min(), vmax=zi.max(), zorder=4) 

# add colour bar and title 
# add colour bar, title, and scale 
cbar = plt.colorbar(conf, orientation='horizontal', fraction=.057, pad=0.05) 
cbar.set_label("Mean Rainfall - mm") 

m.drawmapscale(
    24., -9., 28., -13, 
    100, 
    units='km', fontsize=10, 
    yoffset=None, 
    barstyle='fancy', labelstyle='simple', 
    fillcolor1='w', fillcolor2='#000000', 
    fontcolor='#000000', 
    zorder=5) 

plt.title("Mean Rainfall") 
plt.savefig("rainfall.png", format="png", dpi=300, transparent=True) 
plt.show() 

usando il metodo di matplotlib griddata è conveniente, ma può anche essere lento. In alternativa, è possibile utilizzare griddata metodi di SciPy, che sono sia più veloce e più flessibile:

from scipy.interpolate import griddata as gd 

zi = gd(
    (data[['projected_lon', 'projected_lat']]), 
    data.Z.values, 
    (xi, yi), 
    method='linear') 

Se si utilizza il metodo di SciPy griddata, dovrete anche per determinare quale dei metodi (nearest, linear, cubic) forniscono la trama risultante migliore.

Devo aggiungere che i metodi di interpolazione dimostrati e discussi sopra sono i più semplici possibili e non sono necessariamente validi per l'interpolazione dei dati sulle precipitazioni. This article offre una buona panoramica degli approcci e delle considerazioni validi per l'uso in idrologia e modellistica idrologica. L'implementazione di questi (probabilmente usando Scipy) è lasciata come esercizio & c.

Answer 2

non ho tutto installato qui per eseguire il codice, ma si dovrebbe provare la stampa alla mappa di base m creata come:

# fig = plt.figure(figsize=(10,8)) # omit this at line 28 

(...) 

m.contourf(xi, yi, zi) 
m.scatter(data.Lon, data.Lat, c=data.Z, s=100, 
    vmin=zi.min(), vmax=zi.max()) 

(informi se questo non funziona)