Projection from Lambert Conformal Conic to Lat/Lon

Specific Domains Geo
,
1

Hello Community

After taking a look at packages like ArchGDAL I have still no idea how to project coordinates from an arbitrary Lambert Conformal Conic projection, defined by standard_parallel, latitude_of_projection_origin and longitude_of_central_meridian, to geographical latitude / longitude coordinates.
Unfortunately, there is no EPSG number available.

Is there a “function” like (lat,lon) = Lambert_Conformal_Conical_to_Lat_Lon(standard_parallel,latitude_of_projection_origin,longitude_of_central_meridian)?

If not, I would have to implement some equations like that Link to equations

Thanks for any replies, Dieter

2

Please check this post using GMT.
As your projection uses only one standard parallel, believe you have to set lat_1 equal to lat_2.

It should be possible to do this also using proj4.jl

3

Thanks rafael.guerra

I have also seen this post but is GMT not only usable for ploting maps instead of converting points from one projection into another? But I will take a closer look at this post.

4

All numerical transformations should be available as arrays in GMT.
Copying @joa-quim for expert advice.

NB:
Is the data from North America, Europe or some other region?

5

The data is from Europe, this is the NetCDF output from the EuroCordex Gepotential Climate Forecast:

Lambert_Conformal                    standard_parallel                    49.5                                                                      
Lambert_Conformal                    latitude_of_projection_origin        49.5                                                                      
Lambert_Conformal                    grid_mapping_name                    lambert_conformal_conic                                                   
Lambert_Conformal                    longitude_of_central_meridian        10.5
6

We may need to collect a bit more info, like provided here, to define the right projection.

7

Off course GMT does it. Both xy2lonlat and mapproject will do the job but the latter has a more complicated manual (it does a lot more than coordinate conversions) that was not yet converted do .md (big job).

Just find the proj4 string that defines Lambert data and follow the instructions of xy2lonlat

help?> xy2lonlat
search: xy2lonlat

  xy2lonlat(xy::Matrix{<:Real}, s_srs::String; t_srs=::String="+proj=longlat +datum=WGS84")

  or

  xy2lonlat(D::GMTdataset, s_srs::String; t_srs=::String="+proj=longlat +datum=WGS84")

  Computes the inverse projection from XY to LonLat in the given projection. The output is assumed to be in WGS84. If that isn't right, pass the
  appropriate projection info via the t_srs option.

  Parameters
  ––––––––––––

    •  xy: The input data. It can be a Matrix, or a GMTdataset (or vector of it)

    •  s_srs: The data projection system. This can be a PROJ4 or a WKT string

    •  t_srs: The target SRS. If the default is not satisfactory, provide a new projection info (PROJ4 or WKT)

  Returns
  –––––––––

  A Matrix if input is a Matrix or a GMTdadaset if input had that type
2 Likes
8

Here is an example using Proj4.jl directly. Like GMT and ArchGDAL, which use PROJ as well, you can create the proj or WKT string yourself, see also the documentation here: Lambert Conformal Conic — PROJ 9.1.0 documentation.

using Proj4

lcc = "+proj=lcc +lon_0=-90 +lat_1=33 +lat_2=45"
trans = Proj4.Transformation(lcc, "EPSG:4326", always_xy=true)

julia> trans((0.0, 0.0))
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 -90.0
   1.2722218725854067e-14
1 Like
9

xy2lonlat && lonlat2xy use PROJ too :slight_smile: but using EPSGs is not implemented in those functions yet.

EDIT: which is silly because there is a epsg2proj function (probably added later)

julia> epsg2proj(4326)
"+proj=longlat +datum=WGS84 +no_defs"
1 Like
10

Unfortunately, the NetCDF File does not contain more information as provided in

Lambert_Conformal                    standard_parallel                    49.5                                                                      
Lambert_Conformal                    latitude_of_projection_origin        49.5                                                                      
Lambert_Conformal                    grid_mapping_name                    lambert_conformal_conic                                                   
Lambert_Conformal                    longitude_of_central_meridian        10.5
11

Thanks joa-quim

xy2lonlat([-2825,-2825],"+proj=lcc +lat_1=49.5 +lat_2=49.5 +lat_0=49.5 +lon_0=10.5 +x_0=0 +y_0=0 +a=6371000 +units=km +no_defs",t_srs="+proj=longlat +datum=WGS84")

gives the correct result which is also identical to the self programmed set of equations:

function XY_Lambert_Conformal_Conical_2_Lat_Lon(x,y,λ0,ϕ0,ϕ1)
    
    R = 6371 # Radius of Earth

    λ0_rad = λ0*π/180
    ϕ0_rad = ϕ0*π/180
    ϕ1_rad = ϕ1*π/180

    x_skaliert = x/R
    y_skaliert = y/R

    n = sin(ϕ1_rad)

    F = (cos(ϕ1_rad)*(tan(π/4 + ϕ1_rad/2))^n)/n
    ρ0 = F*(cot(π/4 + ϕ0_rad/2))^n

    θ = atan(x_skaliert/(ρ0 - y_skaliert))
    ρ = sign(n)*sqrt(x_skaliert^2 + (ρ0 - y_skaliert)^2)

    ϕ_rad = 2atan((F/ρ)^(1/n)) - π/2
    λ_rad = λ0_rad + θ/n

    ϕ = ϕ_rad*180/π
    λ = λ_rad*180/π

    return λ,ϕ

end
2 Likes