Julia - 在兰伯特等角投影中转换纬度、经度

Posted 2023-03-27

【中文标题】Julia - 在兰伯特等角投影中转换纬度、经度

【英文标题】：Julia - converting Latitude, Longitude in Lambert conformal projection

【发布时间】：2022-01-23 11:47:57

【问题描述】：

我已从 HRRR 下载数据，类似于此笔记本中的 grib2 文件：

https://nbviewer.org/github/microsoft/AIforEarthDataSets/blob/main/data/noaa-hrrr.ipynb

我现在希望将数据用于特定的经度、纬度。但是我不知道如何将我的（经度，纬度）转换为矩阵数据中的网格坐标。

笔记本提到“HRRR 数据来自 Lambert 保形投影”。 （见单元格 8）。

我查看了 GMT 包，他们似乎可以处理 Lambert 保形投影：https://docs.juliahub.com/GMT/EoU0j/0.30.1/proj_examples/。

但是如何转换坐标呢？

下面的代码好像可以转换，但是我觉得这不是给Lambert的，看了GMT documentation之后，无法调整命令中的设置。

lat=37.0; lon=-119.0;

gmt("mapproject -J+proj=merc", [lat;lon])
VectorGMTdataset with 1 segments
First segment DATA
Global BoundingBox:    [-1.3247019404399555e7, 4.118821159351122e6]
First seg BoundingBox: [-1.3247019404399555e7, 4.118821159351122e6]
2×1 MatrixFloat64:
  4.118821159351122e6
 -1.3247019404399555e7

【参考方案1】：

发现经纬度其实在grib文件里，所以不用转换：

  using GRIB
  f = GribFile(grb2_filename)
  
  lons, lats, values = data(Message(f))
  # lons in range [225.90452026573686, 299.0828072281622]   = [-134.09547973426314, -60.91719277183779]
  # lats in range [21.138123000000018, 52.61565330680793]

所以我们可以只查找最近的经纬度索引，读取values中对应的值。

由于纬度和经度都接近110公里，我将距离最小化如下：

(min_error, coord) = findmin(abs.(lats .- lat) .+ abs.(lons .- 360 .- lon))
(0.020456228700048484, CartesianIndex(269, 548))

values[coord]
294.8936767578125

虽然这实际上并没有回答标题问题，但它满足了我当前的需求，也许对其他人有用。