[Mapserver-dev] Rotated Map Rendering

[Frank Warmerdam]

Folks,

Tydac AG would like to be able to render rotated map views with MapServer.
I have prepared the following proposal for changes I would like to make
to make this possible.

The approach is essentially to hold the rotational transformation down
within the projectionObj for the map.  The map will "think" it has a
north up coordinate system still so it can continue to use the mapObj->extents
as it does now in all the rendering code.

Any objections to be pursuing this in the 4.3 tree?

It is hoped that support for non-square pixels will also fall out of this
work, allowing us to complete WMS conformance testing.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent




	MapServer Rotated Map Support Proposal
	======================================


Overall Approach
----------------

We will provide a new directive, and corresponding MapScript calls to set
a rotated extent based on a center point, width/height and rotation angle.

Internally a "fake" rotated coordiante system will be maintained for the
mapObj and all projection routines will know how to apply the extra step to
map from true map projections coordinates to this fake coordinate system. This
fake coordinate system will internally still have a seemingly north-up map.



Caveats and Notes
-----------------

  o The rotation effect will only cause an alteration in function in
    msDrawMap().  It should have no effect on query mode or other operations.

  o The rotated coordinate system should not be used with OGC WMS protocol
    as there is no way to communicate an appropriate BBOX.

  o There is some confusion still in my mind about how what means we will
    have for interacting with the rotation mode via CGI if any.  The
    emphasis is to make it easily used from MapScript and raw .map files
    for now.

  o There may be some confusing issues with scale testing.

  o Text angles will be relative to the rotated map coordinate system.

  o Reference map drawing will not currently be updated to handle proper
    rotated reference area display.


New Map File Directive
----------------------

EXTENT_ROT <center_x> <center_y> <width> <height> <rotation_angle>

<center_x>,<center_y> - center of view window in projected coordinates.
<width><height> - width and height of window in projected units(meters/degrees)
<rotation_angle> - counter-clockwise rotation of view in degrees.


This keyword and the setExtentRot() method in MapScript will set the mapObj
extents rectObj to be exactly the same size as the rotated view, but
unrotated.  This unrotated extent will cause a correct cellsize to be
set, but in fact, will not include quite all the data that will be needed
for the rotated render.


MapScript
---------

New methods:

  map->setExtentRot( center_x, center_y, width, height, rotation_angle )


Documentation
-------------

Updates will be made to the mapfile reference, and PHP and non-PHP MapScript
documentation briefly explaining the new rotation mode.  A mini-HOWTO on
generation of rotated maps will be added to the Wiki.


Data Structure Changes
----------------------

The mapObj will include following extra parameters:

  // EXTENT_ROT
  int  using_ext_rot;
  double ext_rot_center_x;
  double ext_rot_center_y;
  double ext_rot_width;
  double ext_rot_height;
  double ext_rot_rotation_angle;

  // This is the affine transformation from pixel/line coordinates to
  // pixel/line coordinates.
  double extent_geotransform[6];
  double extent_invgeotransform[6];

The projectionObj will include the following extra parameters:

  int    using_geotransform;
  double pixel_to_proj_geotransform[6];
  double proj_to_pixel_geotransform[6];


Computations
------------

mapObj->extent_geotransform[] computation:

   gt[0] = ext_rot_center_x
        - 0.5 * (cos(rot_angle)*ext_rot_width - sin(rot_angle)*ext_rot_height)
   gt[1] = cos(rot_angle) * ext_rot_width / mapObj->width
   gt[2] = sin(rot_angle) * ext_rot_height / mapObject->height
   gt[3] = ext_rot_center_y
        + 0.5 * (sin(rot_angle)*ext_rot_width + cos(rot_angle)*ext_rot_height)
   gt[4] = sin(rot_angle) * ext_rot_width / mapObj->width
   gt[5] = -cos(rot_angle) * ext_rot_height / mapObj->height

mapObj->extent_invgeotransform[] is computed fro extent_geotransform[] using
InvGeoTransform() in mapresample.c (invert affine transform).

Note that this definition of geotransform[] maps from a coordinate system
with x going from 0 to width (in pixels) and y going from height to 0 (in
pixels).   This is a nice well behaved coordinate system from the point
of view of the rendering.  Note that the mapObj->cellsize will be 1.0.

A "rotated map coordinate system" coordinate would be converted to real map
projected coordinates with the transformation:

  x_real = gt[0] + x_fake * gt[1] + y_fake * gt[2]
  y_real = gt[3] + x_fake * gt[4] + y_fake * gt[5]

The reverse just uses the inverted geotransform.


msDrawMap() Alterations
-----------------------

Within msDrawMap() (used for WMS renders, regular cgi renders and full
renders from MapScript) we will take the following steps:

  o Copy the extent_geotransform[] and extent_invgeotransform[] to the
    mapObj->projection pixel_to_proj_geotransform[] and
    proj_to_pixel_geotransform[] fields, and set projection "using_geotransform"
    flag.
  o Update the mapObj->extent to have coordinates that produce the desired
    region in the altered (rotated) coordinate system (set to 0 0 width height).
  o On exit, the mapObj->projection->using_geotransform will be disabled and
    mapObj->extent values will be restored to it's original value.
  o Note that because this is being done at the drawMap() level, explicit
    MapScript drawLayer() calls will ignore our rotation.


mapproject.c Alterations
------------------------

In msProjectPoint() if using_geotransform is TRUE, apply the appropriate
geotransform when going to or from the map's projection.

In msProjectionsDiffer(), ensure that TRUE is returned if a geotransform
is in effect.


Raster Renderer
---------------

The raster renderer does not go through msProjectPoint() and friends in
the resampler.  Extra code will need to be added to apply the geotransform
when enabled.  Note that because the map projection will differ from the
file projection (due to the geotrasform) rotated maps will always be
processed through the mapresample.c code rather than through the more
direct code.  This does imply some performance reduction if reprojection
of the image was not otherwise going to take place.