[mapserver-users] raster reprojection

[Kent S Ridl]

greetings list...

i've been sifting through the message archives for the last couple days,
so i know that this topic has been beat to death many times... please
forgive me...

my script was working joyously before i dove into reprojecting anything;
everything (tifs, shapefiles) displayed as it should.  when i began
experiementing with reprojection, i got the shapefiles to do so with
little problem, but the rasters won't.  i'm sure of my extents, the
original projections of my test images, etc.  i can't think of anything
else that could be screwing it up, so a vector from anyone out there would
be most helpful.

i'm running MS 3.5 (compiled --without-tiff), using gdal, perl
mapscript and rosa applet, using world files for my tifs... i've attached
my script if anyone wants to look at it.  and if anyone needs more info,
just let me know and i'll get it to you right away.

thanks in advance for your help!
kent
-------------- next part --------------
#!/usr/bin/perl -T

use CGI qw(:standard :html escape);
use mapscript;

# Image paths - absolute and relative
# First 2 are where the mapscript-generated images are stored before being displayed
# Last one is where the input images are taken from
my $img_path = '/www/html/testdata/tmp/';
my $img_rel_path = '/testdata/tmp/';
my $data_path = '/www/html/testdata/';

# Input images
#my $nd = 'nd37wk_ndvi2001.tif';
#my $sd = 'sd37wk_ndvi2001.tif';
#my $mn = 'mn37wk_ndvi2001.tif';
#my $mt = 'mt37wk_ndvi2001.tif';
#my $id = 'id37wk_ndvi2001.tif';
#my $wy = 'wy37wk_ndvi2001.tif';
#my $jpg = '3027d.jpg';
my $refbase = 'refbase.png';
my $tif1 = '3027.tif';
my $tif2 = '3028.tif';
#my $test_tif = '3028_082400.tif';

# Output image and reference map image
# Each new image generated has a timestamp in its name
#   Required for mapserver to display "refreshed" images after a pan/zoom action
my $mapimg = 'testimg'.time().'.png';
my $refmap = 'ref'.time().'.png';

# Uncomment next statement and comment out next 3 if you want to execute this script from the command line
#print "Content type: text/html\n\n";
my $r = shift;
$r->content_type("text/html");
$r->send_http_header;

#---------------------------------------------------------------
# Subroutine set_extent() will calculate a new image's extent.
#   Variables are read in from the HTML form via a param() call.
#---------------------------------------------------------------
sub set_extent() {
  # Variables to hold coordinate conversion values/parameters and new extent
  my ($zoom, at imgext);
  my ($x,$y);
  my ($cx,$cy);
  my ($dx,$dy);

  # If old image extent saved in the hidden form variable "extent", then calculate a new extent
  # An interactive interface is assumed at this point
  # Calculation of new image extents is required to provide the pan/zoom effect
  if(param('extent')) {

    # Parse current extent
    # Resulting array values:
    #   imgext[0] = old minx
    #   imgext[1] = old miny
    #   imgext[2] = old maxx
    #   imgext[3] = old maxy
    @imgext = split(' ',param('extent'));

    # Parse the mouse input coordinates
    # Coordinates are recived in the following forms:
    #   x1,xy;x2,y2 -> opposing corner coordinates of selection box for zoom in
    #   x,y -> single-point mouse click - for map recentering
    @mouse_coords = split(/[,;]/,param('mouse_coords'));

    # Determine which action to take - zoom in, zoom out, recenter
    # Set toolbar buttons and zoom factor accordingly
    if (param('action') == 1) { # zoom in
      $last_cmd = "zin"; # preserve the pressed button on the toolbar
      $dx = $mouse_coords[2] - $mouse_coords[0]; # width of selection box
      $dy = $mouse_coords[3] - $mouse_coords[1]; # height of selection box
      if ($dx <= 0) {$dx = 1;} # prevent divide-by-0 error
      $zoom = $map->{width} / $dx; # zoom factor is map width proportionate - determined by width of selection box
      $x = $mouse_coords[0] + ($dx / 2); # x-coord of box center - offset from upper left corner of box
      $y = $mouse_coords[1] + ($dy / 2); # y-coord of box center - offset from upper left corner of box
    } # end if
    elsif (param('action') == 0) { # recenter the map
      $last_cmd = "recenter"; # preserve the pressed button on the toolbar
      $zoom = 1;
      $x = $mouse_coords[0];
      $y = $mouse_coords[1];
    } #end elsif
    elsif (param('action') == -1) { # zoom out
      $zoom = -2;

      # Use the map center coordinate as center of zoom out function
      $x = $map->{width} / 2;
      $y = $map->{height} / 2;
    } #end elsif

    $zoom = 1 / abs($zoom) if $zoom < 0; # tweak zoom factor for calculations

    # Calculate cellsize in x and y
    $cx = ($imgext[2] - $imgext[0]) / ($map->{width}); # scale cell by extent width w/ respect to map width
    $cy = ($imgext[3] - $imgext[1]) / ($map->{height}); # scale cell by extent height w/ respect to map height

    # Convert mouse coordinate values to map coordinates w/ respect to map extents ($ci * $i)
    # Offset scaled coordinates from upper left corner of previous image
    $x = $imgext[0] + $cx * $x;
    $y = $imgext[3] - $cy * $y;

    # Calculate new extent - not sure about the reasoning behind the math, but it works
    $map->{extent}->{minx} = $x - .5 * (($imgext[2] - $imgext[0]) / $zoom);
    $map->{extent}->{miny} = $y - .5 * (($imgext[3] - $imgext[1]) / $zoom);
    $map->{extent}->{maxx} = $x + .5 * (($imgext[2] - $imgext[0]) / $zoom);
    $map->{extent}->{maxy} = $y + .5 * (($imgext[3] - $imgext[1]) / $zoom);
  } # end if
} # end sub set_extent()

