Code Reorganization and Geodesic Routines

REQUIRE

@@ -1 +1 @@
-julia 0.3
+julia 0.4

gen/generate_projection_codes.jl

@@ -15,7 +15,7 @@ for filename in FILENAMES
     lines = readlines(proj_file)
     close(proj_file)
 
-    println("$filename = Compat.@compat Dict(")
+    println("$filename = Dict(")
     for i=1:length(lines)
         m = parse_projection(lines[i])
         if is(m, nothing)

src/Proj4.jl

@@ -1,171 +1,31 @@
 module Proj4
 
-using Compat
-
-export transform, transform!, Projection, is_latlong, is_geocent
-
-## Types
-
-# Projection context.  TODO: Will this be exposed?
-#type Context
-#    rep::Ptr{Void} # Pointer to internal projCtx struct
-#end
-
-"""
-Cartographic projection type
-"""
-type Projection
-    #ctx::Context   # Projection context object
-    rep::Ptr{Void} # Pointer to internal projPJ struct
-end
-
-## Low-level interface pieces
 const libproj = "libproj"
 
-# The following functions are generally named after the associated C API
-# functions, but without the pj prefix.
-
-"Free C datastructure associated with a projection.  For internal use only!"
-function _free(proj::Projection)
-    @assert proj.rep != C_NULL
-    ccall((:pj_free, libproj), Void, (Ptr{Void},), proj.rep)
-    proj.rep = C_NULL
-end
-
-"Get human readable error string from proj.4 error code"
-function _strerrno(code::Cint)
-    bytestring(ccall((:pj_strerrno, libproj), Cstring, (Cint,), code))
-end
-
-"Get global errno string in human readable form"
-function _strerrno()
-    _strerrno(unsafe_load(ccall((:pj_get_errno_ref, libproj), Ptr{Cint}, ())))
-end
-
-"Get projection definition string in the proj.4 plus format"
-function _get_def(proj::Projection)
-    @assert proj.rep != C_NULL
-    opts = 0 # Apparently obsolete argument, not used in current proj source
-    bytestring(ccall((:pj_get_def, libproj), Cstring, (Ptr{Void}, Cint), proj.rep, opts))
-end
-
-"""Low level interface to libproj transform, allowing user to specify strides"""
-function _transform!(src::Projection, dest::Projection, point_count, point_stride, x, y, z)
-    @assert src.rep != C_NULL && dest.rep != C_NULL
-    ccall((:pj_transform, libproj), Cint,
-          (Ptr{Void}, Ptr{Void}, Clong, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
-          src.rep, dest.rep, point_count, point_stride, x, y, z)
-end
-
-## High level interface
-
-"""
-Construct a projection from a string in proj.4 "plus format"
-
-The projection string `proj_string` is defined in the proj.4 format,
-with each part of the projection specification prefixed with '+' character.
-For example:
+export Projection, # proj_types.jl
+       transform, transform!,  # proj_functions.jl
+       is_latlong, is_geocent, compare_datums, spheroid_params,
+       xy2lonlat, xy2lonlat!, lonlat2xy, lonlat2xy!
 
-    `wgs84 = Projection("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")`
-"""
-function Projection(proj_string::ASCIIString)
-    proj = Projection(ccall((:pj_init_plus, libproj), Ptr{Void}, (Cstring,), proj_string))
-    if proj.rep == C_NULL
-        # TODO: use context?
-        error("Could not parse projection string: \"$proj_string\": $(_strerrno())")
-    end
-    finalizer(proj, _free)
-    proj
-end
-
-# Pretty printing
-Base.print(io::IO, proj::Projection) = print(io, strip(_get_def(proj)))
-Base.show(io::IO, proj::Projection) = print(io, "Projection(\"$proj\")")
-
-"""
-Return true if the projection is a geographic coordinate system (lon,lat)
-"""
-function is_latlong(proj::Projection)
-    @assert proj.rep != C_NULL
-    ccall((:pj_is_latlong, libproj), Cint, (Ptr{Void},), proj.rep) != 0
-end
-
-"""
-Return true if the projection is a geocentric coordinate system
-"""
-function is_geocent(proj::Projection)
-    @assert proj.rep != C_NULL
-    ccall((:pj_is_geocent, libproj), Cint, (Ptr{Void},), proj.rep) != 0
-end
-
-
-"""
-Transform between geographic or projected coordinate systems in place.
+include("projection_codes.jl") # ESRI and EPSG projection strings
+include("proj_capi.jl") # low-level C-facing functions (corresponding to src/proj_api.h)
 
-See transform() for details
-
-"""
-function transform!(src::Projection, dest::Projection, position::Array{Float64,2}; radians::Bool=false)
-    npoints = size(position,1)
-    ncomps = size(position,2)
-    if ncomps != 2 && ncomps != 3
-        error("position must be Nx2 or Nx3")
-    end
-    if !radians && is_latlong(src)
-        position[:,1:2] = deg2rad(position[:,1:2])
-    end
-    P = pointer(position)
-    x = P
-    y = P + sizeof(Cdouble)*npoints
-    z = (ncomps < 3) ? C_NULL : P + 2*sizeof(Cdouble)*npoints
-    err = _transform!(src, dest, npoints, 1, x, y, z)
-    if err != 0
-        error("transform error: $(_strerrno(err))")
-    end
-    if !radians && is_latlong(dest)
-        position[:,1:2] = rad2deg(position[:,1:2])
-    end
-    position
+__m = match(r"(\d+).(\d+).(\d+),.+", _get_release())
+version_release = parse(Int, __m[1])
+version_major  = parse(Int, __m[2])
+version_minor = parse(Int, __m[3])
+if version_release >= 4 && version_major >= 9 && version_minor >= 1
+    export geod_direct, geod_inverse, destination, ellps_distance
+    include("proj_geodesic.jl") # low-level C-facing functions (corresponding to src/geodesic.h)
 end
 
-transform!(src::Projection, dest::Projection, position::Vector{Float64}; radians::Bool=false) =
-    transform!(src, dest, reshape(position,(1,length(position))), radians=radians)
-
+include("proj_types.jl") # type definitions for proj objects
+include("proj_functions.jl") # user-facing proj functions
 
-"""
-Transform between geographic or projected coordinate systems
-
-Args:
-
-    src      - Source coordinate system definition
-    dest     - Destination coordinate system definition
-    position - An Nx2 or Nx3 array of coordinates to be transformed in place.
-               For geographic coordinate systems, the first two columns are
-               the *longitude* and *latitude*, in that order.
-    radians  - If true, treat geographic lon,lat coordinates as radians on
-               input and output.
-
-Returns:
-
-    position - Transformed position
-"""
-transform(src::Projection, dest::Projection, position::Array{Float64,2}; radians::Bool=false) =
-    transform!(src, dest, copy(position), radians=radians)
-
-transform(src::Projection, dest::Projection, position::Vector{Float64}; radians::Bool=false) =
-    transform!(src, dest, copy(reshape(position,(1,length(position)))), radians=radians)
-
-transform{T <: Real}(src::Projection, dest::Projection, position::Array{T,2}; radians::Bool=false) =
-    transform!(src, dest, @compat(map(Float64,position)), radians=radians)
-
-transform{T <: Real}(src::Projection, dest::Projection, position::Vector{T}; radians::Bool=false) =
-    transform!(src, dest, reshape(@compat(map(Float64,position)),(1,length(position))), radians=radians)
-
-"Get a string describing the underlying version of libproj in use"
-function libproj_version()
-    bytestring(ccall((:pj_get_release, libproj), Cstring, ()))
-end
+@doc "Get a global error string in human readable form" ->
+error_message() = _strerrno()
 
-include("projection_codes.jl")
+@doc "Get a string describing the underlying version of libproj in use" ->
+version() = _get_release()
 
 end # module

src/proj_capi.jl

@@ -0,0 +1,163 @@
+# The following functions are generally named after the associated C API
+# functions, but without the pj prefix.
+
+immutable ProjUV
+    u::Cdouble
+    v::Cdouble
+end
+
+@doc "forward projection from Lat/Lon to X/Y (only supports 2 dimensions)" ->
+function _fwd!(lonlat::Vector{Cdouble}, proj_ptr::Ptr{Void})
+    xy = ccall((:pj_fwd, libproj), ProjUV, (ProjUV, Ptr{Void}), ProjUV(lonlat[1], lonlat[2]), proj_ptr)
+    _errno() == 0 || error("forward projection error: $(_strerrno())")
+    lonlat[1] = xy.u; lonlat[2] = xy.v
+    lonlat
+end
+
+@doc "Row-wise projection from Lat/Lon to X/Y (only supports 2 dimensions)" ->
+function _fwd!(lonlat::Array{Cdouble,2}, proj_ptr::Ptr{Void})
+    for i=1:size(lonlat,1)
+        xy = ccall((:pj_fwd, libproj), ProjUV, (ProjUV, Ptr{Void}),
+                   ProjUV(lonlat[i,1], lonlat[i,2]), proj_ptr)
+        lonlat[i,1] = xy.u; lonlat[i,2] = xy.v
+    end
+    _errno() == 0 || error("forward projection error: $(_strerrno())")
+    lonlat
+end
+
+@doc "inverse projection from X/Y to Lat/Lon (only supports 2 dimensions)" ->
+function _inv!(xy::Vector{Cdouble}, proj_ptr::Ptr{Void})
+    lonlat = ccall((:pj_inv, libproj), ProjUV, (ProjUV, Ptr{Void}),
+                   ProjUV(xy[1], xy[2]), proj_ptr)
+    _errno() == 0 || error("inverse projection error: $(_strerrno())")
+    xy[1] = lonlat.u; xy[2] = lonlat.v
+    xy
+end
+
+@doc "Row-wise projection from X/Y to Lat/Lon (only supports 2 dimensions)" ->
+function _inv!(xy::Array{Cdouble,2}, proj_ptr::Ptr{Void})
+    for i=1:size(xy,1)
+        lonlat = ccall((:pj_inv, libproj), ProjUV, (ProjUV, Ptr{Void}),
+                       ProjUV(xy[i,1], xy[i,2]), proj_ptr)
+        xy[i,1] = lonlat.u; xy[i,2] = lonlat.v
+    end
+    _errno() == 0 || error("inverse projection error: $(_strerrno())")
+    xy
+end
+
+function _init_plus(proj_string::ASCIIString)
+    proj_ptr = ccall((:pj_init_plus, libproj), Ptr{Void}, (Cstring,), proj_string)
+    if proj_ptr == C_NULL
+        # TODO: use context?
+        error("Could not parse projection: \"$proj_string\": $(_strerrno())")
+    end
+    proj_ptr
+end
+
+@doc "Free C datastructure associated with a projection. For internal use!" ->
+function _free(proj_ptr::Ptr{Void})
+    @assert proj_ptr != C_NULL
+    ccall((:pj_free, libproj), Void, (Ptr{Void},), proj_ptr)
+end
+
+@doc "Get human readable error string from proj.4 error code" ->
+function _strerrno(code::Cint)
+    bytestring(ccall((:pj_strerrno, libproj), Cstring, (Cint,), code))
+end
+
+@doc "Get global errno string in human readable form" ->
+function _strerrno()
+    _strerrno(_errno())
+end
+
+@doc "Get error number"
+function _errno()
+    unsafe_load(ccall((:pj_get_errno_ref, libproj), Ptr{Cint}, ()))
+end
+
+@doc "Get projection definition string in the proj.4 plus format" ->
+function _get_def(proj_ptr::Ptr{Void})
+    @assert proj_ptr != C_NULL
+    opts = 0 # Apparently obsolete argument, not used in current proj source
+    bytestring(ccall((:pj_get_def, libproj), Cstring, (Ptr{Void}, Cint), proj_ptr, opts))
+end
+
+@doc "Low level interface to libproj transform, C_NULL can be passed in for z, if it's 2-dimensional" ->
+function _transform!(src_ptr::Ptr{Void}, dest_ptr::Ptr{Void}, point_count::Int, point_stride::Int,
+                     x::Ptr{Cdouble}, y::Ptr{Cdouble}, z::Ptr{Cdouble})
+    @assert src_ptr != C_NULL && dest_ptr != C_NULL
+    err = ccall((:pj_transform, libproj), Cint, (Ptr{Void}, Ptr{Void}, Clong, Cint, Ptr{Cdouble}, Ptr{Cdouble},
+                Ptr{Cdouble}), src_ptr, dest_ptr, Clong(point_count), Cint(point_stride), x, y, z)
+    err != 0 && error("transform error: $(_strerrno(err))")
+end
+
+function _transform!(src_ptr::Ptr{Void}, dest_ptr::Ptr{Void}, position::Vector{Cdouble})
+    @assert src_ptr != C_NULL && dest_ptr != C_NULL
+    ndim = length(position)
+    @assert ndim >= 2
+
+    x = pointer(position)
+    y = x + sizeof(Cdouble)
+    z = (ndim == 2) ? Ptr{Cdouble}(C_NULL) : 2*sizeof(Cdouble)
+
+    _transform!(src_ptr, dest_ptr, 1, 1, x, y, z)
+    position
+end
+function _transform!(src_ptr::Ptr{Void}, dest_ptr::Ptr{Void}, position::Array{Cdouble,2})
+    @assert src_ptr != C_NULL && dest_ptr != C_NULL
+    npoints, ndim = size(position)
+    @assert ndim >= 2
+
+    x = pointer(position)
+    y = x + sizeof(Cdouble)*npoints
+    z = (ndim == 2) ? Ptr{Cdouble}(C_NULL) : 2*sizeof(Cdouble)*npoints
+    #z = Ptr{Cdouble}((ndim == 2) * (x + 2*sizeof(Cdouble)*npoints))
+
+    _transform!(src_ptr, dest_ptr, npoints, 1, x, y, z)
+    position
+end
+
+function _is_latlong(proj_ptr::Ptr{Void})
+    @assert proj_ptr != C_NULL
+    ccall((:pj_is_latlong, libproj), Cint, (Ptr{Void},), proj_ptr) != 0
+end
+
+function _is_geocent(proj_ptr::Ptr{Void})
+    @assert proj_ptr != C_NULL
+    ccall((:pj_is_geocent, libproj), Cint, (Ptr{Void},), proj_ptr) != 0
+end
+
+@doc "Get a string describing the underlying version of libproj in use" ->
+function _get_release()
+    bytestring(ccall((:pj_get_release, libproj), Cstring, ()))
+end
+
+@doc """
+Fetch the internal definition of the spheroid as a tuple (a, es), where
+
+    a = major_axis
+    es = eccentricity squared
+
+""" ->
+function _get_spheroid_defn(proj_ptr::Ptr{Void})
+    major_axis = Ref{Cdouble}()
+    eccentricity_squared = Ref{Cdouble}()
+    ccall((:pj_get_spheroid_defn, libproj), Void, (Ptr{Void}, Ptr{Cdouble}, Ptr{Cdouble}),
+          proj_ptr, major_axis, eccentricity_squared)
+    major_axis[], eccentricity_squared[]
+end
+
+@doc "Returns true if the two datums are identical, otherwise false." ->
+function _compare_datums(p1_ptr::Ptr{Void}, p2_ptr::Ptr{Void})
+    Bool(ccall((:pj_compare_datums, libproj), Cint, (Ptr{Void}, Ptr{Void}), p1_ptr, p2_ptr))
+end
+
+# Unused/untested
+
+# @doc """
+# Return the lat/long coordinate system on which a projection is based.
+# If the coordinate system passed in is latlong, a clone of the same will be returned.
+# """ ->
+# function _latlong_from_proj(proj_ptr::Ptr{Void})
+#     ccall((:pj_latlong_from_proj, libproj), Ptr{Void}, (Ptr{Void},), proj_ptr)
+# end