Example 1.6: Defining Projections and Extents

Whoa! What happened here?

You've just experienced on-the-fly reprojection. When a PROJECTION object is defined in the mapfile, MapServer passes the layer and projection info to the PROJ library (formerly 'PROJ.4') which does the reprojection. For more information about the PROJ library, please visit https://proj.org/ .

This example attempts to shed some light on the projection support in MapServer.

This is what the mapfile looks like (Example1-6.map):

The first thing you might have noticed with our mapfile is the original EXTENT has been commented out and the new EXTENT values don't look anything like latitude and longitude values. Also, we have added a PROJECTION object near the top and within each of our layers.

Let's have a look at the new object and parameter:

EXTENT 201621.496941 -294488.285333 1425518.020722 498254.511514
The extent we provide to MapServer needs to be in the same units as the output projection. Since Lambert Azimuthal Equal-Area's units are in meters, we have to give the new extent in meters. We can calculate the new extent using QGIS or some other GIS package, or you can use PROJ's cs2cs utility. The command below is what can be used to reproject the original extent values:

cs2cs +proj=latlong +datum=WGS84 +to +proj=laea +ellps=clrk66 +lat_0=45 +lon_0=-100

After typing the command, enter the southwestern coordinate pair (the lower left coordinates), separated by a space:  -97.5 41.619778
to which the "cs2cs" utility gives the values:

            208398.01       -372335.44 0.000
          

I then type the northeastern coordinate pair (upper right corner coordinates), again separated by space:  -82.122902 49.38562 and the following values are returned:

            1285308.08      632638.93 0.000
          

You will notice that "cs2cs" returns a set of three values. You can ignore the third value, 0.000, as it means to represent altitute (which we're not using). Anyway, we can not present the EXTENT as:

            EXTENT 208398.01 -372335.44 1285308.08 632638.93
          

But it doesn't match up the extent above, you say. Well, that's what usually happens when you reproject--the map doesn't necessarily get centered as you'd expect. You can fudge around it, taking a few thousand meters off to the left and adding a few to the right. Or, you can use a graphical GIS package to give you the extent. Here are the instructions for getting EXTENT using QGIS.

MapServer can take projection definitions in two ways. The first is shown in the mapfile's output PROJECTION object (the first PROJECTION block). This is the traditional way of passing parameters to the PROJ library. The other way is to use the EPSG codes. These codes are standard projection codes (or spatial reference identifiers) as defined by the European Petroleum Survey Group (EPSG). In the case of our Lambert Azimuthal Equal-Area projection, it has a code of "2163". If you comment out the four lines after PROJECTION and uncomment the line "init=epsg:2163", this will provide the same information to PROJ. If you want to learn more about EPSG codes, have a look at "/usr/local/share/proj/epsg" or "C:/proj/nad/epsg" ("/ms4w/proj/nad" in MS4W).

You can search for EPSG codes and see MapServer usage examples at http://epsg.io/ Also, check out the EPSG web site at http://www.epsg.org .

If you need to learn more about projections, please look through some of the links below:

http://www.geography.hunter.cuny.edu/mp/
http://en.wikipedia.org/wiki/Map_projection
https://store.usgs.gov/assets/mod/storefiles/PDF/16573.pdf

You can find more by searching the web for "map projection".

Back to Example 1.6 | Back to the Section 1 | Back to the Sections Page | Proceed to Example 1.7