Introduction

There are many Geographic Information System tools, but the most basic is a simple conversion between geographic and projection coordinates. There are several ways to do this, but most that we'll discuss rely upon the proj library and tools. You can access the proj library through the following methods

1. the command line interface,
2. the pyproj python package,
3. and the GRASS program.

The most basic is the command line interface, so start there. The rest of the tools rely on knowledge there in.

proj package

The most basic access to the proj library is through the command line interface. The proj package provides command line tools to convert between geographic coordinates (e.g., longitude/latitude) and map projections (e.g. Lambert Conformal Conic) coordinates. The core of the proj package or the proj and invproj commands.

proj/invproj cli

The proj (invproj) command line utility converts from a geographic (projected) coordinates to projected (geographic) coordinates. The conversion between any two map projection systems can be performed by first converting to an intermediate geographic projection. Further documentation of the proj library is provided at http://trac.osgeo.org/proj/.

Tutorial

In air quality work, the three most common projections are Geographic, Lambert Conformal Conic and Universal Transverse Mercator.

This tutorial will start will provide a basic usage of the proj and invproj commands. Then, it will conclude with an example converting from the coordinates in a nested domain to the coordinates of an outer domain.

proj/invproj args and options

proj and invproj have the same arguments and options. First, we'll discuss basic usage. The projection is defined by a series of arguments in the form +key=value, where key indicates the type of value being supplied. The first key, value pair are "proj" and the projection type: lcc or utm. lonlat/latlon are not relevant for proj and invproj which convert to and from Geographic to map projections.

Convert Geographic (lon/lat) to Map Projection

For a UTM and LCC projection, the following tutorials convert from the longitude latitude to the map projection. The coordinates will be for Gainesville, FL (-82d19'29.97" 29d39'7.19"). Each example will make use of the -f option that is detailed in the "Relevant Options" section.

UTM

This example will convert the longitude/latitude of Gainesville, FL to the kilometer definition for UTM in zone 17.

$ echo "-82d19'29.97\" 29d39'7.19\"" | proj +proj=utm +zone=17 +to_meter=1000 -f "%.6f"

Output: 371.758225 3280.958525

LCC

This example will convert the longitude/latitude of Gainesville, FL to the kilometer definition in the Lambert Conformal Conic projection used by the AQMEII Continental United States domain.

$ echo "-82d19'29.97\" 29d39'7.19\"" | proj +ellps=sphere +proj=lcc +lat_1=33. +lat_2=45. +lat_0=40 +lon_0=-97. +x_0=2556000. +y_0=1728000 -f "%.6f"

Output: 3978657.607076 694180.092449

The data can also be output in grid units using the +to_meter option

$ echo "-82d19'29.97\" 29d39'7.19\"" | proj +ellps=sphere +proj=lcc +lat_1=33. +lat_2=45. +lat_0=40 +lon_0=-97. +x_0=2556000. +y_0=1728000 +to_meter=12000. -f "%.6f"

Output: 331.554801 57.848341

Converting Between Map Projections

LCC to UTM

The commands below convert between the arbitrary projections used in the example.

$ echo 1422657.607076 3722031.663847 | cs2cs +ellps=sphere +proj=lcc +lat_1=33. +lat_2=45. +lat_0=40 +lon_0=-97. +x_0=2556000. +y_0=1728000 +to +proj=utm +zone=17 -f "%.6f"

Output: 371758.224684 3280958.525460 0.000000

To input the data in grid units, the radius of the sphere must be adjusted.

$ echo 118.554801 310.169305 | cs2cs +ellps=sphere +proj=lcc +lat_1=33. +lat_2=45. +lon_0=-97. +to_meter=12000 +to +proj=utm +zone=17 -f "%.6f"

Output: 371758.228849 3280958.520873 0.000000

To convert the output unit, the +to_meter argument must also be specified in the +to projection.

LCC Grid to LCC Grid

These examples are meant to be practical and are derived from the Houston-Galveston-Brazoria 2005/2006 SIP modeling. The projection parameters are summarized here to prevent loss of data.

Lambert Conformal Conic Parameters

Lambert Grid Definitions

Grid Conversions

Finding the lower left corner is a special where the units and/or +to_meter property are irrelevant for the input project. In this tutorial, it is specified in both input and +to projection for completeness.

$ echo 0 0 | cs2cs +proj=lcc +a=6370000.0m +lon_0=100w +lat_1=30n +lat_2=60n +lat_0=40n +x_0=-394000m +y_0=1154000m +to_meter=2000 +to +proj=lcc +a=6370000.0m +lon_0=100w +lat_1=30n +lat_2=60n +lat_0=40n +x_0=-356000m +y_0=1228000m +to_meter=4000

Output: 9.50 18.50 0.00

Note that the corner has a 1 2-km grid-cell overlap with the 4-km domain. This is standard and expected. In a nested domain, this grids values will always be zero.

Relevant options

 -m mult

mult will be multiplied to the output of proj or the input to invproj will be divided by mult. This option is not available in cs2cs and so is not used in the tutorial, but has the same effect as +to_meter when mult in the form 1:X where X is the grid +to_meter value.

Effectively, this lets you convert from meter units to grid units fractional.

-f format

Format is a printf format string to control the form of the output values. For inverse projections, the output will be in degrees when this option is employed. The default format is "%.2f" for forward projection and DMS for inverse.

pyproj Python

Easy Example

# Import Proj from pyproj from pyproj import Proj # Define projection for 2-km and 4-km Texas domains p2k = Proj('+proj=lcc +a=6370000.0m +lon_0=100w +lat_1=30n +lat_2=60n +lat_0=40n +x_0=-394000m +y_0=1154000m +to_meter=2000', preserve_units = True) p4k = Proj('+proj=lcc +a=6370000.0m +lon_0=100w +lat_1=30n +lat_2=60n +lat_0=40n +x_0=-356000m +y_0=1228000m +to_meter=4000', preserve_units = True) # Convert from 2-km grid lower left corder to 4-km space p4k(*p2k(0, 0, inverse = True))

Output: (9.500000000000174, 18.499999999999535)