Как получить данные контурного графика для каждой точки разброса?

Question

Я построил контурный участок как фон, который представляет собой высоту области.

И были установлены 100 точек рассеяния, представляющих реальный источник выбросов загрязняющих веществ. Есть ли способ получить высоту каждой точки?

Это мой код:

%matplotlib inline 

fig=plt.figure(figsize=(16,16)) 
ax=plt.subplot() 
xi,yi = np.linspace(195.2260,391.2260,50),   
np.linspace(4108.9341,4304.9341,50) 
height=np.array(list(csv.reader(open("/Users/HYF/Documents/SJZ_vis/Concentration/work/terr_grd.csv","rb"),delimiter=','))).astype('float') 
cmap = cm.get_cmap(name='terrain', lut=None) 
terrf = plt.contourf(xi, yi, height,100, cmap=cmap) 
terr = plt.contour(xi, yi, height, 100, 
      colors='k',alpha=0.5 
      ) 
plt.clabel(terr, fontsize=7, inline=20) 
ax.autoscale(False) 
point= plt.scatter(dat_so2["xp"], dat_so2["yp"], marker='o',c="grey",s=40) 


ax.autoscale(False) 
for i in range(0,len(dat_so2["xp"]),1): 
    plt.text(dat_so2["xp"][i], dat_so2["yp"][i],  
    str(i),color="White",fontsize=16) 

ax.set_xlim(225,275) 
ax.set_ylim(4200,4260) 

plt.show() 

Answer 1

Вы можете сделать это с scipy.interpolate.interp2d

Например, вы могли бы добавить в код:

from scipy import interpolate 
hfunc = interpolate.interp2d(xi,yi,height) 

pointheights = np.zeros(dat_so2["xp"].shape) 
for i,(x,y) in enumerate(zip(dat_so2["xp"],dat_so2["yp"])): 
    pointheights[i]=hfunc(x,y) 

Сведя вместе с остальной частью вашего скрипта, и некоторый образец da та, дает это (я упростил пару вещей, но вы получите идею):

import matplotlib.pyplot as plt 
import matplotlib.cm as cm 
import numpy as np 
from scipy import interpolate 

fig=plt.figure(figsize=(8,8)) 
ax=plt.subplot() 
#xi,yi = np.linspace(195.2260,391.2260,50),np.linspace(4108.9341,4304.9341,50) 
xi,yi = np.linspace(225,275,50),np.linspace(4200,4260,50) 

# A made up function of height (in place of your data) 
XI,YI = np.meshgrid(xi,yi) 
height = (XI-230.)**2 + (YI-4220.)**2 
#height=np.array(list(csv.reader(open("/Users/HYF/Documents/SJZ_vis/Concentration/work/terr_grd.csv","rb"),delimiter=','))).astype('float') 

cmap = cm.get_cmap(name='terrain', lut=None) 
terrf = plt.contourf(xi, yi, height,10, cmap=cmap) 
terr = plt.contour(xi, yi, height, 10, 
      colors='k',alpha=0.5 
      ) 
plt.clabel(terr, fontsize=7, inline=20) 
ax.autoscale(False) 

# Some made up sample points 
dat_so2 = np.array([(230,4210),(240,4220),(250,4230),(260,4240),(270,4250)],dtype=[("xp","f4"),("yp","f4")]) 

point= plt.scatter(dat_so2["xp"], dat_so2["yp"], marker='o',c="grey",s=40) 

# The interpolation function 
hfunc = interpolate.interp2d(xi,yi,height) 

# Now, for each point, lets interpolate the height 
pointheights = np.zeros(dat_so2["xp"].shape) 
for i,(x,y) in enumerate(zip(dat_so2["xp"],dat_so2["yp"])): 
    pointheights[i]=hfunc(x,y) 
print pointheights 


ax.autoscale(False) 
for i in range(0,len(dat_so2["xp"]),1): 
    plt.text(dat_so2["xp"][i], dat_so2["yp"][i],  
      str(i),color="White",fontsize=16) 
    # We can also add a height label to the plot 
    plt.text(dat_so2["xp"][i], dat_so2["yp"][i],  
      "{:4.1f}".format(pointheights[i]),color="black",fontsize=16,ha='right',va='top') 


ax.set_xlim(225,275) 
ax.set_ylim(4200,4260) 

plt.show()