.. index::
   single: Raster data

.. _raster:

*****************************************************************************
 Raster Data 
*****************************************************************************

:Author:        Frank Warmerdam
:Contact:       warmerdam at pobox.com
:Last Updated: 2013/07/04

.. contents:: Table of Contents
    :depth: 2
    :backlinks: top


Introduction
============

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
:ref:`mapfile`, but does explain the raster specific aspects of mapfiles.

How are rasters added to a Map file?
====================================

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.

.. index::
   single: Raster classification

Classifying Rasters
-------------------

Rasters can be classified in a manner similar to vectors, with a few
exceptions.

.. index::
   single: pixel

.. index::
   single: red

.. index::
   single: green

.. index::
   single: blue

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, equivalent 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).

Classifying Non-8bit Rasters
****************************

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

.. index::
   single: Raster formats

Supported Formats
=================

Since version 6.2, MapServer raster input support is through the GDAL
raster library only.

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 if GDAL support is built into a particular MapServer
executable, use the -v flag 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 6.2.0 OUTPUT=GIF OUTPUT=PNG OUTPUT=JPEG
    SUPPORTS=PROJ SUPPORTS=GD SUPPORTS=AGG SUPPORTS=FREETYPE
    SUPPORTS=CAIRO SUPPORTS=OPENGL SUPPORTS=ICONV SUPPORTS=WMS_SERVER
    SUPPORTS=WMS_CLIENT SUPPORTS=WFS_SERVER SUPPORTS=WFS_CLIENT
    SUPPORTS=WCS_SERVER SUPPORTS=SOS_SERVER SUPPORTS=FASTCGI
    SUPPORTS=THREADS SUPPORTS=GEOS INPUT=JPEG INPUT=POSTGIS INPUT=OGR
    INPUT=GDAL INPUT=SHAPEFILE


    warmerda@gdal2200[18]% gdalinfo --formats
    Supported Formats:
      VRT (rw+v): Virtual Raster
      GTiff (rw+vs): GeoTIFF
      NITF (rw+vs): National Imagery Transmission Format
      RPFTOC (rovs): Raster Product Format TOC format
      ECRGTOC (rovs): ECRG TOC format
      ...

.. index::
   pair: Raster; Tile indexing

Rasters and Tile Indexing
=========================

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 Mapfile 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).

Starting with MapServer 6.4, the files can have different coordinate
system (projection). This requires specifying the :ref:`TILESRS
<tilesrs>` keyword and generating the tileindex with a few additional
options. See :ref:`tileindex_mixed_srs`.

::

  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

Tile Index Notes
-------------------

- The shapefile (index_file) will be created if it doesn't already exist.
- The default tile index field is 'location'.
- Simple rectangular polygons must be generated in the same
  coordinate system as the raster layer. If the files in the tileindex are
  not in the same projection as the raster layer, or are in heterogeneous
  projections, the TILESRS keyword must be specified in the LAYER definition.
  See :ref:`tileindex_mixed_srs`
- Raster filenames will be put in the file exactly as they are specified
  on the commandline.
- Many problems with tile indexes relate to how relative paths in the tile
  index are evaluated.  They should be evaluated relative to the `SHAPEPATH` 
  if one is set, otherwise relative to the tileindex file.  When in doubt 
  absolute paths may avoid path construction problems.

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.

.. seealso::
    :ref:`tileindex`
    

.. index::
   pair: Raster; Warping

Raster Warping
==============

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.

.. index::
   pair: Raster; 24bit

24bit RGB Rendering
===================

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. Rendering of
24bit imagery is supported via the GDAL renderer.

