Beefy Boxes and Bandwidth Generously Provided by pair Networks
Do you know where your variables are?
 
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 studying the Monastery: (4)
As of 2014-08-29 01:20 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The best computer themed movie is:











    Results (275 votes), past polls