【问题标题】:Facing weird problem when trying to plot simple lat/ lon points尝试绘制简单的纬度/经度点时遇到奇怪的问题
【发布时间】:2019-08-21 06:12:26
【问题描述】:

我有以下数据框(相应的 csv 托管在这里:http://www.sharecsv.com/s/3795d862c1973efa311d8a770e978215/t.csv

            lat     lon
count   6159.000000     6159.000000
mean    37.764859   -122.355491
std     0.028214    0.038874
min     37.742200   -122.482783
25%     37.746317   -122.360133
50%     37.746417   -122.333717
75%     37.785825   -122.331300
max     37.818133   -122.331167

以下代码图正确:

    test_ax = plt.axes(projection=ccrs.Mercator())
    test_ax.plot(test_df['lon'], test_df['lat'], color="blue", linewidth=4, alpha=1.0,
            transform=ccrs.Geodetic())
    plt.show()

但如果我取一个子集,它不会:

    test_ax = plt.axes(projection=ccrs.Mercator())
    test_ax.plot(test_df['lon'][:1001], test_df['lat'][:1001], color="blue", linewidth=4, alpha=1.0,
            transform=ccrs.Geodetic())
    plt.show()

但对另一个子集这样做。

    test_ax = plt.axes(projection=ccrs.Mercator())
    test_ax.plot(test_df['lon'][:3501], test_df['lat'][:3501], color="blue", linewidth=4, alpha=1.0,
            transform=ccrs.Geodetic())
    plt.show()

我很确定我在做一些愚蠢的事情,但我无法弄清楚这种行为的原因。

编辑:

在进一步的实验中,我发现如果我手动设置地图范围以包含0 子午线,则子集:1001 的图(之前未显示)开始显示(旧金山附近的蓝点)。

    test_ax = plt.axes(projection=ccrs.Mercator())
    test_ax.plot(test_df['lon'][:1001], test_df['lat'][:1001], color="blue", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
    test_ax.coastlines()
    test_ax.set_extent([-130, 0, 30, 40], crs=ccrs.Geodetic())
    test_ax.gridlines(draw_labels=True)
    plt.show()

编辑:带有可重现的示例

(对于 jupyter 笔记本)

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pandas as pd

df_csv_url = 'http://www.sharecsv.com/dl/76dd767525a37180ca54cd1d9314b9dc/t1.csv'
test_df = pd.read_csv(df_csv_url)
figure_params = { 'width': 9.6, 'height': 5.4 }

fig = plt.figure(
        figsize=(figure_params["width"], figure_params["height"])        
    )
test_ax = fig.add_axes((0, 0.5, 0.5, 0.5), projection=ccrs.Mercator(), label="map1")
test_ax.plot(test_df['lon'], test_df['lat'], color="blue", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax.coastlines()
test_ax.gridlines(draw_labels=True)
test_ax.set_title("Path doesn\'t show", y=1.5)

# Including 0 meridian in extent shows the path
test_ax1 = fig.add_axes((0, 0, 0.5, 0.5), projection=ccrs.Mercator(), label="map2")
test_ax1.plot(test_df['lon'], test_df['lat'], color="blue", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax1.set_extent([-130, 0, 30, 40], crs=ccrs.Geodetic())
test_ax1.coastlines()
test_ax1.gridlines(draw_labels=True)
test_ax1.set_title("Path shows (blue dot near San Francisco)", y=1.1)

plt.show()

编辑

(带有简化的可重现示例)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

test_df = pd.DataFrame()
test_df['lon'] = np.linspace(-120, -60, num=1000)
test_df['lat'] = 38

test_df1 = pd.DataFrame()
test_df1['lon'] = np.linspace(-120, -60, num=1001)
test_df1['lat'] = 38


fig = plt.figure()

meridian=0

test_ax = fig.add_axes((0, 0, 1, 0.6), projection=ccrs.Mercator())
test_ax.plot(test_df['lon'], test_df['lat'], color="blue", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax.coastlines()
test_ax.set_extent((-125, meridian, 36, 38))
gl = test_ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_left = False
test_ax.set_title('Path with {} points, eastern edge={}'.format(len(test_df),meridian))


test_ax1 = fig.add_axes((0, 0.7, 1, 0.6), projection=ccrs.Mercator())
test_ax1.plot(test_df1['lon'], test_df1['lat'], color="red", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax1.coastlines()
test_ax1.set_extent((-125, meridian, 36, 38))
gl1 = test_ax1.gridlines(draw_labels=True)
gl1.xlabels_top = False
gl1.ylabels_left = False
test_ax1.set_title('Path with {} points, eastern edge={}'.format(len(test_df1),meridian))


meridian=-10

test_ax2 = fig.add_axes((0, 1.4, 1, 0.6), projection=ccrs.Mercator())
test_ax2.plot(test_df['lon'], test_df['lat'], color="black", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax2.coastlines()
test_ax2.set_extent((-125, -10, 36, 38))
gl2 = test_ax2.gridlines(draw_labels=True)
gl2.xlabels_top = False
gl2.ylabels_left = False
test_ax2.set_title('Path with {} points, eastern edge={}'.format(len(test_df),meridian))


test_ax3 = fig.add_axes((0, 2.1, 1, 0.6), projection=ccrs.Mercator())
test_ax3.plot(test_df1['lon'], test_df1['lat'], color="green", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax3.coastlines()
test_ax3.set_extent((-125, -10, 36, 38))
gl3 = test_ax3.gridlines(draw_labels=True)
gl3.xlabels_top = False
gl3.ylabels_left = False
test_ax3.set_title('Path with {} points, eastern edge={}'.format(len(test_df1),meridian))

plt.show()

【问题讨论】:

  • 它可以与transform=ccrs.PlateCarree() 一起使用吗? (老实说,我不确定大地坐标变换是否对纬度/经度坐标有意义)
  • 我认为应该......认为我的知识肯定存在差距,因为我是新手。尽管形状(和位置)正确,但它在两种情况下都能正常工作。
  • 我刚刚尝试过,如果我的分数超过 1000 分,并且设置了test_ax1.set_extent([-130, 0, 30, 40], crs=ccrs.Geodetic()),它就可以工作。如果我把它拿出来,或者如果我将第二个值设置为小于 0 的任何值,它就不会显示数据。这可能是一个错误吗?
  • 感谢@MikeSperry 的试用和编辑。这似乎是一个错误。我添加了一个 Github 问题:github.com/SciTools/cartopy/issues/1357
  • @virtualmic 只是为了进一步确认这似乎是一个 cartopy 问题,而不是更一般的 pyplot 问题,当我将代码减少到 test_ax3 = fig.add_axes((0, 0.8, 1, 0.15))test_ax3.plot(test_df1['lon'], test_df1['lat'], color="green", linewidth=4, alpha=1.0) 时,数据绘制完全正常跨度>

标签: matplotlib cartopy


【解决方案1】:

鉴于 cartopy 似乎存在一些问题,我认为最好的解决方法是将您的数据分成小于 1000 个的块,然后绘制它的所有部分。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

test_df = pd.DataFrame()
test_df['lon'] = np.linspace(-120, -60, num=1001)
test_df['lat'] = 38

fig = plt.figure()

meridian=0

test_ax = fig.add_axes((0, 0.05, 1, 0.3), projection=ccrs.Mercator())
test_ax.plot(test_df['lon'], test_df['lat'], color="red", 
                       linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax.coastlines()
test_ax.set_extent((-125, meridian, 36, 38))
gl = test_ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_left = False
test_ax.set_title('Path with {} points, eastern edge={}'.format(len(test_df),meridian))

meridian=-10

test_ax3 = fig.add_axes((0, 0.55, 1, 0.3), projection=ccrs.Mercator())
# plot first 500
test_ax3.plot(test_df['lon'][:500], test_df['lat'][:500], color="green", 
                        linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
# plot to cover the gap
test_ax3.plot(test_df['lon'][499:501], test_df['lat'][499:501], color="blue", 
                        linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
# plot last 501
test_ax3.plot(test_df['lon'][500:], test_df['lat'][500:], color="yellow", 
                        linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax3.coastlines()
test_ax3.set_extent((-125, -10, 36, 38))
gl3 = test_ax3.gridlines(draw_labels=True)
gl3.xlabels_top = False
gl3.ylabels_left = False
test_ax3.set_title('Path with {} points, eastern edge={}'.format(len(test_df),meridian))

plt.show()

对于1001分的情况,我只是分成500分的一段和501分的一段。

由于您正在绘制线,因此我还添加了绘图以覆盖间隙,放大时显示为蓝色。

如果您也在绘制点,那么设置间隙填充而不是重叠两个部分的原因就出现了,如下所示:

test_ax3.plot(test_df['lon'][:500], test_df['lat'][:500], color="green",
                        linewidth=1, alpha=1.0, transform=ccrs.Geodetic(), marker='.')
# plot to cover the gap
test_ax3.plot(test_df['lon'][499:501], test_df['lat'][499:501], color="blue",
                        linewidth=1, alpha=1.0, transform=ccrs.Geodetic(), marker=None)
# plot last 501
test_ax3.plot(test_df['lon'][500:], test_df['lat'][500:], color="yellow",
                        linewidth=1, alpha=1.0, transform=ccrs.Geodetic(), marker='.')

通过分离填充物,您可以确保没有重复点,如果您的 alpha 值小于1.0,这可能是个问题。

将此应用于您的原始数据,您可以创建一个函数来循环遍历数据帧,块的大小等于您想要的任何大小。:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pandas as pd

PLOT_LIMIT = 1000

df_csv_url = 'http://www.sharecsv.com/dl/76dd767525a37180ca54cd1d9314b9dc/t1.csv'
test_df = pd.read_csv(df_csv_url)
figure_params = { 'width': 9.6, 'height': 5.4 }

fig = plt.figure(
        figsize=(figure_params["width"], figure_params["height"])
    )

print(len(test_df['lon']))

def ax_plot(test_ax, test_df):
    # this function will loop over the dataframe in chunks equal to PLOT_LIMIT
    len_df = len(test_df)
    n=0
    for i in range(len_df//PLOT_LIMIT):
        test_ax.plot(test_df['lon'][1000*i:1000*(i+1)], test_df['lat'][1000*i:1000*(i+1)], color="blue",
                        linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
        if (len_df-((n+1)*PLOT_LIMIT)) != 0:
            test_ax.plot(test_df['lon'][(1000*i)-1:(1000*(i+1))+1], test_df['lat'][(1000*i)-1:(1000*(i+1))+1], color="blue",
                            linewidth=4, alpha=1.0, transform=ccrs.Geodetic(), marker='None')
        n+=1

    test_ax.plot(test_df['lon'][1000*n:], test_df['lat'][1000*n:], color="blue",
                    linewidth=4, alpha=1.0, transform=ccrs.Geodetic())


test_ax1 = fig.add_axes((0, 0.55, 1, 0.45), projection=ccrs.Mercator(), label="map1")
ax_plot(test_ax1, test_df)
test_ax1.coastlines()
test_ax1.gridlines(draw_labels=True)
test_ax1.set_title("Path shows", y=1.5)

# Including 0 meridian in extent shows the path
test_ax2 = fig.add_axes((0, 0.1, 1, 0.45), projection=ccrs.Mercator(), label="map2")
ax_plot(test_ax2, test_df)
test_ax2.set_extent([-130, -30, 30, 40], crs=ccrs.Geodetic())
test_ax2.coastlines()
test_ax2.gridlines(draw_labels=True)
test_ax2.set_title("Path shows (blue dot near San Francisco)", y=1.1)

plt.show()

如您所见,您现在应该可以灵活地设置地图上的查看窗口。我还没有检查过像穿越反子午线这样的边缘情况,但在提出的情况下它是有效的。

【讨论】:

  • 非常感谢您提供详细的答案和解决方法。我会尝试一下,让你知道它是如何进行的,因为我的实际应用程序比这些场景要复杂一些。
  • 我已经能够找出另一种解决方法(我仍然需要在实际场景中尝试,但在当前情况下有效)。我已将其添加为另一个答案。
  • 我为我的实际用例尝试了另一种解决方法(我必须绘制超过 7000 个图)并且效果很好。现在使用它...感谢@Mike Sperry 的所有帮助!
【解决方案2】:

我已经找到了另一种解决方法。如果在使用plot 函数之前对点进行了转换(而不是传递transform 参数),则它可以工作。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

test_df = pd.DataFrame()
test_df['lon'] = np.linspace(-120, -60, num=1000)
test_df['lat'] = 38

test_df1 = pd.DataFrame()
test_df1['lon'] = np.linspace(-120, -60, num=1001)
test_df1['lat'] = 38


fig = plt.figure()

meridian=0

test_ax = fig.add_axes((0, 0, 1, 0.6), projection=ccrs.Mercator())
test_ax.plot(test_df['lon'], test_df['lat'], color="blue", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax.coastlines()
test_ax.set_extent((-125, meridian, 36, 38))
gl = test_ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_left = False
test_ax.set_title('Path with {} points, eastern edge={}'.format(len(test_df),meridian))


test_ax1 = fig.add_axes((0, 0.7, 1, 0.6), projection=ccrs.Mercator())
test_ax1.plot(test_df1['lon'], test_df1['lat'], color="red", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax1.coastlines()
test_ax1.set_extent((-125, meridian, 36, 38))
gl1 = test_ax1.gridlines(draw_labels=True)
gl1.xlabels_top = False
gl1.ylabels_left = False
test_ax1.set_title('Path with {} points, eastern edge={}'.format(len(test_df1),meridian))


meridian=-10

test_ax2 = fig.add_axes((0, 1.4, 1, 0.6), projection=ccrs.Mercator())
test_ax2.plot(test_df['lon'], test_df['lat'], color="black", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax2.coastlines()
test_ax2.set_extent((-125, -10, 36, 38))
gl2 = test_ax2.gridlines(draw_labels=True)
gl2.xlabels_top = False
gl2.ylabels_left = False
test_ax2.set_title('Path with {} points, eastern edge={}'.format(len(test_df),meridian))


test_ax3 = fig.add_axes((0, 2.1, 1, 0.6), projection=ccrs.Mercator())
test_ax3.plot(test_df1['lon'], test_df1['lat'], color="green", linewidth=4, alpha=1.0, transform=ccrs.Geodetic())
test_ax3.coastlines()
test_ax3.set_extent((-125, -10, 36, 38))
gl3 = test_ax3.gridlines(draw_labels=True)
gl3.xlabels_top = False
gl3.ylabels_left = False
test_ax3.set_title('Path with {} points, eastern edge={}'.format(len(test_df1),meridian))


test_ax4 = fig.add_axes((0, 2.8, 1, 0.6), projection=ccrs.Mercator())
# Instead of transforming within the plot function, transform and then plot
transformed_points = ccrs.Mercator().transform_points(ccrs.Geodetic(), test_df1['lon'].values, test_df1['lat'].values)
test_ax4.plot([p[0] for p in transformed_points], [p[1] for p in transformed_points], color="green", linewidth=4, alpha=1.0)
test_ax4.coastlines()
test_ax4.set_extent((-125, -10, 36, 38))
gl3 = test_ax4.gridlines(draw_labels=True)
gl3.xlabels_top = False
gl3.ylabels_left = False
test_ax4.set_title('Path with {} prior transformed points, eastern edge={}'.format(len(test_df1),meridian))


plt.show()

【讨论】:

    猜你喜欢
    • 2023-04-02
    • 1970-01-01
    • 2020-11-18
    • 1970-01-01
    • 2012-04-19
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2016-08-15
    相关资源
    最近更新 更多