Using the PROJ.4 projection library in VB/VBA projects

by

Eric G. Miller
GIS Analyst
Wildlife & Habitat Data Analysis Branch
Department of Fish and Game

DFG GIS Mini-Conference
October 23-4, 2003
Norden, California

Outline

  1. Introduction
    1. Audience
    2. Purpose
    3. Background on PROJ.4
  2. Setup
    1. proj446_win32_bin.zip
      1. PATH=%PATH%;C:\proj\bin
      2. PROJ_LIB=C:\proj\nad
    2. nad2bin.exe
    3. proj_api.zip
  3. Demonstration
    1. Example.mdb
    2. Basic Usage
      1. pjInitPlus()
      2. pjTransform()
      3. pjFree()
    3. More Details and Pitfalls
      1. pjIsLatLong() & DEG_TO_RAD/RAD_TO_DEG
      2. pjErrNum() & pjErrorString()
      3. pjFree() > 1 → Boom!
      4. Projection limits checking, UTM zones
      5. Initialization strings: +proj +datum +nadgrids +nodefs (see docs for more)
  4. Additional Resources
    1. PROJ.4 project page (http://www.remotesensing.org/proj)
      1. Semi-out-of-date PDF docs.
      2. Man page for proj command
      3. Mailing list
    2. The proj_api repository on maphost (ftp://maphost.dfg.ca.gov/outgoing/Nafwb/proj4).
      1. proj446_win32_bin.zip
      2. proj_api.zip
      3. norden03/
      4. 00README.txt
      5. SETUP.TXT
      6. nad2bin.txt

Notes

  1. Introduction
    1. Audience

      The audience of this discussion is directed to people with at least a moderate understanding of Visual Basic programming. However, people without any programming experience may still benefit from this discussion.

    2. Purpose

      The purpose of this discussion is to demonstrate to Visual Basic programmers how to use the PROJ.4 projection library within their projects.

    3. Background on PROJ.4

      PROJ.4 is a cartographic projection library and a set of utility programs. The library was originally written by Gerald Evenden when he was with the USGS in the early 1990's. It is written in the ANSI C programming language. Although it was originally written for the UNIX platform, it was easily ported to other platforms. The current maintainer of the library is Frank Warmerdam, a programmer/consultant living in the Toronto area.

      The library currently supports over 100 projections and supports datum transformations (3-parameter, 7-parameter and grid shift methods). The API is relatively small and easy to use with three basic functions handling most of the work.

  2. Setup
    1. proj446_win32_bin.zip

      First, install the precompiled binaries. Using WinZip, PowerArchiver, or InfoZip's unzip, unpack the ZIP archive into the C:\ drive, retaining folders. The result should be a directory structure like:
      C:\proj\
      C:\proj\bin\
      C:\proj\bin\cs2cs.exe
      C:\proj\bin\geod.exe
      C:\proj\bin\nad2nad.exe
      C:\proj\bin\proj.dll
      C:\proj\bin\proj.exe
      C:\proj\nad\
      C:\proj\nad\alaska.lla
      C:\proj\nad\conus.lla



      C:\proj\nad\world
      C:\proj\README.TXT

      1. PATH=%PATH%;C:\proj\bin

        Now, modify the PATH environment variable. Right-click on "My Computer" and select "Properties" from the context menu. In the "System Properties" window, select the "Advanced" tab, then click on the "Environment Variables" button. You may either modify the PATH variable for the system as a whole or just your user by adding ;C:\proj\bin to the end. Select the "Edit..." button to edit an existing variable. If you'd like to just modify your user environment and you don't have a PATH variable, create it with the "New..." button.

        setting the PATH variable
      2. PROJ_LIB=C:\proj\nad

        In a similar fashion to setting the PATH variable, create the new environment variable PROJ_LIB and set its value to C:\proj\nad. Click "Ok" on the "Environment Variables" window and then "Ok" on the "System Properties" window to finish basic PROJ.4 installation.

    2. nad2bin.exe

      In order to perform datum transformations using the more accurate NADCON method, we must convert the "conus.lla" ASCII file to a binary representation with the nad2bin.exe program. Launch an MS-DOS (Command.com) window and change directory to C:\proj\nad. Then run the program as follows:
      C:\proj\nad>nad2bin conus < conus.lla

    3. proj_api.zip

      Now unpack the proj_api.zip file. A good location would be C:\proj\proj_api\ but anywhere will do. Then copy the file proj_api.dll to C:\proj\bin\. Now we are ready to use the wrapper library in our VB/VBA projects.

  3. Demonstration
    1. Example.mdb

      PROJ.4 example

      In the example Access database, several coordinate systems have been predefined. The user need only select the input and output coordinate systems, enter the coordinates to convert, and click the Convert button.

      When the output coordinate system is geographic, the user may select one of three formats (Degrees/Minutes/Seconds, Degrees/Minutes, Degrees). On input, all three formats are handled if the coordinate system is geographic. Currently, the only recognized delimiter is the space character.

    2. Basic Usage

      In order to program with the library, there are three basic functions that do most of the work. Since the library wrapper is not a COM component, we must declare each function we wish to use.

      1. pjInitPlus()

        This function initializes a coordinate system definition taking a string of parameters as its only argument and returning a pointer or "handle" to an opaque C structure. The declaration for this function looks like:

        Declare Function pjInitPlus Lib "proj_api.dll" (ByVal definition As String) As Long

        The syntax basically tells the VB interpreter/compiler that there is a function named "pjInitPlus" in the library "proj_api.dll" that takes a string argument by value and returns a long integer. VB understands that its idea of a string passed by value to an external library will convert it to a C string. The return value is really a C pointer, but on the Wintel architechture, we can safely cast a pointer to a long integer and back since they are the same size (32 bits).

        To use the function we need to create a string with the parameters defining our coordinate system. The following is an example:

        Dim strCoordSys As String
        Dim lngCoordSys As Long

        strCoordSys = "+proj=utm +zone=11 +datum=NAD83"
        lngCoordSys = pjInitPlus(strCoordSys)

      2. pjTransform()

        The pjTransform function encapsulates all the steps necessary to transform from one coordinate system to another. It is declared as follows:

        Declare Function pjTransform Lib "proj_api.dll" _
        (ByVal srcPointer As Long, ByVal dstPointer As Long, ByVal count As Long, _
        ByVal offset As Long, xs As Double, ys As Double, zs As Double) As Long

        This declaration is somewhat more complex. The first two arguments are our "pointers" to coordinate system definitions as returned by pjInitPlus. The first is the input coordinate system and the second is the output coordinate system. The count parameter indicates how many coordinates will be transformed. The offset parameter indicates how many array elements the function should skip to fetch the next X, Y, or Z value. The last three arguments are the X, Y and Z values and actually use a trick in order to pass arrays to the C function. VB arguments are ByRef by default which is similar to passing a pointer to the object. By calling the function with a ByRef argument, and either using simple Double arguments or the first element of an array of Double, a pointer to the first element is passed rather than the actual value. The return value of the function is a status code, and is zero if the function was successful. Perhaps an example will clarify the matter.

        Dim lngInputSys as Long
        Dim lngOutputSys as Long
        Dim dblPoints(3) as Double
        Dim lngStatus as Long

        lngInputSys = pjInitPlus("+proj=utm +zone=11 +datum=NAD27")
        lngOutputSys = pjInitPlus("+proj=latlong +datum=NAD83")
        dblPoints(0) = 728878.58  ' X(0)
        dblPoints(1) = 3659676.74 ' Y(0)
        dblPoints(2) = 726892.55  ' X(1)
        dblPoints(3) = 3678124.22 ' Y(1)

        lngStatus = pjTransform(lngInputSys, lngOutputSys, 2, 2, dblPoints(0), dblPoints(1), vbNull)

        First, we create two coordinate system definitions. Then we populate a four element double array with two coordinates, alternating the X and Y UTM values. When we invoke pjTransform we pass 2 for the count and 2 for the offset to indicate successive X or Y values are two elements apart in the array. Notice how we pass the first element for the xs argument and the second element for the ys argument. Because of the ByRef nature of the argument, the C function is actually able to access the whole array. Since we don't have any height values, we pass vbNull for the zs argument. The underlying C function knows to check for this possibility with the height argument.

        Upon succesful completion (a zero return value), our array now holds the transformed coordinates. However, see below for handling geographic (lat/long) coordinates...

      3. pjFree()

        Because the library wrapper is not a COM component, there is no way for VB to tell the library when to release resources that are no longer needed. In our case, we have these "pointers" to coordinate system definitions which are taking up a little bit of memory. Should we lose the "pointer" at the end of a function, that memory will continue to be used, but we can no longer reference it (a "memory leak"). Therefore, we must call pjFree for each coordinate system definition we have created. This is similar to setting a COM object to Nothing explicity within VB code. The subroutine is declared as:

        Declare Sub pjFree Lib "proj_api.dll" (ByVal pjPointer As Long)

        The subroutine is simple to use: pjFree(lngCoordSys)

    3. More Details and Pitfalls
      1. pjIsLatLong() & DEG_TO_RAD/RAD_TO_DEG

        Above, an example of transforming UTM coordinates to geographic coordinates was provided. However, we left out a crucial detail. Geographic coordinates must be passed in and are returned in radians rather than degrees. To handle these cases we can use the function pjIsLatLong and the constants DEG_TO_RAD and RAD_TO_DEG.

        Before transforming coordinates, you can check if the input coordinates are geographic by calling pjIsLatLong with the "pointer" to the input system definition. If so, just multiply each coordinate value with the DEG_TO_RAD constant. Similarly, we can check the output system definition and transform coordinates with RAD_TO_DEG if needed. For example:

        If pjIsLatLong(lngInputCoordSys) <> 0 Then
           dblPoint(0) = dblPoint(0) * DEG_TO_RAD
           dblPoint(1) = dblPoint(1) * DEG_TO_RAD
        End If
        lngStatus = pjTransform(lngInputCoordSys, lngOutputCoordSys, 1, 2, dblPoint(0), dblPoint(1), vbNull)
        If pjIsLatLong(lngOutputCoordSys) <> 0 Then
           dblPoint(0) = dblPoint(0) * RAD_TO_DEG
           dblPoint(1) = dblPoint(1) * RAD_TO_DEG
        End If

      2. pjErrNum() & pjErrorString()

        Because we are accessing a plain C library, errors can not be raised and caught as is typically done in VB code. Instead we must check the return values of functions for error indications and possibly fetch the error code with another function call.

        For instance, if pjInitPlus returns zero, we know the call wasn't successful, but we don't know why. We can attempt to handle the error condition by then fetching a copy of a global error number with pjErrNum. This code number can be used to look up standard errors from PROJ.4. The VB function pjErrorString in the example database attempts to return a meaningful error string. As an example:

        lngPointer = pjInitPlus("+proj=utm +zone=10")
        If lngPointer = 0 Then
           lngErr = pjErrNum()
           goto HandleErrors
        End If

        HandleErrors:
           If lngErr <> 0 Then
             MsgBox "Projection Error #" & lngErr & " " & pjErrorString(lngErr)
           End If

        In the above example, an error would be created because the coordinate system lacks a datum or the minimum ellipsoid parameters. Do not just fetch the error number without having called a function that indicates an error occurred. Doing so will likely result in a false positive with a meaningless error number.

      3. pjFree() > 1 → Boom!

        A word of caution when using pjFree is in order. Calling the function with a bogus value, such as a previously freed "pointer" will almost certainly crash your application. To protect yourself, you might want to set free'd pointers to zero after calling pjFree. A zero value translates to a NULL pointer in C, which is always valid in a call to the underlying C free function.

      4. Projection limits checking, UTM zones

        Unfortunately, this version of the library does not specify an easy means of detecting if input coordinates are well outside the reasonable limits of the output coordinate system. Depending on the scenario, you may want to proactively code boundary checks for input data.

        Additionally, there is no facility for the library reporting the appropriate output UTM zone given an "arbitrary" input. The only real solution at this point is to project to geographic coordinates with the appropriate output datum, and then determine which zone is the appropriate output zone. UTM zones are six degrees wide and California's zones 10 and 11 meet at 120 degrees west...

      5. Initialization strings: +proj +datum +nadgrids +nodefs +init (see docs for more)

        So far, this discussion has not really gone into detail on the arguments to initialize a projection definition. In the example database, the table "projections_lut" has both the displayed values and the initilization strings used in the form interface. This table holds several of the more common projections our organization might deal with.

        +proj
        Most initilization strings will have a +proj argument specifying the projection. A list of abbreviations and their full names can be output with the proj command line utility using the argument -l (hyphen ell).
        +datum
        The easiest and typically best way to specify most of the underlying parameters of the geographic reference frame is by specifying the datum. Typical values we would use are NAD27 and NAD83. You may also use WGS84 though NAD83 is nearly the same.
        +nadgrids
        This argument specifies a list of grid files to use in datum transformations. It can be left off the argument list and fairly sensible defaults will be used. However, you may want to specify that only the "conus" file should be consulted. Never use this argument with datums/ellipsoids that are essentially equivalent to WGS84/NAD83. Doing so will give you erroneous results.
        +nodefs
        This argument simply specifies that when the library creates a coordinate system definition it should not use any default values for unspecified/incomplete parameters. Generally, it's probably a good idea to specify this so you don't get surprised.
        +init
        It is possible to predefine coordinate systems and reference them by filename and id number. Such definitions should be placed in the PROJ_LIB directory. The State Plane coordinate system definitions are referenced this way in the example database

        There are a number of other arguments that haven't been covered here. Please consult the documentation for specifics on the recognized arguments for the different projections.

  4. Additional Resources
    1. PROJ.4 project page (http://www.remotesensing.org/proj)
      1. Semi-out-of-date PDF docs.
      2. Man page for proj command
      3. Mailing list
    2. The proj_api repository on maphost (ftp://maphost.dfg.ca.gov/outgoing/Nafwb/proj4).
      1. proj446_win32_bin.zip - Precompiled binary PROJ.4 distribution
      2. proj_api.zip - Wrapper library and example
      3. norden03/ - Norden presentation
      4. 00README.txt - A README and non-warranty
      5. SETUP.TXT - Some notes on setup
      6. nad2bin.txt - Notes on using nad2bin.exe