Python - 从绘制的地图中使用 'onclick(event)' 抓取网格数据
Python - Grab grid data with 'onclick(event)' from plotted map
对于缺少可重现的示例提前表示歉意。我的数据结构方式使这变得困难。
我正在构建一个程序来帮助我收集数据（在本例中为网格坐标）以训练分类算法。我希望能够从我点击的情节中拉出特定的网格单元格。到目前为止，我已经使用 fig.canvas.mpl_connect('button_press_event', onclick) 并且能够使用 event.x 和 event.y 提取像素坐标（我也尝试过使用 event.xdata 和 event.ydata 但这似乎给了我一些我不认识的其他坐标集？）但我一直无法提取我想要的实际数据。
在我的代码中，我使用了以下模块：
from netCDF4 import Dataset # reads netCDF file
from os import listdir
from os.path import isfile, join
import matplotlib.pyplot as plt
import pylab
from mpl_toolkits.basemap import Basemap # basemap tools
from datetime import datetime, timedelta # for working with datetimes
from random import randint
import numpy as np
我的数据来自存储在 NetCDF 文件中的数据，该文件是更大模型 运行 的子集，变量存储为 165*116 数组（这是网格维度）。为了构建我的情节，我使用以下函数：
def grab_sst_time(time_idx):
"""
This reads the datetime value in my NetCDF files into local time
"""
dtcon_days = time[time_idx]
dtcon_start = datetime(1990,1,1) # This is the "days since" part
dtcon_delta = timedelta(dtcon_days/24/60/60) # Create a time delta object from the number of days
dtcon_offset = dtcon_start + dtcon_delta # Add the specified number of days to 1990
frame_time = dtcon_offset
return frame_time
def plot_temp(temp, time_idx, fig_no):
"""
Plot temperature
"""
# make frame_idx an integer to avoid slicing errors
frame_idx = int(time_idx)
# get 'frame_time'
frame_time = grab_sst_time(frame_idx)
# map setup
fig = plt.figure()
fig.subplots_adjust(left=0., right=1., bottom=0., top=0.9)
# Setup the map
m = Basemap(projection='merc', llcrnrlat=-38.050653, urcrnrlat=-34.453367,\
llcrnrlon=147.996456, urcrnrlon=152.457344, lat_ts=20, resolution='h')
# draw stuff
m.drawcoastlines()
m.fillcontinents(color='black')
# plot temp
cs = m.pcolor(lons,lats,np.squeeze(temp), latlon = True ,vmin=temp_min, vmax=temp_max, cmap='plasma')
# plot colourbar
plt.colorbar()
# datetime title
plt.title('Regional - Temperature (Celcius)\n' +
frame_time.strftime("%Y-%m-%d %H:%M:%S") + ' | ' + str(fname) + '_idx: ' + str(frame_idx))
# stop axis from being cropped
plt.tight_layout()
return fig
当使用正确的数据调用时，会生成如下图：
现在我开始感到困惑了。我从一个随机文件沙时间制作我的图，并尝试制作它们 return 我点击的数据并将其保存到数据库中。下面是调用上述函数的代码：
def onclick(event):
"""
On click function for selecting data
"""
global ix, iy
ix, iy = event.xdata, event.ydata
print('x='+str(ix), 'y='+str(iy), 'class='+str(class_value), time_value)
global train_data
train_data.append((ix, iy, class, time_value))
# set colour scale variables
temp_min = 14
temp_max = 24
# list to hold data collected
train_data = []
# __Random Data__
# get list of files in data directory
directory = "/Users/directory/with/my/data"
file_ls = [f for f in listdir(directory) if isfile(join(directory, f))]
file_ls = list(filter(lambda x:'naroom_avg' in x, file_ls))
# set randomness seed
plot_num = 2
np.random.seed(1010)
rnd_file = np.random.randint(len(file_ls), size=plot_num)
rnd_times = np.random.randint(29, size=plot_num)
# __Make plots__
for i in range(0, plot_num):
# grab file
file_no = rnd_file[i]
file_path = directory + "/" + file_ls[file_no]
fname = str(file_ls[i])[11:16]
# grab time
time_idx = rnd_times[i]
fh = Dataset(file_path, mode='r')
# extract data
lats = fh.variables['lat_rho'][:]
lons = fh.variables['lon_rho'][:]
time = fh.variables['ocean_time'][:]
temp = fh.variables['temp'][time_idx,29,:,:]
# time output
time_value = grab_sst_time(time_idx)
# make interactive plot
fig = plot_temp(temp, time_idx, i)
cid = fig.canvas.mpl_connect('button_press_event', onclick)
# Select class A
class = 'A'
print('Select data...')
input("Press Enter to continue...")
plt.show()
这似乎工作正常，但我无法提取我真正想要的数据（NetCDF 数据）。理想情况下，我想要 return 网格单元（即 165*116 网格的坐标，类似于 x=112，y=154）。相反，单击会给我如下输出：
x=432047.683555 y=449210.017634 class=A 2004-06-17 12:00:00
x=448214.625063 y=430733.513053 class=A 2004-06-17 12:00:00
x=448214.625063 y=408792.663863 class=A 2004-06-17 12:00:00
x=448792.015832 y=441703.937648 class=A 2004-06-17 12:00:00
知道如何获得我想要的输出吗？
补充资料
以下是我的 NetCDF 文件的相关变量的结构：
dimensions:
s_rho = 30;
eta_rho = 116;
xi_rho = 165;
ocean_time = UNLIMITED; // (30 currently)
variables:
double lat_rho(eta_rho=116, xi_rho=165);
:long_name = "latitude of RHO-points";
:units = "degree_north";
:standard_name = "latitude";
:field = "lat_rho, scalar";
:_CoordinateAxisType = "Lat";
double lon_rho(eta_rho=116, xi_rho=165);
:long_name = "longitude of RHO-points";
:units = "degree_east";
:standard_name = "longitude";
:field = "lon_rho, scalar";
:_CoordinateAxisType = "Lon";
double ocean_time(ocean_time=30);
:long_name = "averaged time since initialization";
:units = "seconds since 1990-01-01 00:00:00";
:calendar = "gregorian";
:field = "time, scalar, series";
:_CoordinateAxisType = "Time";
float temp(ocean_time=30, s_rho=30, eta_rho=116, xi_rho=165);
:long_name = "time-averaged potential temperature";
:units = "Celsius";
:time = "ocean_time";
:coordinates = "lon_rho lat_rho s_rho ocean_time";
:field = "temperature, scalar, series";
:_FillValue = 1.0E37f; // float
要从经度、纬度转换为投影坐标，您可以使用 lon,lat 值调用底图实例，
m = Basemap(...)
xpt,ypt = m(lon,lat)
要将投影坐标转换回经纬度坐标，请使用逆函数
lon,lat = m(xpt,ypt, inverse=True)
我们可能会修改 here 中的示例，以在屏幕上显示点击的坐标。
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# setup Lambert Conformal basemap.
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.drawmapboundary(fill_color='azure')
m.fillcontinents(color='sandybrown',lake_color='azure')
parallels = np.arange(0.,81,10.)
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
m.drawmeridians(meridians,labels=[True,False,False,True])
# plot blue dot on Boulder, colorado and label it as such.
lon, lat = -104.237, 40.125 # Location of Boulder
# convert to map projection coords.
# Note that lon,lat can be scalars, lists or numpy arrays.
xpt,ypt = m(lon,lat)
# convert back to lat/lon
lonpt, latpt = m(xpt,ypt,inverse=True)
point, = m.plot(xpt,ypt,'bo') # plot a blue dot there
# put some text next to the dot, offset a little bit
# (the offset is in map projection coordinates)
annotation = plt.annotate('Boulder (%5.1fW,%3.1fN)' % (lon, lat), xy=(xpt,ypt),
xytext=(20,35), textcoords="offset points",
bbox={"facecolor":"w", "alpha":0.5},
arrowprops=dict(arrowstyle="->", connectionstyle="arc3"))
def onclick(event):
ix, iy = event.xdata, event.ydata
xpti, ypti = m(ix, iy,inverse=True)
string = '(%5.1fW,%3.1fN)' % (xpti, ypti)
print(string)
annotation.xy = (ix, iy)
point.set_data([ix], [iy])
annotation.set_text(string)
plt.gcf().canvas.draw_idle()
cid = plt.gcf().canvas.mpl_connect("button_press_event", onclick)
plt.show()
对于缺少可重现的示例提前表示歉意。我的数据结构方式使这变得困难。
我正在构建一个程序来帮助我收集数据（在本例中为网格坐标）以训练分类算法。我希望能够从我点击的情节中拉出特定的网格单元格。到目前为止，我已经使用 fig.canvas.mpl_connect('button_press_event', onclick) 并且能够使用 event.x 和 event.y 提取像素坐标（我也尝试过使用 event.xdata 和 event.ydata 但这似乎给了我一些我不认识的其他坐标集？）但我一直无法提取我想要的实际数据。
在我的代码中，我使用了以下模块：
from netCDF4 import Dataset # reads netCDF file
from os import listdir
from os.path import isfile, join
import matplotlib.pyplot as plt
import pylab
from mpl_toolkits.basemap import Basemap # basemap tools
from datetime import datetime, timedelta # for working with datetimes
from random import randint
import numpy as np
我的数据来自存储在 NetCDF 文件中的数据，该文件是更大模型 运行 的子集，变量存储为 165*116 数组（这是网格维度）。为了构建我的情节，我使用以下函数：
def grab_sst_time(time_idx):
"""
This reads the datetime value in my NetCDF files into local time
"""
dtcon_days = time[time_idx]
dtcon_start = datetime(1990,1,1) # This is the "days since" part
dtcon_delta = timedelta(dtcon_days/24/60/60) # Create a time delta object from the number of days
dtcon_offset = dtcon_start + dtcon_delta # Add the specified number of days to 1990
frame_time = dtcon_offset
return frame_time
def plot_temp(temp, time_idx, fig_no):
"""
Plot temperature
"""
# make frame_idx an integer to avoid slicing errors
frame_idx = int(time_idx)
# get 'frame_time'
frame_time = grab_sst_time(frame_idx)
# map setup
fig = plt.figure()
fig.subplots_adjust(left=0., right=1., bottom=0., top=0.9)
# Setup the map
m = Basemap(projection='merc', llcrnrlat=-38.050653, urcrnrlat=-34.453367,\
llcrnrlon=147.996456, urcrnrlon=152.457344, lat_ts=20, resolution='h')
# draw stuff
m.drawcoastlines()
m.fillcontinents(color='black')
# plot temp
cs = m.pcolor(lons,lats,np.squeeze(temp), latlon = True ,vmin=temp_min, vmax=temp_max, cmap='plasma')
# plot colourbar
plt.colorbar()
# datetime title
plt.title('Regional - Temperature (Celcius)\n' +
frame_time.strftime("%Y-%m-%d %H:%M:%S") + ' | ' + str(fname) + '_idx: ' + str(frame_idx))
# stop axis from being cropped
plt.tight_layout()
return fig
当使用正确的数据调用时，会生成如下图：
现在我开始感到困惑了。我从一个随机文件沙时间制作我的图，并尝试制作它们 return 我点击的数据并将其保存到数据库中。下面是调用上述函数的代码：
def onclick(event):
"""
On click function for selecting data
"""
global ix, iy
ix, iy = event.xdata, event.ydata
print('x='+str(ix), 'y='+str(iy), 'class='+str(class_value), time_value)
global train_data
train_data.append((ix, iy, class, time_value))
# set colour scale variables
temp_min = 14
temp_max = 24
# list to hold data collected
train_data = []
# __Random Data__
# get list of files in data directory
directory = "/Users/directory/with/my/data"
file_ls = [f for f in listdir(directory) if isfile(join(directory, f))]
file_ls = list(filter(lambda x:'naroom_avg' in x, file_ls))
# set randomness seed
plot_num = 2
np.random.seed(1010)
rnd_file = np.random.randint(len(file_ls), size=plot_num)
rnd_times = np.random.randint(29, size=plot_num)
# __Make plots__
for i in range(0, plot_num):
# grab file
file_no = rnd_file[i]
file_path = directory + "/" + file_ls[file_no]
fname = str(file_ls[i])[11:16]
# grab time
time_idx = rnd_times[i]
fh = Dataset(file_path, mode='r')
# extract data
lats = fh.variables['lat_rho'][:]
lons = fh.variables['lon_rho'][:]
time = fh.variables['ocean_time'][:]
temp = fh.variables['temp'][time_idx,29,:,:]
# time output
time_value = grab_sst_time(time_idx)
# make interactive plot
fig = plot_temp(temp, time_idx, i)
cid = fig.canvas.mpl_connect('button_press_event', onclick)
# Select class A
class = 'A'
print('Select data...')
input("Press Enter to continue...")
plt.show()
这似乎工作正常，但我无法提取我真正想要的数据（NetCDF 数据）。理想情况下，我想要 return 网格单元（即 165*116 网格的坐标，类似于 x=112，y=154）。相反，单击会给我如下输出：
x=432047.683555 y=449210.017634 class=A 2004-06-17 12:00:00
x=448214.625063 y=430733.513053 class=A 2004-06-17 12:00:00
x=448214.625063 y=408792.663863 class=A 2004-06-17 12:00:00
x=448792.015832 y=441703.937648 class=A 2004-06-17 12:00:00
知道如何获得我想要的输出吗？
补充资料
以下是我的 NetCDF 文件的相关变量的结构：
dimensions:
s_rho = 30;
eta_rho = 116;
xi_rho = 165;
ocean_time = UNLIMITED; // (30 currently)
variables:
double lat_rho(eta_rho=116, xi_rho=165);
:long_name = "latitude of RHO-points";
:units = "degree_north";
:standard_name = "latitude";
:field = "lat_rho, scalar";
:_CoordinateAxisType = "Lat";
double lon_rho(eta_rho=116, xi_rho=165);
:long_name = "longitude of RHO-points";
:units = "degree_east";
:standard_name = "longitude";
:field = "lon_rho, scalar";
:_CoordinateAxisType = "Lon";
double ocean_time(ocean_time=30);
:long_name = "averaged time since initialization";
:units = "seconds since 1990-01-01 00:00:00";
:calendar = "gregorian";
:field = "time, scalar, series";
:_CoordinateAxisType = "Time";
float temp(ocean_time=30, s_rho=30, eta_rho=116, xi_rho=165);
:long_name = "time-averaged potential temperature";
:units = "Celsius";
:time = "ocean_time";
:coordinates = "lon_rho lat_rho s_rho ocean_time";
:field = "temperature, scalar, series";
:_FillValue = 1.0E37f; // float
要从经度、纬度转换为投影坐标，您可以使用 lon,lat 值调用底图实例，
m = Basemap(...)
xpt,ypt = m(lon,lat)
要将投影坐标转换回经纬度坐标，请使用逆函数
lon,lat = m(xpt,ypt, inverse=True)
我们可能会修改 here 中的示例，以在屏幕上显示点击的坐标。
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# setup Lambert Conformal basemap.
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
m.drawmapboundary(fill_color='azure')
m.fillcontinents(color='sandybrown',lake_color='azure')
parallels = np.arange(0.,81,10.)
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
m.drawmeridians(meridians,labels=[True,False,False,True])
# plot blue dot on Boulder, colorado and label it as such.
lon, lat = -104.237, 40.125 # Location of Boulder
# convert to map projection coords.
# Note that lon,lat can be scalars, lists or numpy arrays.
xpt,ypt = m(lon,lat)
# convert back to lat/lon
lonpt, latpt = m(xpt,ypt,inverse=True)
point, = m.plot(xpt,ypt,'bo') # plot a blue dot there
# put some text next to the dot, offset a little bit
# (the offset is in map projection coordinates)
annotation = plt.annotate('Boulder (%5.1fW,%3.1fN)' % (lon, lat), xy=(xpt,ypt),
xytext=(20,35), textcoords="offset points",
bbox={"facecolor":"w", "alpha":0.5},
arrowprops=dict(arrowstyle="->", connectionstyle="arc3"))
def onclick(event):
ix, iy = event.xdata, event.ydata
xpti, ypti = m(ix, iy,inverse=True)
string = '(%5.1fW,%3.1fN)' % (xpti, ypti)
print(string)
annotation.xy = (ix, iy)
point.set_data([ix], [iy])
annotation.set_text(string)
plt.gcf().canvas.draw_idle()
cid = plt.gcf().canvas.mpl_connect("button_press_event", onclick)
plt.show()