Trying to use ESRI shapefiles

johnmc has asked for the wisdom of the Perl Monks concerning the following question:

Greetings!
(1st time poster...)
(2nd update to the 1st post - forgive my ignorance...) I am trying to compare two shapefiles (one with a countywide zoning map, the other parcel points) to identify the zoning area parcels fall into.

Trying to use the Math::Geometry::Planar; to load the polygon for zoning and the $polygon->isinside($point) for existence check.

In the code below, I load the zoning shapes points into an array POLYPTS, then reading the parcel points I want to check which zoning area the parcel falls into.

I am not handling the polypts correctly for the creation of the polygon in the function.

Any and all comments will be greatly appreciated (even if you want to kick a rookie around a bit...)

 use Geo::ShapeFile;
 use Math::Geometry::Planar;
 
 
    $zf = new Geo::ShapeFile("PLA_CDD_Zoning_1006");
    $apn = new Geo::ShapeFile("CAD_AO_ParcelPoints_113006") ;
#    $shapefile = new Geo::ShapeFile("CAD_AO_Parcel_Polygon_1006") ;

    print "Load the Zoning Polygons\n" ;
    undef @zontext ;
    for $zidx (1 .. $zf->shapes()) {
        $zshape = $zf->get_shp_record($zidx);
        %db = $zf->get_dbf_record($zidx);

        $zoning = $db{"ZONING"} ;
        ($work = $db{"OVERLAY"}) =~ s/\_deleted//g ;
        $work =~ tr/\t\n\r//d ;
        $overlay = $work ;
        push @zontext, $zoning.$overlay ;
        undef @pts1 ;
        for $numprt (1 .. $zshape->num_parts) {
                @parts = $zshape->get_segments($numprt);
            $arrayline = '[' ;
                foreach $partarray (@parts) {
                $j = 0 ;
                foreach $partpoint (@$partarray) {
                        ($j1,$xval,$j2,$yval) = split/\=|\,|\)/,$partp
+oint ;
                        $arrayline .= "," if (++$j>1) ;
                    $arrayline .= '[' ;
                        $arrayline .= $xval .','. $yval.']' ;
                 }
             }
             $arrayline .= ']' ;
             eval '$polypts[$zcnt] = $arrayline'; warn $@ if $@;
            ++$zcnt ;
              }
          }

     $outfn = "CCZone.txt" ;
    open (Outfile, ">$outfn") || die "Cannot open output file $outfn" 
+;
    $| = 1, select $_ for select Outfile ;    
    
# need to collect all shapes within all zoning area
    $i = 0 ;
    print "Zoning file version: " ;
    print $zf->file_version(),"\n" ;
    print "Zoning file records: " ;
    print $zf->records(),"\n" ;
    print "Zoning file shapes: " ;
    print $zf->shapes(),"\n" ;
#    $temp = <STDIN> ;


    print "Parcel file shapes: " ;
    print $apn->shapes(),"\n" ;

    print "Main processing...\n" ;
    for $parcels (1 .. $apn->shapes()) {
        
            $pshape = $apn->get_shp_record($parcels);
        $ppoint = $pshape->get_part(1) ;
        ($j1,$px,$j2,$py) = split/\=|\,|\)/,$ppoint ;
            %pids = $apn->get_dbf_record($parcels);

#print $pids{"APN"},"\t" ;
#print $px,",",$py,"\n" ;
            
            $thispoint = [$px, $py] ;
            
        print $parcels,"\n" if ($parcels%10==0) ;
        $zcnt = $zhit = 0 ;
        for $thispoly (@polypts) {
print $thispoly,"\n" ;
# how do I get the string to be an array reference???????
            $polyg1 = Math::Geometry::Planar->new; # creates a new pol
+ygon object;
            $polyg1->points($thispoly) ;
            if ($polyg1->isinside($thispoint)) {
                printf ("%s\t%s\n",$pids{"APN"},$zontext[$zcnt]) 
                    if ($pids{"APN"}) ;
                printf (Outfile "%s\t%s\n",$pids{"APN"},$zontext[$zcnt
+]) 
                    if ($pids{"APN"}) ;
                ++$zhit ;
            }
            ++$zcnt ;
        }
        if (!$zhit) {
             print "No:",$pids{"APN"}," $px, $py\n" ;
        }
    }

exit ;
[download]

Results:

