雷达极坐标数据处理为栅格数据

原始文件是一个个扇形组织的，
扇形1：起始角度(AngleStart)、角度步长(AngleDelta)、随距离(r)变化的值；
扇形2：......（同上）
....
....
扇形N

如何将上述表达方式转换为M*N的栅格图像呢？有没有方便高效的算法，能够判断给定距离下，某一个角度范围（AngleStart到AngleStart+AngleDelta）包含了哪些像素呢？

——AngleDelta基本上为1deg左右，如果忽略AngleDelta，直接按照x=r*sin(AngleStart)，y=r*cos(AngleStart)来计算的话，会有很多像素没有值。

这个不难。给个样例数据，我帮你写个demo

收集够了信息，可以写一个demo了。题主有360个扇区，每个扇区100个库（题主提供的资料误写为1000个库），总计36000个库，每个库对应一个数值，没错吧？那我先生成3个数组：一个是100行360列的tg数组（theta_grid），一个是100行360列的rg数组（r_grid），一个是100行360列的dg数组（data_grid）。

>>> import numpy as np
>>> t = np.arange(360)
>>> r = np.arange(100)
>>> tg, rg = np.meshgrid(t, r)
>>> tg.shape, rg.shape
((100, 360), (100, 360))
>>> dg = np.empty(tg.shape)
>>> dg.fill(np.nan)
>>> dg.shape
(100, 360)

现在，dg数组填充值是np.nan（not a number）。从excele读出数值，写入dg数组对应位置上，题主应该没有问题吧？接下来，我就用随机数来代替真正的数据，以便继续后面的代码。

>>> dg = np.random.random(tg.shape)

至此，数据读取工作完成，tg/rg/dg三个数组构成了一个高效完整的数据集。假如题主想检索6km~9km距离内、0°~5°范围内的库，可以这样写：

>>> w = np.where((rg>=40)&(rg<60)&(tg>=0)&(tg<5))
>>> d = dg[w] # 符合条件的库的数值
>>> t = tg[w] # 符合条件的库的theta
>>> r = rg[w] # 符合条件的r
>>> result = np.stack((t,r,d), axis=1) # 这是三个结果组合在一起的写法
>>> result.shape # 总共找到100个库
(100, 3)
>>> result[:5] # 只显示前5个库的theta/r/data
array([[0.00000000e+00, 4.00000000e+01, 5.66644807e-01],
       [1.00000000e+00, 4.00000000e+01, 8.72776055e-01],
       [2.00000000e+00, 4.00000000e+01, 3.82545134e-02],
       [3.00000000e+00, 4.00000000e+01, 8.24244873e-01],
       [4.00000000e+00, 4.00000000e+01, 8.73632778e-01]])

你在2021-06-28 20:53的消息中说，“数据包括360行（表示360个径向），1000列（表示每个径向上有1000个库）”，也就是总共有36万个数据。请问这些数据和你最新给出的网格数据之间是什么关系？如果这36万个数据不是你说的网格数据，哪这些网格数据从何而来？如果这36万个数据就是你说的网格数据，那我理解的没有问题——excel表给出的数据虽然是网格的，但只是形式上的网格，或者说逻辑上的网格，并非你想象的平面上的网格。我给你提供的算法，返回的就是对应范围内所有库的扇区索引、距离索引和数值，没有遗漏。

这是最终的答案。约定长度单位使用千米（km），角度单位使用度（°），目标网格精度是0.15km。关于0°定义和角度的正负是逻辑上的，如何使用则取决于题主导入数据时的预处理，比如，excel文件中第i扇区是正北方向的，就需要将该扇区角度设置为0°。

考虑到1°夹角的扇面在150km处开口宽度约5.2km，相当于17.5个网格宽度，故将1个扇区分为20等分，足可覆盖最大库区内的所有网格。

>>> import numpy as np
>>> k, cell = 20, 0.15 # 扇区分割数，目标网格精度
>>> t = np.arange(0, 360, 1/k)
>>> r = np.arange(1000) * 0.15 + 0.075
>>> tg, rg = np.meshgrid(t, r)
>>> dg = np.random.random((1000,360))

上面代码的最后一行，是模拟的excel文件读入，请题主用实际数据替换此行。此时dg仍然是1000行360的数组，首行最近末行最远，首列为0°（即正北方向），末列为359°。

>>> xs = rg * np.cos(np.radians(tg))
>>> ys = rg * np.sin(np.radians(tg))
>>> xs = xs.ravel()
>>> ys = ys.ravel()
>>> ds = np.repeat(dg.ravel(), k)
>>> xs.shape, ys.shape, ds.shape
((7200000,), (7200000,), (7200000,))

这段代码完成了极坐标系到平面直角坐标系的转换。然后，再生成目标网格，大小是2001x2001。

>>> cells = complex(0, int(1+(300/0.15)))
>>> cells
2001j
>>> lats, lons = np.mgrid[150:-150:cells, -150:150:cells]
>>> data = np.empty(lats.shape)
>>> data.fill(np.nan)
>>> lats.shape, lons.shape, data.shape
((2001, 2001), (2001, 2001), (2001, 2001))

接下来就比较魔幻了，请细细揣摩。

>>> xi = np.int32(xs/cell) + 1000
>>> yj = np.int32(ys/cell) + 1000
>>> data[(yj, xi)] = ds

data就是题主需要的数据了。检索数据的方法，和这个类似：在给定的角度范围内，按照1/k°分割，转为平面直角坐标系上的坐标之后，再除以网格精度转为数据的行列号，然后就可以提取数据了。

老师这个应该是读出了对应范围内库的值嘛。那怎么把这个库的值赋给创建栅格网呢