#-------------------
# Construct the map.
#-------------------

# Map object and parameters
$map = new mapObj(undef) or die("Unable to create mapfile... $!");
$map->{width} = 500;
$map->{height} = 400;
$old_proj = new projectionObj("datum=wgs84"); # old projection system
$map->setProjection("proj=utm,ellps=GRS80,zone=14"); # target map projection
$map->{transparent} = $mapscript::MS_ON;

# Layer for shapefile and parameters
$layer = new layerObj($map) or die("Unable to create new layer... $!");
$layer->{type} = $mapscript::MS_LAYER_POLYGON;
$layer->setProjection("datum=wgs84"); # original shapefile projection

# Class for shapefile
# Used to change outline color, fill color, etc. of the shapefile
$class = new classObj($layer) or die("Unable to create class object... $!");
$class->{outlinecolor} = $map->addColor(0,0,255);

# Layer for raster images and parameters
$layer_zappa = new layerObj($map) or die("No zappa for you!... $!");
$layer_zappa->{type} = $mapscript::MS_LAYER_RASTER;
$layer_zappa->setProjection("datum=wgs84"); # original projection for raster data
$layer_zappa->{status} = $mapscript::MS_ON;
$layer_zappa->{offsite} = 0; # color in this index is drawn transparent

# Open the shapefile
$shapefile = new shapefileObj("/var/www/html/testdata/umac_states",-1) or die("Unable to open shapefile... $!");

# Set original map extents to be the full extent of the entire shapefile
$map->{extent} = $shapefile->{bounds}; # bounding box of shapefile is original map extents
$map->{extent}->project($old_proj,$map->{projection}); # reproject map extents to target projection

# Check if maps original extents need to be restored
if (!param("restore")) { # check if "Restore" button clicked
  &set_extent(); # set appropriate extents if pan/zoom
} # end if

$img = $map->prepareImage(); # prepare an image to be written to

#
# The raster draws will be replaced later on with a loop to draw the appropriate LandSat data
#
# Draw raster data of the UMAC states
# For some reason, a die clause after drawing a raster will crash the program... I don't know why.
$layer_zappa->{data} = $data_path . $tif1;
$layer_zappa->draw($map,$img);
$layer_zappa->{data} = $data_path . $tif2;
$layer_zappa->draw($map,$img);
#$layer_zappa->{data} = $data_path . $test_tif;
#$layer_zappa->draw($map,$img);
#$layer_zappa->{data} = $data_path . $mn;
#$layer_zappa->draw($map,$img); # draw MN raster
#$layer_zappa->{data} = $data_path . $mt;
#$layer_zappa->draw($map,$img); # draw MT raster
#$layer_zappa->{data} = $data_path . $id;
#$layer_zappa->draw($map,$img); # draw ID raster
#$layer_zappa->{data} = $data_path . $wy;
#$layer_zappa->draw($map,$img); # draw WY raster

# Create shape object to hold outline of each state
$shape = new shapeObj(-1) or die("Unable to create shape... $!");

# Draw outlines for UMAC states
# Draw order - MT, ND, SD, WY, ID, MN
# Die clause also crashes the program after drawing a shape
for ($i = 0; $i <= 5; $i++) {
  $shapefile->get($i,$shape); # get a shape from the shapefile, referenced by index number
  $shape->draw($map,$layer,$img); # draw the shape
} # end for

# Create a reference map
# Used to maintain a global sense of location whilst zoomed in
$map->{reference}->{extent} = $shapefile->{bounds}; # bounding box of shapefile is original refmap extents
$map->{reference}->{extent}->project($old_proj,$map->{projection}); # reproject refmap extents to target projection
#$map->{reference}->{color}->{red} = 0;
#$map->{reference}->{color}->{green} = 0;
#$map->{reference}->{color}->{blue} = 255;
#$map->{reference}->{color} = $map->addColor(0,0,255);
$map->{reference}->{outlinecolor}->{red} = 255;
$map->{reference}->{outlinecolor}->{green} = 0;
$map->{reference}->{outlinecolor}->{blue} = 0;
$map->{reference}->{image} = $refbase; # base image for reference map (needs to be rebuilt with target projection)
$map->{reference}->{status} = $mapscript::MS_ON;
$ref = $map->drawReferenceMap(); # draw the reference map

