JuliaGeo / Proj.jl

Julia wrapper for the PROJ cartographic projections library
MIT License
48 stars 9 forks source link

Optional +x_0, +y_0 projection parameters #62

Closed graziano-giuliani closed 2 years ago

graziano-giuliani commented 2 years ago

I am having a problem with optional projection parameters. I am using this proj4 strings:

pj4string = "+proj=ob_tran +o_proj=longlat +ellps=sphere +no_defs +to_meter=0.017453292519943295 +o_lon_p=198.0 +o_lat_p=39.25 +x_0=2.997602990300261 +y_0=1.7528341677779053 +R=6371229" reference = "+proj=latlong +R=6371229"

Using it on the proj command line utility the above proj4 string gives me the below result:

proj +proj=ob_tran +o_proj=longlat +ellps=sphere +no_defs +to_meter=0.017453292519943295 +o_lon_p=198.0 +o_lat_p=39.25 +x_0=2.997602990300261 +y_0=1.7528341677779053 +R=6371229 0 0 9.75 49.68

I am getting different answer using the Julia implementation:

using Proj dirtrans = Proj4.Transformation(pj4string, reference) dirtrans((0.0,0.0)) 2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2): 152.81767357346513 -47.4333686809606

I can get the above from command line by using the inverse projection and removing the optional parameters:

invproj +proj=ob_tran +o_proj=longlat +ellps=sphere +no_defs +to_meter=0.017453292519943295 +o_lon_p=198.0 +o_lat_p=39.25 +R=6371229 -f "%f" 152.817674 -47.433369

What am I missing?

visr commented 2 years ago

I'm not sure! I checked if it may be due to different versions, but using the bundled proj from PROJ_jll I get the same result:

import PROJ_jll

# from test/applications.jl
function read_cmd(cmd)
    bytes = read(cmd)
    str = String(bytes)
    # on Windows there are CFLR line endings
    lf_only = replace(str, "\r" => "")
    no_tabs = replace(lf_only, "\t" => " ")
    return String(strip(no_tabs))
end

PROJ_jll.proj() do proj
    read_cmd(pipeline(IOBuffer("0 0"),
        `$proj +proj=ob_tran +o_proj=longlat +ellps=sphere +no_defs +to_meter=0.017453292519943295 +o_lon_p=198.0 +o_lat_p=39.25 +x_0=2.997602990300261 +y_0=1.7528341677779053 +R=6371229`))
end
# -> "9.75 49.68"
visr commented 2 years ago

I assume the difference probably lies in that Transformation uses proj_create_crs_to_crs and proj_trans, which corresponds to the cs2cs application, which does a different thing than the proj application (https://proj.org/apps/index.html).

graziano-giuliani commented 2 years ago

Thanks, this has solved the problem for me. I was using the wrong application.