CODEWITHSUNDEEP
julianoaa
Timothée HENRY
Timothée HENRY

Reputation: 14614

Julia - converting Latitude, Longitude in Lambert conformal projection

I have dowloaded data from the HRRR, similar to the grib2 file from this notebook:

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

I now wish to use the data for specific Longitude, Latitude. But I do not know how to convert my (Longitude, Latitude) to the grid coordinates in the matrix data.

The notebook mentions that “the HRRR data comes in the Lambert conformal projection.” (See cell 8).

I have looked at the GMT package, and they seem to handle the Lambert conformal projection: https://docs.juliahub.com/GMT/EoU0j/0.30.1/proj_examples/.

But how can I convert the coordinates?

The following code seems to convert, but I don't think this is for Lambert, and after looking at the GMT documentation, I am not able to adjust the settings in the command.

lat=37.0; lon=-119.0;

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

Upvotes: 2

Views: 246

Answers (1)

Timothée HENRY
Timothée HENRY

Reputation: 14614

I found out that the longitude and latitude were actually in the grib file, so there is no need to convert:

  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]

So we can just look for the indexes of the closest longitude and latitude, and read the corresponding value in values.

Since lat and lon degree both approximate to 110 kilometers, I will just minimize the distance as follows:

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

values[coord]
294.8936767578125

While this actually does not answer the title question, it answers my current need and perhaps will be useful for someone else.

Upvotes: 1

Related Questions