# Save current extent
$imgext = join(' ',$map->{extent}->{minx},$map->{extent}->{miny},$map->{extent}->{maxx},$map->{extent}->{maxy});

$map->save("foo.map"); # create a mapfile
mapscript::msSaveImage($img,$img_path.$mapimg,$mapscript::MS_PNG,1,1,0); # save image to file
mapscript::msFreeImage($img); # free image from memory
mapscript::msSaveImage($ref,$img_path.$refmap,$mapscript::MS_PNG,1,1,0); # save reference image to file
mapscript::msFreeImage($ref); # free reference image from memory

#----------------
# Print out HTML.
#----------------

print << "EOF";
   <html>

   <head>
   <title>Test image print</title>
   </head>

   <body>
   <p><h1>HI AMANDA!!  :-)</h1></p>
   <p><form name="mapserv" action="testrosa.pl" method="GET">
      <table>
        <tr>
	<td>
          <applet name="RosaApplet" archive="rosa/class/rosa.jar" code="Rosa2000" width="$map->{width}" height="$map->{height}" MAYSCRIPT>
	  <param name="img_url" value="$img_rel_path$mapimg">
	  <param name="inp_form_name" value="mapserv">
	  <param name="inp_type_name" value="input_type">
	  <param name="inp_coord_name" value="mouse_coords">
          <!--<param name="bg_color" value="#bcbcbc">-->
	  <param name="tb_buttons" value="zin|zout|recenter|fullextent">
	  <param name="tb_position" value="top">
          <param name="tb_selected_button" value="$last_cmd">
	  <param name="tb_but_zin_img" value="zin1.gif">
	  <param name="tb_but_zin_img_pr" value="zin2.gif">
	  <param name="tb_but_zin_input" value="auto_rect">
	  <param name="tb_but_zin_name" value="action">
	  <param name="tb_but_zin_value" value="1">
          <param name="tb_but_zin_hint" value="Click this button and drag a|box around your area of interest.">
	  <param name="tb_but_zout_img" value="zout1.gif">
	  <param name="tb_but_zout_img_pr" value="zout2.gif">
	  <param name="tb_but_zout_input" value="submit">
	  <param name="tb_but_zout_name" value="action">
	  <param name="tb_but_zout_value" value="-1">
          <param name="tb_but_zout_hint" value="Click to zoom out 2x.">
	  <param name="tb_but_recenter_img" value="recenter1.gif">
	  <param name="tb_but_recenter_img_pr" value="recenter2.gif">
	  <param name="tb_but_recenter_input" value="auto_point">
	  <param name="tb_but_recenter_name" value="action">
	  <param name="tb_but_recenter_value" value="0">
          <param name="tb_but_recenter_hint" value="Click this button and then choose the|map's new center with a mouse click.">
	  <param name="tb_but_fullextent_img" value="fullextent1.gif">
	  <param name="tb_but_fullextent_img_pr" value="fullextent2.gif">
	  <param name="tb_but_fullextent_input" value="submit">
	  <param name="tb_but_fullextent_name" value="restore">
	  <param name="tb_but_fullextent_value" value="restore">
          <param name="tb_but_fullextent_hint" value="Click to restore the map|to it's original state.">
	  </applet>
	</td>

        <td valign="top"><img src="$img_rel_path$refmap" border="0"></td>
        </tr>
      </table>
      <input type=hidden name="extent" value="$imgext">
      <input type=hidden name="input_type" value="">
      <input type=hidden name="mouse_coords" value="">
      <input type=hidden name="action" value="">
      <input type=hidden name="restore" value="">
      </form></p>
   </body>

   </html>
EOF

#----------------------------------
# In case I need any of this again.
#----------------------------------

# print "$mapscript::ms_error->{code} $mapscript::ms_error->{message} $mapscript::ms_error->{routine}\n<p>";

#$aspect = ($map->{extent}->{maxy} - $map->{extent}->{miny}) / ($map->{extent}->{maxx} - $map->{extent}->{minx});
#$map->{reference}->{width} = $map->{width} / 3;
#$map->{reference}->{height} = $map->{reference}->{height} * $aspect;

#$cellx = ($map->{extent}->{maxx} - $map->{extent}->{minx}) / $map->{reference}->{width};
#$celly = ($map->{extent}->{maxy} - $map->{extent}->{miny}) / $map->{reference}->{height};
#if ($cellx > $celly) {
#  $cellsize = $cellx;
#}
#else {
#  $cellsize = $celly;
#}
#$map->{reference}->{extent}->{minx} = $map->{extent}->{minx};
#$map->{reference}->{extent}->{miny} = $map->{extent}->{maxy} - ($cellsize * ($map->{height} + 1));
#$map->{reference}->{extent}->{maxx} = $map->{extent}->{minx} + ($cellsize * ($map->{width} + 1));
#$map->{reference}->{extent}->{maxy} = $map->{extent}->{maxy};