#!/usr/bin/perl

use strict;
use mapscript;

my $pointx = 429960.850000; 
my $pointy = 3109004.900000;
#my $extents = "429960.56 3109004.3 429962.39 3109006.1";  # all pixels
my $extents = "429960.732 3109004.790 429961.255 3109005.026";  # just two pixels
my($minx,$miny,$maxx,$maxy) = split(/ /, $extents);

my($mapcode,$map,$qshp,$line,$point,$rect);

$mapcode = join('',<DATA>);
    
my $mapfile_scratch = "./work_map.map";
open(MAP,">$mapfile_scratch") || "Couldn't create temp mapfile $mapfile_scratch:$!";
print MAP $mapcode;
close MAP;

$map = new mapscript::mapObj($mapfile_scratch); 

########################
print "-" x 60, "\n";
print "Query by point\n";

$point = new mapscript::pointObj();
$point->{x} = $pointx;
$point->{y} = $pointy;
my $ret = $map->queryByPoint($point, $mapscript::MS_MULTIPLE, 0.0);

print_results($map);

########################
print "-" x 60, "\n";
print "Query by rectangle\n";

$map = new mapscript::mapObj($mapfile_scratch); 
my $rect = new mapscript::rectObj( $minx, $miny, $maxx, $maxy );
my $ret = $map->queryByRect($rect);

print_results($map);

###############
print "-" x 60, "\n";
print "Query by shape\n";

$map = new mapscript::mapObj($mapfile_scratch); 

$line = new mapscript::lineObj();
$point = new mapscript::pointObj();
$point->{x} = $minx;
$point->{y} = $miny;
$line->add($point);
$point = new mapscript::pointObj();
$point->{x} = $maxx;
$point->{y} = $miny;
$line->add($point);
$point = new mapscript::pointObj();
$point->{x} = $maxx;
$point->{y} = $maxy;
$line->add($point);
$point = new mapscript::pointObj();
$point->{x} = $minx;
$point->{y} = $maxy;
$line->add($point);
$point = new mapscript::pointObj();
$point->{x} = $minx;
$point->{y} = $miny;
$line->add($point);
$qshp = new mapscript::shapeObj($mapscript::MS_SHAPE_POLYGON);
$qshp->add($line);
my $ret = $map->queryByShape($qshp);

print_results($map);

exit;



#################################################
sub print_results
{
    my $map = shift;
    
    my $numlayers = $map->{numlayers};
    
    my $layer = $map->getLayer(0);
    
    my $numresults = $layer->getNumResults();
    
    $ret = $layer->open();
    if ($ret == $mapscript::MS_FAILURE)
    {
        die "Layer open failed\n";
    }
    
    for(my $i=0;$i<$numresults; $i++)
    {
        my $record = $layer->getResult($i);
        
        my $shapeidx = $record->{shapeindex};
        my $tileidx  = $record->{tileindex};
        
        my $shp = $layer->getFeature($shapeidx,$tileidx);
        
        my $numitems = $layer->{numitems};
        my %rec = ();
        for (my $i=0; $i<$numitems; $i++)
        {
            my $item = $layer->getItem($i);
            my $value = $shp->getValue($i);
            
            $rec{$item} = $value;
            
            #print "item=$item, value=$value\n";
        }
        
        my $ptr = "";
        $ptr = "<---" if $rec{'x'} == $pointx  && $rec{'y'} == $pointy;
        
        printf "%15.6f %15.6f %4d %4d %4d %s\n", $rec{'x'}, $rec{'y'}, $rec{'red'}, $rec{'green'}, $rec{'blue'}, $ptr;
        
        $shp = undef;  # try to free it
    }
}



__END__
# Start of map file
#
NAME RQ_TEST

SIZE 200 200
UNITS METERS

LAYER
  STATUS DEFAULT
  DATA b008_9929_26917.tif
  TEMPLATE "x"
  TOLERANCE 0
  TOLERANCEUNITS meters
  #NAME "RASTER"
  TYPE RASTER
END
END # Map File
