从Python中的大数据集(IMS雪盖)中对二维纬度、经度和数据进行子集设置

Question

我心中有一个"简单"的目标，但很难实现。我有三个二维阵列lat、lon和data，它们的尺寸都是(2457624576)。前两个是变量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)具有相同的维度，然后我可以使用它来屏蔽原始数据，但不能对其进行子集化

在CCD_ 10中使用相同的条件会产生一个阵列形状的CCD_。我认为这里的问题是，在2-D阵列的每个点上都需要满足多个条件，并且不能保证这会导致连续的矩形子矩阵。然而，为了绘制这些数据，我真的需要对它们进行子集划分，或者找到另一种方法来详细说明它们。

我是不是错过了一些显而易见的东西？有没有其他聪明的方法可以在不出现任何错误的情况下绘制原始数据？

数据来源：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公里分辨率数据时遇到这篇文章，我有一个穷人解决方案，它会很好地工作。

内存错误是由plot例程引起的，但数组仍然可以存储到python中。因此，我没有使用where函数进行子集设置，而是尝试使用像这样的显式索引

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

然后用plot.imshow()绘制结果，而不需要绘制等高线或使用地图投影来了解数据的外观。这让我可以选择指数跨度，我感兴趣的区域就在里面。现在我可以在没有记忆错误的情况下绘制轮廓图了。

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