Python（IMS snowcover）中的大数据集的二维纬度、经度和数据子集

Question

我有一个 "simple" objective 的想法，结果证明很难实现。我有三个二维数组 lat、lon 和 data，它们的维度均为 (24576, 24576)。前两个是我试图在地图上绘制的变量 data 的纬度和经度坐标。我正在从二进制文件和文本文件的组合中读取所有这些数据，因此任何预处理操作实际上是不可能的，需要在 Python 脚本中完成。

考虑到数组的维数，由于内存限制，即使在地球上的一小块区域上选择底图投影，也几乎不可能直接绘制数据。我已经尝试过并在尝试使用 basemap.contourf.

绘图时出现内存错误

因此，我需要在将数组传递给轮廓函数之前对其进行子集化。我尝试了很多东西，但似乎没有任何效果。我的想法是做这样的事情

lat_bnds, lon_bnds = [35, 50], [5, 20]
condition=((lats > lat_bnds[0]) & (lats < lat_bnds[1])) & (lons > lon_bnds[0]) & (lons < lon_bnds[1])

其中lats和lons是二维坐标数组。这会生成一个二维布尔数组，它具有与 lats（或等效的 lons）相同的维度，然后我可以使用它来屏蔽原始数据但不对其进行子集化。

在 numpy.where 中使用相同的条件会生成一个我无法使用的形状为 (2, 2564856) 的数组。我认为这里的问题是二维数组的每个点都需要满足多个条件，并且不能保证这会导致连续的矩形子矩阵。然而，为了绘制这些数据，我真的需要对它们进行子集化，或者找到另一种方法来详细说明它们。

我是不是遗漏了什么明显的东西？有没有其他聪明的方法来绘制原始数据而不会遇到任何错误？

数据来源： http://nsidc.org/data/G02156

读取数据的简短脚本假定您下载了必要的文件：

import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit

with open('IMS1kmLats.24576x24576x1.double', 'rb') as f:
 data = np.fromfile(f, dtype='<d', count=24576*24576)
 lats = np.reshape(data, [24576, 24576], order='F')

with open('IMS1kmLons.24576x24576x1.double', 'rb') as f:
 data = np.fromfile(f, dtype='<d', count=24576*24576)
 lons = np.reshape(data, [24576, 24576], order='F')

widths=np.full((24576), 1, dtype=int).tolist()
data=np.array(pd.read_fwf('ims2017312_1km_v1.3.asc', skiprows=30, 
widths=widths, lineterminator='\n', header=None))

Answer 1

我不确定你是如何得到那个形状的数组的，但是考虑一下，如果你想要原始数组的一个子集并且你有你的条件，你可以从原始列表中切片。

a = np.random.randint(10,50, size=(50,2))  # ---- generate coordinates

array([[44, 11],
       [40, 36],
       [19, 26],
       ..., 
       [33, 26],
       [42, 12],
       [15, 25]])

lats = a[:, 1]  # ---- latitude is Y-axis
lons = a[:, 0]  # ---- longitude is X-axis
lat_bnds, lon_bnds = [35, 50], [5, 20]  # ---- your desired bounds

condition =((lats > lat_bnds[0]) & (lats < lat_bnds[1])) & 
            (lons > lon_bnds[0]) & (lons < lon_bnds[1])

condition
array([False, False, False, ..., False, False, False], dtype=bool)

condition.shape  # => (50,)

a[condition]     # ---- slice the original array

array([[11, 45],
       [15, 40],
       [15, 43],
       [15, 49]])

希望这能让你朝着你想去的方向前进。

Answer 2

万一有人在尝试使用 Python 绘制 IMS 1 公里分辨率数据时遇到这个 post，我有一个穷人解决方案，它可以正常工作。

绘图例程引发内存错误，但数组仍可存储到 python。因此，我没有使用 where 函数进行子集化，而是尝试使用像

这样的显式索引

lat_subset=lat[imin:imax, jmin:jmax]

然后使用 plot.imshow() 绘制结果而不进行等值线图或使用地图投影以了解数据的外观。这使我能够以我感兴趣的区域在内部的方式选择索引跨度。现在我可以在没有内存错误的情况下绘制等高线图了。

我有一个带有笔记本的小型存储库，其中显示了如何绘制这些数据：https://github.com/guidocioni/snow_ims/blob/mistral/plot_ims.ipynb。它的优点是可以直接在线读取文件，但考虑到坐标文件的大小，仍然需要下载它们。

Python（IMS snowcover）中的大数据集的二维纬度、经度和数据子集

Subsetting 2-D latitudes, longitudes and data from a big dataset in Python (IMS snowcover)

python

arrays

numpy

matplotlib-basemap