If the output is still 8bit pseudo-colored (the `IMAGEMODE` is PC256 in the
associated `OUTPUTFORMAT` 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.

.. index::
   pair: Raster; Processing directives

Special Processing Directives
=============================

As of MapServer 4.0, the PROCESSING parameter was added to the :ref:`LAYER` of
the :ref:`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.

.. index::
   triple: LAYER; PROCESSING; BANDS

**BANDS=red_or_grey[,green,blue[,alpha]]**
    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"
      
.. index::
   triple: LAYER; PROCESSING; COLOR_MATCH_THRESHOLD

**COLOR_MATCH_THRESHOLD=n**
    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"
  
.. index::
   triple: LAYER; PROCESSING; DITHER

**DITHER=YES**
    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"
      
.. index::
   triple: LAYER; PROCESSING; EXTENT_PRIORITY

**EXTENT_PRIORITY=WORLD**
    Override GDAL with a world file.

    Example:
  
    ::

       PROCESSING "EXTENT_PRIORITY=WORLD"
      
.. index::
   triple: LAYER; PROCESSING; LOAD_FULL_RES_IMAGE

**LOAD_FULL_RES_IMAGE=YES/NO**
    This option affects how image data is loaded for the resampler when
    reprojecting or otherwise going through complex resampling (as opposed 
    to the fast default image decimation code path).  This forces the source
    image to be loaded at full resolution if turned on (default is NO).  This
    helps work around problems with default image resolution selection in 
    when radical warping is being done.  It can result in very slow processing
    if the source image is large.

.. index::
   triple: LAYER; PROCESSING; LOAD_WHOLE_IMAGE

**LOAD_WHOLE_IMAGE=YES/NO**
   This option affects how image data is loaded for the resampler (as above).
   This option, if turned on, will cause the whole source image to be loaded
   and helps make up for problem identifying the area required, usually due
   to radical image reprojection near a dateline or projection "horizon". 
   The default is NO.  Turning this on can dramatically affect rendering 
   performance and memory requirements.

.. index::
   triple: LAYER; PROCESSING; LUT[_n]

**LUT[_n]=<lut_spec>**
    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.

    Starting with MapServer 7.2, the range of input values can go up to 65535.
    The last point of the curve must be explicitly listed, e.g. for a 12-bit
    raster, 4095:255. This extended LUT syntax is only properly supported if
    no SCALE directive is specified, since the purpose of extended LUT is mainly to
    avoid pre-scaling.

    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"
      
.. index::
   triple: LAYER; PROCESSING; OVERSAMPLE_RATIO

**OVERSAMPLE_RATIO=double**
    Default is 2.5.  Rendering time will increase with increasing
    OVERSAMPLE_RATIO.

    Example:
  
    ::
  
        PROCESSING "OVERSAMPLE_RATIO=1.0"
      
.. index::
   triple: LAYER; PROCESSING; RESAMPLE

**RESAMPLE=NEAREST/AVERAGE/BILINEAR**
    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
    :ref:`rfc4`.

    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"

.. index::
   triple: LAYER; PROCESSING; SCALE

**SCALE[_n]=AUTO or min,max**
    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"
      
.. index::
   triple: LAYER; PROCESSING; WORLDFILE

**WORLDFILE=<file>**
    Specifies an alternative world file (for georeferencing).  If a
    path only is specified, the base name of the dataset will be
    appended.  The suffix (`.wld` / `.tfw` / ...) can be omitted.
  
    Example:
  
    ::
  
       PROCESSING "WORLDFILE=/path/"
        or
       PROCESSING "WORLDFILE=/path/file.wld"
        or
       PROCESSING "WORLDFILE=/path/file"
      
.. index::
   pair: Raster; Query

Raster Query
============

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:

.. index::
   pair: Raster query; x

:x: georeferenced X location of pixel.

.. index::
   pair: Raster query; y

:y: georeferenced Y location of pixel.

.. index::
   pair: Raster query; value_list

:value_list: a comma separated list of the values of all selected bands at the
  target pixel. 
  
.. index::
   pair: Raster query; value_n

: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.

.. index::
   pair: Raster query; class

:class:
  Name of the class this pixel is a member of (classified layers only).

.. index::
   pair: Raster query; red

:red: red component of the display color for this pixel.

.. index::
   pair: Raster query; green

:green: green component of the display color for this pixel.

.. index::
   pair: Raster query; blue

: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 a lot 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:

.. code-block:: python

  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
        
    s = layer.getShape( result )
    for i in range(layer.numitems):
      print '%s: %s' % (layer.getItem(i), s.getValue(i))
            
  layer.close()

This following is a simple example query :ref:`TEMPLATE <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.

.. index::
   pair: Raster; Performance tips

Raster Display Performance Tips
===============================

- Build overview levels for large rasters to ensure only a reasonable
  amount of data needs to be touched to display an overview of a large
  layer. Overviews can be implemented as a group of raster layers at
  different resolutions, using `MINSCALEDENOM`, and `MAXSCALEDENOM` to
  control which layers are displayed at different
  resolutions. Another, perhaps easier way, is to build overviews for
  GDAL supported formats using the `gdaladdo`_ utility.

- When using tileindexes to manage many raster files as a single file,
  it is especially important to have an overview layer that kicks in
  at high scales to avoid having to open a large number of raster
  files to fulfill the map request.

- Preprocess RGB images to eightbit with a colormap to reduce the
  amount of data that has to be read, and the amount of computation to
  do on the fly.

- For large images use tiling to reduce the overhead for loading a
  view of a small area. This can be accomplished using the `TILEINDEX`
  mechanism of the mapfile, or by creating a tiled format file
  (ie. TIFF with GDAL).

- Ensure that the image is kept on disk in the most commonly requested
  projection to avoid on-the-fly image warping which is fairly
  expensive.

- If you are getting debug output from MapServer in your web server
  log file, check to see if the message msResampleGDALToMap in effect
  appears. If so, the raster layer is being resampled. If you don't
  think it should be resampled carefully review your map file to
  ensure that the layer projection exactly matches the map projection
  or that the layer has no projection definition.


.. index::
   pair: Raster; Preprocessing

Preprocessing Rasters
=====================

The following operations use GDAL commandline utilities, some of which are
python scripts.  They are generally available on any GDAL installation with
python support.  

.. index::
   pair: Raster; Tiled datasets

Producing Tiled Datasets
------------------------

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


.. index::
   pair: Raster; 8bit

Reducing RGB to 8bit
--------------------

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.


.. index::
   pair: Raster; Internal overviews

Building Internal Overviews
---------------------------

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.


.. index::
   pair: Raster; External overviews

Building External Overviews
---------------------------

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 `MINSCALEDENOM` / `MAXSCALEDENOM`
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.


.. index::
   pair: Raster; Georeferencing

.. index::
   pair: Raster; World files

Georeference with World Files
=============================

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 are 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.


.. #### rST Link Section ####

.. _`supported formats list`: http://www.gdal.org/formats_list.html
.. _`gdaltindex`: http://www.gdal.org/gdal_utilities.html
.. _`FWTools`: http://fwtools.maptools.org/
.. _`MS4W`: http://www.ms4w.com
.. _`OSGeo4W`: http://trac.osgeo.org/osgeo4w/
.. _`RasterQuery`: https://github.com/mapserver/mapserver/wiki/RasterQuery
.. _`gdaladdo`: http://www.gdal.org/gdal_utilities.html
.. _`gdal_translate`: http://www.gdal.org/gdal_utilities.html
