[Mapserver-users] Point in polygon ?

[Eric Bridger]

--=-UxW5cG3mjojAWQSuwjrO
Content-Type: text/plain
Content-Transfer-Encoding: 7bit

On Tue, 2003-10-28 at 17:17, Stephen Clark wrote:
> Yes that would be nice.
> 
> thanks,
> Stephen
> 

The attached file contains Perl fragments and a map file fragment, just
to demonstrate my approach to determining which of a set of points is
inside of the polygons a user clicked.  This is for Mapserver 3.6.

HTH,
Eric




--=-UxW5cG3mjojAWQSuwjrO
Content-Disposition: attachment; filename=polyquery.txt
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain; name=polyquery.txt; charset=ISO-8859-15

# Perl mapscript for Mapserver 3.6
# Determine which of a set of points are contained in the polygons the user=
 clicked on.

	use CGI ":cgi";
	# Map Server
	use mapscript;
	use XBase;

# A hash of points. we want to determine which of these points are in the p=
olygons clicked by the user.
my %points =3D (
	ONE =3D>{'longitude' =3D> -67.0173,
			'latitude'  =3D> 44.8911,
			},
	TWO =3D>{'longitude' =3D>  -66.0146,
			'latitude' =3D> 45.2045,
			},
	THREE =3D>{'longitude' =3D> -68.3578,
			'latitude' =3D> 43.7148,
			},
	FOUR =3D>	{'longitude' =3D> -66.5528,
			'latitude' =3D> 43.6243,
			},
	FIVE =3D>{'longitude' =3D> -68.9983,
			'latitude' =3D> 44.0555,
			},
	SIX =3D>{'longitude' =3D> -67.8800,
			'latitude' =3D> 43.4900,
			},
	SEVEN =3D>{'longitude' =3D> -70.5665,
			'latitude' =3D> 42.5185,
			},
);

	# form input
	my $q =3D new CGI;

	my $map =3D new mapscript::mapObj('zone.map');

	my ($lon, $lat) =3D get_click($q, $map);

	my $click_pt =3D new mapscript::pointObj();
	$click_pt->{x} =3D $lon;
	$click_pt->{y} =3D $lat;

# this really should be a sub-routine.
	my $layerObj =3D $map->getLayerByName('zones');
	$layerObj->{status} =3D $mapscript::MS_ON;

	# 1) Get all the polygons the user clicked.=20
	$layerObj->queryByPoint($map,$clk_point,$mapscript::MS_MULTIPLE,undef);

	# Open the dbf file. to get the polygon id. This is optional
	my $dbfile =3D $layerObj->{data};
	$dbfile .=3D '.dbf';
	my $hDBF =3D new XBase "$dbfile" or warn XBase->errstr . "\n";

	# 2) Open the shape file to get the shape polygon. must do after each quer=
yByPoint

	if( $layerObj->open($layerObj->{data}) ){
		warn "Could not open shapefile: $layerObj->{data}\n";
	}

	# the results go here
	my @points_found =3D ();
	my @polys_found =3D ();

	my $results =3D $layerObj->{resultcache};
	my $num_results =3D $results->{numresults};
	# for each polygon in which the user clicked.
	# Usually only 1 polygon has been clicked, but not always. Some of
	# our zones overlap.
	for my $i (0..$num_results-1){

		my $rslt =3D $layerObj->getResult($i);
		# shpidx is the record number in the polygon shape dbf file.
		my $shpidx =3D $rslt->{shapeindex};

		# Get the polygon id.  This part is optional, but perhaps we want a list =
of
		# points and the polygon ids they were in.
  		my @row =3D $hDBF->get_record_nf($shpidx, $fld) or warn $hDBF->errstr .=
 "\n";
		my $poly_id =3D $row[1];
		# get the shape polygon
  		my $shape =3D new mapscript::shapeObj(-1);
		if($layerObj->getShape($shape, undef, $shpidx)){
			warn "Error in layer->getShape(): $shpidx\n";
			next;=20
		}

		# 4) get the list of points contained in this polygon.
	=09
		my $point =3D new mapscript::pointObj(); # we re-use this point
		foreach my $point_id ( keys %points) {
			my $lat =3D sprintf("%.4f", $points{$point_id}{latitude});
			my $lon =3D sprintf("%.4f", $points{$point_id}{longitude});
			$point->{x} =3D $lon;
			$point->{y} =3D $lat;
			# Does the shape contain the point?
			my $ret =3D $shape->contains($point);
			if($ret =3D=3D 1){ # yes
				push @points_found, $point_id;
				push @polys_found, $poly_id;
			}

		} # end foreach point

	} # end for num_results

	$layerObj->close();

	$hDBF->close();


############################################################
# SUBROUTINES
############################################################
# translate user click point (img.x and img.y) to lat/long
sub get_click {
	my ($q, $map) =3D @_;
	my @current_extent =3D ();
	my($x_map, $y_map) =3D (0,0);

	if($q->param('img.x') && $q->param('imgext')) { # Make sure we got a click
		@current_extent =3D split(' ', $q->param('imgext'));
		my $click_x =3D $q->param('img.x');
		my $click_y =3D $q->param('img.y');
		my $x_pct =3D ($click_x / $map->{width});
		my $y_pct =3D 1 - ($click_y / $map->{height});

		$x_map =3D $current_extent[0] + ( ($current_extent[2] - $current_extent[0=
]) * $x_pct);
		$y_map =3D $current_extent[1] + ( ($current_extent[3] - $current_extent[1=
]) * $y_pct);
	}
	return ($x_map, $y_map);
}

############################################################
# MAP FILE FRAGMENT
############################################################

# zones
LAYER
  NAME "zones"
  TYPE POLYGON
  STATUS OFF=20
  DATA "../../data/mapserver/vector/zones"
  TEMPLATE 'bogus.html'
  TOLERANCE 1
  CLASS
    NAME "Zone"
    OUTLINECOLOR 0 0 255=20
    # to get thicker lines
    SYMBOL 'circle'
    SIZE 2
  END # class
  PROJECTION
    "proj=3Dlatlong"
  END
END # end zone layer

--=-UxW5cG3mjojAWQSuwjrO--