Load the Zoning Polygons
Zoning file version:
Zoning file records: 1019
Zoning file shapes: 1019
Parcel file shapes: 350660
Main processing...
[[6233436.1909606,2194112.25395507],[6233761.94604823,2194103.2760573]
+[6233761.9
4604823,2194103.2760573],[6233762.83704991,2194253.27401428][6233762.8
+3704991,21
94253.27401428],[6233862.31905238,2194250.97103715][6233862.31905238,2
+194250.971
03715],[6233862.11796135,2194125.32895491][6233862.11796135,2194125.32
+895491],[6
233853.10995097,2194050.84399622][6233853.10995097,2194050.84399622],[
+6233913.26
002859,2194048.73994099][6233913.26002859,2194048.73994099],[6233915.6
+0294322,21
93851.74997113][6233915.60294322,2193851.74997113],[6233893.64995314,2
+193869.222
05086][6233893.64995314,2193869.22205086],[6233850.98701156,2193903.17
+798063][62
33850.98701156,2193903.17798063],[6233844.27100507,2193908.52296993][6
+233844.271
00507,2193908.52296993],[6233815.7489989,2193931.22405435][6233815.748
+9989,21939
31.22405435],[6233789.86006027,2193951.83103602][6233789.86006027,2193
+951.831036
02],[6233674.96001573,2193955.16103198][6233674.96001573,2193955.16103
+198],[6233
546.54502092,2193958.87802609][6233546.54502092,2193958.87802609],[623
+3525.83404
865,2193959.47900245][6233525.83404865,2193959.47900245],[6233482.0790
+1437,21939
60.74602904][6233482.07901437,2193960.74602904],[6233429.35195173,2193
+962.273989
1][6233429.35195173,2193962.2739891],[6233404.12906333,2193963.0039648
+4][6233404
.12906333,2193963.00396484],[6233378.90298503,2193963.73700288][623337
+8.90298503
,2193963.73700288],[6233328.37299498,2193965.19797513][6233328.3729949
+8,2193965.
19797513],[6233360.91503536,2194114.33006681][6233360.91503536,2194114
+.33006681]
,[6233361.95698291,2194213.73602217][6233361.95698291,2194213.73602217
+],[6233362
.47195967,2194262.98903315][6233362.47195967,2194262.98903315],[623343
+7.67196539
,2194261.58598941][6233437.67196539,2194261.58598941],[6233436.1909606
+,2194112.2
5395507]]
Can't use string ("[[6233436.1909606,2194112.253955") as an ARRAY ref 
+while "str
ict refs" in use at C:/Perl/site/lib/Math/Geometry/Planar.pm line 923.
[download]

Thanks!

Comment on Trying to use ESRI shapefiles Select or Download Code

Replies are listed 'Best First'.
Re: Trying to use ESRI shapefiles by Albannach (Monsignor) on Apr 03, 2007 at 21:41 UTC
Are you certain your sample point is inside the polygon? In a case like this it is useful to provide a simple case to test your logic: `use strict; use warnings; use Math::Geometry::Planar; my $polygon = Math::Geometry::Planar->new; my $points = [[0,0], [1,0], [1,1], [0,1], ]; $polygon->points($points) ; print "Perimeter:",$polygon->perimeter,"\n" ; print "Area: ",$polygon->area,"\n" ; for my $point ( [.725,.396], [2,3], [1,1] ) { if ($polygon->isinside($point)) { print "point @{$point} is inside\n" ; }else{ print "point @{$point} is outside\n" ; } }` [download] produces the expected `Perimeter:4 Area: 1 point 0.725 0.396 is inside point 2 3 is outside point 1 1 is outside` [download] -- I'd like to be able to assign to an luser	[reply] [d/l] [select]
Re: Trying to use ESRI shapefiles by roboticus (Chancellor) on Apr 03, 2007 at 21:28 UTC
You didn't provide enough code to tell what the problem would be. Did you give `$polygon->points()` what it was expecting? ...roboticus	[reply] [d/l]
Re: Trying to use ESRI shapefiles by thezip (Vicar) on Apr 03, 2007 at 21:34 UTC
Could you please describe how it is not working? ie. What indication do you see? Where do you want them* to go today?*	[reply]
Re: Trying to use ESRI shapefiles by thezip (Vicar) on Apr 03, 2007 at 21:43 UTC
BTW, if you reduce your data set to a trivial data set, and then test for a point that must be inside, does this also produce the same symptom? Where do you want them* to go today?*	[reply]
Re^2: Trying to use ESRI shapefiles by johnmc (Novice) on Apr 03, 2007 at 21:58 UTC
Thanks for the replies - I have restricted the data sets and the behaviour does work nominally. I will run the larger data sets through and see if I can get matching results.	[reply]


"be consistent"
	PerlMonks