| Author: | Frank Warmerdam |
|---|---|
| Contact: | warmerdam at pobox.com |
| Revision: | $Revision: 8480 $ |
| Date: | $Date: 2009-01-29 12:04:25 -0800 (Thu, 29 Jan 2009) $ |
| Last Updated: | 2007/12/09 |
Table of Contents
MapServer supports rendering a variety of raster file formats in maps. The following describes some of the supported formats, and what capabilities are supported with what formats.
This document assumes that you are already familiar with setting up MapServer Mapfile, but does explain the raster specific aspects of mapfiles.
A simple raster layer declaration looks like this. The DATA file is interpreted relative to the SHAPEPATH, much like shapefiles.
LAYER
NAME "JacksonvilleNC_CIB"
DATA "Jacksonville.tif"
TYPE RASTER
STATUS ON
END
Though not shown rasters can have PROJECTION, METADATA, PROCESSING, MINSCALE, and MAXSCALE information. It cannot have labels, CONNECTION, CONNECTIONTYPE, or FEATURE information.
Rasters can be classified in a manner similar to vectors, with a few exceptions.
There is no need to specify a CLASSITEM. The raw pixel value itself (“[pixel]”) and, for paletted images, the red, green and blue color associated with that pixel value (“[red]”, “[green]” and “[blue]”) are available for use in classifications. When used in an evaluated expression the pixel, red, green and blue keywords must be in lower case.
LAYER
NAME "JacksonvilleNC_CIB"
DATA "Jacksonville.tif"
TYPE RASTER
STATUS ON
CLASSITEM "[pixel]"
# class using simple string comparison, equivelent to ([pixel] = 0)
CLASS
EXPRESSION "0"
STYLE
COLOR 0 0 0
END
END
# class using an EXPRESSION using only [pixel].
CLASS
EXPRESSION ([pixel] >= 64 AND [pixel] < 128)
STYLE
COLOR 255 0 0
END
END
# class using the red/green/blue values from the palette
CLASS
NAME "near white"
EXPRESSION ([red] > 200 AND [green] > 200 AND [blue] > 200)
STYLE
COLOR 0 255 0
END
END
# Class using a regular expression to capture only pixel values ending in 1
CLASS
EXPRESSION /*1/
STYLE
COLOR 0 0 255
END
END
END
As usual, CLASS definitions are evaluated in order from first to last, and the first to match is used. If a CLASS has a NAME attribute it may appear in a LEGEND. Only the COLOR, EXPRESSION and NAME parameters within a CLASS definition are utilized for raster classifications. The other styling or control information is ignored.
Raster classifications always take place on only one raster band. It defaults to the first band in the referenced file, but this can be altered with the BANDS PROCESSING directive. In particular this means that including even a single CLASS declaration in a raster layer will result in the raster layer being rendered using the one band classification rules instead of other rules that might have applied (such as 3 band RGB rendering).
As of MapServer 4.4 support has been added for classifying non-8bit raster inputs. That is input rasters with values outside the range 0-255. Mostly this works transparently but there are a few caveats and options to provide explicit control.
Classifying raster data in MapServer is accomplished by pre-classifying all expected input values and using that table of classification results to lookup each pixel as it is rendered. This is done because evaluating a pixel value against a series of CLASS definitions is relatively expensive to do for the hundreds of thousands of pixels in a typical rendered image.
For simple 8bit inputs, only 256 input values need to be pre-classified. But for non-8bit inputs more values need to be classified. For 16bit integer inputs all 65536 possible input values are pre-classified. For floating point and other input data types, up to 65536 values are pre-classified based on the maximum expected range of input values.
The PROCESSING directive can be used to override the range of values to be pre-classified, or the number of values (aka Buckets) in that range to classify. The SCALE=min,max PROCESSING directive controls the range. The SCALE_BUCKETS PROCESSING directive controls the number of buckets. In some cases rendering can be accelerated considerable by selecting a restricted range of input values and a reduced number of scaling values (buckets).
The following example classifies a floating raster, but only 4 values over the range -10 to 10 are classified. In particular, the values classified would be -7.5, -2.5, 2.5, and 7.5 (the middle of each “quarter” of the range). So those four value are classified, and one of the classification results is selected based on which value is closest to the pixel value being classified.
LAYER
NAME grid1
TYPE raster
STATUS default
DATA data/float.tif
PROCESSING "SCALE=-10,10"
PROCESSING "SCALE_BUCKETS=4"
CLASS
NAME "red"
EXPRESSION ([pixel] < -3)
STYLE
COLOR 255 0 0
END
END
CLASS
NAME "green"
EXPRESSION ([pixel] >= -3 and [pixel] < 3)
STYLE
COLOR 0 255 0
END
END
CLASS
NAME "blue"
EXPRESSION ([pixel] >= 3)
STYLE
COLOR 0 0 255
END
END
END
What raster formats are supported by MapServer is largely controlled by configuration time options. Some formats are considered to be built-in while the remainder are handled by the optional GDAL raster library.
More information on GDAL can be found at http://www.gdal.org, including the supported formats list. Some of the advanced MapServer raster features, such as resampling, RGB color cube generation and automatic projection capture only work with raster formats used through GDAL. GDAL is normally built and installed separately from MapServer, and then enabled during the build of MapServer using the –with-gdal configuration switch.
To find out what is built into a particular MapServer executable, use the -v flags to discover what build options are enabled. To find out what GDAL formats are available, the “gdalinfo –formats” command may be used. For example:
warmerda@gdal2200[124]% mapserv -v
MapServer version 4.4.0-beta2 OUTPUT=GIF OUTPUT=PNG OUTPUT=JPEG OUTPUT=WBMP
SUPPORTS=PROJ SUPPORTS=FREETYPE SUPPORTS=WMS_SERVER SUPPORTS=WMS_CLIENT
SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT SUPPORTS=WCS_SERVER SUPPORTS=FASTCGI
INPUT=EPPL7 INPUT=POSTGIS INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE DEBUG=MSDEBUG
warmerda@gdal2200[18]% gdalinfo --formats
Supported Formats:
GRASS (ro): GRASS Database Rasters (5.7+)
GTiff (rw+): GeoTIFF
NITF (rw+): National Imagery Transmission Format
HFA (rw+): Erdas Imagine Images (.img)
SAR_CEOS (ro): CEOS SAR Image
...
The following formats are potential builtins:
If built with INPUT=TIFF MapServer will have builtin support for reading TIFF or GeoTIFF files. The builtin TIFF support has some limitations with regard to the organization of files that can be read (no tiled, 16bit, RGB, or odd color models). This driver supports world files, or simple builtin GeoTIFF coordinates for georeferencing.
Full featured TIFF/GeoTIFF support is available through GDAL. Note that only GDAL supports tiled TIFF files and TIFF files with overviews. Tiled TIFF files with overviews pre-built are one of the highest performance ways of serving large raster images.
If GD is configured with GIF (OUTPUT=GIF) support, then MapServer will also be able to read GIF files for raster layers. The only way to georeference GIF files is with a world file.
If GD is not configured with GIF support, it may still be available in GDAL.
If GD is configured with PNG (OUTPUT=PNG) support, then MapServer will also be able to read PNG files for raster layers. The only way to georeference PNG files is with a world file.
If GD is not configured with PNG support, it may still be available in GDAL.
If MapServer is built with JPEG (INPUT=JPEG) support then greyscale JPEG files may be rendered in raster layers. RGB files (the more common kind) will not be able to be displayed. Georeferencing is via world files.
If MapServer is not built with native JPEG support, GDAL may still support the format. In this case RGB files are also supported (via the RGB color cube mechanism). Georeferencing is still via world file.
If MapServer is built with GDAL it is generally better to access all possible formats through GDAL rather than via the built-in drivers. The built-in drivers are less featureful, and not as well maintained.
When handling very large raster layers it is often convenient, and higher performance to split the raster image into a number of smaller images. Each file is a tile of the larger raster mosaic available for display. The list of files forming a layer can be stored in a shapefile with polygons representing the footprint of each file, and the name of the files. This is called a TILEINDEX and works similarly to the same feature in vector layers. The result can be represented in the MAP file as one layer, but MapServer will first scan the tile index, and ensure that only raster files overlapping the current display request will be opened.
The following example shows a simple example. No DATA statement is required because MapServer will fetch the filename of the raster files from the Location attribute column in the hp2.dbf file for records associated with polygons in hp2.shp that intersect the current display region. The polygons in hp2.shp should be rectangles representing the footprint of the corresponding file. Note that the files do not have to be all the same size, the formats can vary and they can even overlap (later files will be drawn over earlier ones); however, they must all be in the same coordinate system (projection) as the layer.
LAYER
NAME "hpool"
STATUS ON
TILEINDEX "hp2.shp"
TILEITEM "Location"
TYPE RASTER
END
The filenames in the tileindex are searched for relative to the SHAPEPATH or map file, not relative to the tileindex. Great care should be taken when establishing the paths put into the tileindex to ensure they will evaluate properly in use. Often it is easiest to place the tileindex in the SHAPEPATH directory, and to create the tileindex with a path relative to the SHAPEPATH directory. When all else fails, absolute paths can be used in tileindex, but then they cannot be so easily moved from system to system.
While there are many ways to produce TILEINDEX shapefiles for use with this command, one option is the gdaltindex program, part of the GDAL utility suite. The gdaltindex program will automatically generate a tile index shapefile from a list of GDAL supported raster files passed on the command line.
Usage: gdaltindex [-tileindex field_name] index_file [gdal_file]*
% gdaltindex doq_index.shp doq/*.tif
The gdaltindex program is built as part of GDAL. Prebuilt binaries for GDAL including the gdaltindex program can be downloaded as part of the OSGeo4W, FWTools and MS4W distributions.
See also
MapServer is able to resample GDAL rasters on the fly into new projections. Non-GDAL rasters may only be up or down sampled without any rotation or warping.
Raster warping kicks in if the projection appears to be different for a raster layer than for the map being generated. Warped raster layers are significantly more expensive to render than normal raster layers with rendering time being perhaps 2-4 times long than a normal layer. The projection and datum shifting transformation is computed only at selected points, and generally linearly interpolated along the scanlines (as long as the error appears to be less than 0.333 pixels.
In addition to reprojecting rasters, the raster warping ability can also apply rotation to GDAL rasters with rotational coefficients in their georeferencing information. Currently rotational coefficients won’t trigger raster warping unless the map and layer have valid (though matching is fine) projection definitions.
Traditionally MapServer has been used to produce 8 bit pseudo-colored map displays generated from 8bit greyscale or pseudocolored raster data. However, if the raster file to be rendered is actually 24bit (a red, green and blue band) then additional considerations come into play. Currently rendering of 24bit imagery is only supported via the GDAL renderer. The built-in PNG, JPEG and other drivers do not support 24bit input images.
If the output is still 8bit pseudo-colored (the IMAGEMODE is PC256 in the associated OUTPUT format declaration) then the full 24bit RGB colors for input pixels will be converted to a color in the colormap of the output image. By default a color cube is used. That is a fixed set of 175 colors providing 5 levels of red, 7 levels of green and 5 levels of blue is used, plus an additional 32 greyscale color entries. Colors in the input raster are mapped to the closest color in this color cube on the fly. This substantial degrades color quality, especially for smoothly changing images. It also fills up the colors table, limited to 256 colors, quite quickly.
A variation on this approach is to dither the image during rendering. Dithering selects a color for a pixel in a manner that “diffuses error” over pixels. In an area all one color in the source image, a variety of output pixel colors would be selected such that the average of the pixels would more closely approximate the desired color. Dithering also takes advantage of all currently allocated colors, not just those in the color cube. Dithering requires GDAL 1.1.9 or later, and is enabled by providing the PROCESSING “DITHER=YES” option in the mapfile. Dithering is more CPU intensive than using a simple color cube, and should be avoided if possible in performance sensitive situations.
The other new possibility for handling 24bit input imagery in MapServer 4.0 or later, is to produce 24bit output images. The default “IMAGETYPE png24” or “IMAGETYPE jpeg” declaration may be used to produce a 24bit PNG output file, instead of the more common 8bit pseudo-colored PNG file. The OUTPUTFORMAT declaration provides for detailed control of the output format. The IMAGEMODE RGB and IMAGEMODE RGBA options produce 24bit and 32bit (24bit plus 8bit alpha/transparency) for supported formats.
As of MapServer 4.0, the PROCESSING parameter was added to the LAYER of the Mapfile. It is primarily used to pass specialized raster processing options to the GDAL based raster renderer. The following processing options are supported in MapServer 4.0 and newer.
This turns on error diffusion mode, used to convert 24bit images to 8bit with error diffusion to get better color results.
Example:
PROCESSING "DITHER=YES"
This directive allows a specific band or bands to be selected from a raster file. If one band is selected, it is treated as greyscale. If 3 are selected, they are treated as red, green and blue. If 4 are selected they are treated as red, green, blue and alpha (opacity).
Example:
PROCESSING "BANDS=4,2,1"
This directive instructs the GDAL reader to pre-scale the incoming raster data. It is primarily used to scale 16bit or floating point data to the range 0-255, but can also be used to constrast stretch 8bit data. If an explicit min/max are provided then the input data is stretch (or squished) such that the minimum value maps to zero, and the maximum to 255. If AUTO is used instead, a min/max is automatically computed. To control the scaling of individual input bands, use the SCALE_1, SCALE_2 and SCALE_3 keywords (for red, green and blue) instead of SCALE which applies to all bands.
Example:
PROCESSING "SCALE=AUTO"
or
PROCESSING "SCALE_1=409,1203"
PROCESSING "SCALE_2=203,296"
PROCESSING "SCALE_3=339,1004"
This directive (MapServer 4.9+) instructs the GDAL reader to apply a custom LUT (lookup table) to one or all color bands as a form of on the fly color correction. If LUT is used, the LUT is applied to all color bands. If LUT_n is used it is applied to one color band (n is 1 for red, 2 for green, 3 for blue, 4 for alpha).
The LUT can be specified inline in the form:
<lut_spec> = <in_value>:<out_value>[,<in_value>:<out_value>]*
This essentially establish particular input values which are mapped to particular output values. The list implicitly begins with 0:0, and 255:255. An actual 256 entry lookup table is created from this specification, linearly interpolating between the values. The in values must be in increasing order. The LUT specification may also be in a text file with the <lut_spec> being the filename. The file contents should be in the same syntax, and the file is searched relative to the mapfile.
Example:
PROCESSING "LUT_1=red.lut"
PROCESSING "LUT_2=green.lut"
PROCESSING "LUT_3=blue.lut"
or
PROCESSING "LUT=100:30,160:128,210:200"
As a special case there is also support for GIMP format curve files. That is the text files written out by the Tools->Color->Curves tool. If this is specified as the filename then it will be internally converted into linear segments based on the curve control points. Note that this will not produce exactly the same results as the GIMP because linear interpolation is used between control points instead of curves as used in the GIMP. For a reasonable number of control points the results should be similar. Also note that GIMP color curve files include an overall “value” curve, and curves for red, green, blue and alpha. The value curve and the appropriate color curve will be composed internally to produce the final LUT.
Example:
PROCESSING "LUT=munich.crv"
Alter the precision with which colors need to match an entry in the color table to use it when producing 8bit colormapped output (IMAGEMODE PC256). Normally colors from a raster colormap (or greyscale values) need to match exactly. This relaxes the requirement to being within the specified color distance. So a COLOR_MATCH_THRESHOLD of 3 would mean that an existing color entry within 3 (sum of difference in red, green and blue) would be used instead of creating a new colormap entry. Especially with greyscale raster layers, which would normally use all 256 color entries if available, this can be a good way to avoid “stealing” your whole colormap for a raster layer. Normally values in the range 2-6 will give good results.
Example:
PROCESSING "COLOR_MATCH_THRESHOLD=3"
This option can be used to control the resampling kernel used sampling raster images. The default (and fastest) is NEAREST. AVERAGE will perform compute the average pixel value of all pixels in the region of the disk file being mapped to the output pixel (or possibly just a sampling of them). BILINEAR will compute a linear interpolation of the four pixels around the target location. This topic is discussed in more detail in MS RFC 4: MapServer Raster Resampling.
Resampling options other than NEAREST result in use of the generalized warper and can dramatically slow down raster processing. Generally AVERAGE can be desirable for reducing noise in dramatically downsampled data, and can give something approximating antialiasing for black and white linework. BILINEAR can be helpful when oversampling data to give a smooth appearance.
Example (chose one):
PROCESSING "RESAMPLE=NEAREST"
PROCESSING "RESAMPLE=AVERAGE"
PROCESSING "RESAMPLE=BILINEAR"
A new feature added in MapServer 4.4 is the ability to perform queries on rasters in a manner similar to queries against vector layers. Raster queries on raster layers return one point feature for each pixel matching the query. The point features will have attributes indicating the value of different bands at that pixel, the final rendering color and the class name. The resulting feature can be directly access in MapScript, or processed through templates much like normal vector query results. Only raster layers with a query TEMPLATE associated can be queried, even for the query methods that don’t actually use the query template (much like vector data).
Raster query supports QueryByPoint, QueryByRect, and QueryByShape. QueryByPoint supports single and multiple result queries. Other query operations such as QueryByIndex, QueryByIndexAdd, QueryByAttributes and QueryByFeature are not supported for raster layers.
Raster layers do not support saving queries to disk, nor query maps.
Raster queries return point features with some or all of the following attributes:
| x: | georeferenced X location of pixel. |
|---|---|
| y: | georeferenced Y location of pixel. |
| value_list: | a comma separated list of the values of all selected bands at the target pixel. |
| value_n: | the value for the n’th band in the selected list at this pixel (zero based). There is one value_n entry for each selected band. |
| class: | Name of the class this pixel is a member of (classified layers only). |
| red: | red component of the display color for this pixel. |
| green: | green component of the display color for this pixel. |
| blue: | blue component of the display color for this pixel. |
The red, green and blue attribute are intended to be the final color the pixel would be rendered with, but in some subtle cases it can be wrong (ie. classified floating point results). The selected bands are normally the band that would be used to render the layer. For a pure query-only layer BANDS PROCESSING directive can be used to select more bands than could normally be used in a render operation. For instance for a 7 band landsat scene a PROCESSING “BANDS=1,2,3,4,5,6,7” directive could be used to get query results for all seven bands in results to a query operation.
Care should be taken to avoid providing a large query area (selecting alot of pixels) as each selected pixel requires over 100 bytes of memory for temporary caching. The RASTER_QUERY_MAX_RESULT PROCESSING item can be used to restrict the maximum number of query results that will be returned. The default is one million which would take on the order of 100MB of RAM.
Query results can be returned as HTML via the normal substitution into query template HTML. Query results are also accessible via WMS GetFeatureInfo calls, and from MapScript. The following example shows executing a feature query from Python MapScript and fetching back the results:
map = mapscript.Map('rquery.map')
layer = map.getLayer(0)
pnt = mapscript.Point()
pnt.x = 440780
pnt.y = 3751260
layer.queryByPoint( map, pnt, mapscript.MS_MULTIPLE, 180.0 )
layer.open()
for i in range(1000):
result = layer.getResult( i )
if result is None:
break
print '(%d,%d)' % (result.shapeindex, result.tileindex)
s = layer.getShape( result.shapeindex, result.tileindex )
for i in range(layer.numitems):
print '%s: %s' % (layer.getItem(i), s.getValue(i))
layer.close()
This following is a simple example query TEMPLATE file. The raster pixel attributes will be substituted in before the query result is returned to the user as HTML.
Pixel:<br>
values=[value_list]<br>
value_0=[value_0]<br>
value_1=[value_1]<br>
value_2=[value_2]<br>
RGB = [red],[green],[blue]<p>
Class = [class]<br>
Internally raster query results are essentially treated as a set of temporary features cached in RAM. Issuing a new query operation clears the existing query cache on the layer. The transitory in-memory representation of raster query results is also responsible for the inability to save raster query results since saved query results normally only contain the feature ids, not the entire features. Some addition information is available in the RasterQuery Wiki topic.
The following operations use GDAL commandline utilities, some of which are python scripts. They are generally available on any GDAL installation with python support.
The TIFF and Erdas Imagine formats support internal tiling within files, and will generally give better display speed for local map requests from large images. To produce a GeoTIFF file in internally tiled format using the TILED=YES creation option with the gdal_translate utility:
gdal_translate -co TILED=YES original.tif tiled.tif
Erdas Imagine (HFA) files are always tiled, and can be larger than 4GB (the GeoTIFF limit). Use a command like the following to translate a raster to Imagine format:
gdal_translate -of HFA original.tif tiled.img
Rendering and returning 24bit images (especially as PNG) can be quite expensive in render/compress time and bandwidth. Pre-reducing raster data to 8bit can save disk space, processing time, and bandwidth. However, such a color reduction also implicitly reduces the quality of the resulting map. The color reduction can be done on the fly by MapServer but this requires even more processing. A faster approach is to pre-reduce the colors of 24bit imagery to 8bit. This can be accomplished with the GDAL rgb2pct.py script like this:
rgb2pct.py original.tif 8bit.tif
By default images will be reduced to 256 colors but this can mean there are not enough colors to render other colors in the map. So it may be desired to reduce to even less colors:
rgb2pct.py -n 200 original.tif 8bit.tif
Downsampling to 8bit should be done before internal tiling and overview building. The rgb2pct.py script tries to compute an optimal color table for a given image, and then uses error diffusion during the 24bit to 8bit reduction. Other packages (such as ImageMagick or Photoshop) may have alternative color reduction algorithms that are more appropriate for some uses.
Most GDAL supported raster formats can have overviews pre-built using the gdaladdo utility. However, a few formats, such as JPEG2000, MrSID, and ECW already contain implicit overviews in the format themselves and will not generally benefit from external overviews. For other formats (such as GeoTIFF, and Erdas Imagine format) use a command like the following to build overviews:
gdaladdo tile.tif 2 4 8 16 32 64 128
The above would build overviews at x2 through x128 decimation levels. By default it uses “nearest neighbour” downsampling. That is one of the pixels in the input downsampled area is selected for each output pixel. For some kinds of data averaging can give much smoother overview results, as might be generated with this command:
gdaladdo -r average tiled.tif 2 4 8 16 32 64 128
Note that overview building should be done after translating to a final format. Overviews are lost in format conversions using gdal_translate. Also, nothing special needs to be done to make MapServer use GDAL generated overviews. They are automatically picked up by GDAL when mapserver requests a reduced resolution map.
When working with large collections of raster files using a MapServer tileindex, it is desirable to build reduced resolution overview layers that kick in at high scales (using MINSCALE/MAXSCALE to control which layer activates). Preparing the overviews can be a somewhat complex process. One approach is to use the gdal_merge.py script to downsample and mosaic all the images. For instance if we want to produce an overview of many 1meter ortho photos with 250 meter pixels we might do something like:
gdal_merge.py -o overview.tif -ps 250 250 ortho_*.tif
The gdal_merge.py utility suffers from a variety of issues, including no support for different resampling kernels. With GDAL 1.3.2 or later it should be able to accomplish something similar with the more flexible gdalwarp utility:
gdalwarp -rc -tr 250 250 ortho_*.tif overview.tif
In some cases the easiest way of generating an overview is to let MapServer do it using the shp2img utility. For instance if the tileindex layer is called ‘’orthos’’ we could do something like:
shp2img -m ortho.map -l orthos -o overview.png
Note that the overview will be generated with the extents and size in the .map file, so it may be necessary to temporarily adjust the map extents and size values to match the raster extents and the desired output size. Also, if using this method, don’t leave large files in PNG (or GIF or JPEG) format as they are slow formats to extract subareas from.
World files are a simple mechanism for associating georeferencing (world coordinates) information with raster files. ESRI was the first company to propagate the use of world files, and they often used with TIFF instead of embedding georeferencing information in the file itself.
The world file contents look like the following. The first coefficient is the X pixel size. The second and third are rotational/shear coefficients (and should normally be 0.0). The fourth is the Y pixel size, normally negative indicating that Y decreases as you move down from the top left origin. The final two values are the X and Y location of the center of the top left pixel. This example is for an image with a 2m x 2m pixel size, and a top left origin at (356800E, 5767999N):
2
0.0000000000
0.0000000000
-2
356800.00
5767999.00
The name of the world file is based on the file it relates to. For instance, the world file for aerial.tif might be aerial.tfw. Conventions vary for appropriate endings, but with MapServer the extension .wld is always OK for world files.