Beefy Boxes and Bandwidth Generously Provided by pair Networks
go ahead... be a heretic
 
PerlMonks  

Trying to use ESRI shapefiles

by johnmc (Novice)
on Apr 03, 2007 at 20:23 UTC ( #608146=perlquestion: print w/ replies, xml ) Need Help??
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 ;


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.
Thanks!

Comment on Trying to use ESRI shapefiles
Select or Download Code
Re: Trying to use ESRI shapefiles
by roboticus (Canon) 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

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?
Re: Trying to use ESRI shapefiles
by Albannach (Prior) 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" ; } }
    produces the expected

    Perimeter:4 Area: 1 point 0.725 0.396 is inside point 2 3 is outside point 1 1 is outside

    --
    I'd like to be able to assign to an luser

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?
      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.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://608146]
Approved by pKai
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others musing on the Monastery: (5)
As of 2014-09-17 02:21 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    How do you remember the number of days in each month?











    Results (56 votes), past polls