Re: [UMN_MAPSERVER-USERS] fastest way to convert lon-lat to points

[Manfred Meier]

Hello,

i have had some tasks with many points and polylines (streets) too, but 
not so many points. But the execution time in pure perl had been too 
long, therefore I used some algorithms in C and used swig for the 
coupling. That is very fine.

For the pre-searching, which points are in a given rectangle, i use a 
range searching (tree2D from the book from Sedgewick) in C with swig.

All this I use on Linux, so I think my using C and swig is a lot easier 
than on windows.

Manfred


P Kishor schrieb:
On 7/20/06, Steve Lime <[EMAIL PROTECTED]> wrote:

Puneet: Why is speed such a big deal or will you be doing this often?


I will be doing this a few times a year (so, often, although not
necessarily often enough), but this portion is only a part of a longer
process. So, I have to optimize each step as much as I can for overall
better rate of return.


Anyway, if it were me I'd just write a quick MapScript script using
perl. Create a new shapefileObj, a new XBase table, open the file and
start looping. The loop itself should have like 3 lines of code. The
whole script should be about 15...



yup, here is what I have right now (using Geo::ShapeFile) --

foreach point
  foreach polygon
     pointIsInPoly($point, $poly)


pointIsInPoly {
 my ($x, $y) = @$p;  # point being tested
 my @xy = @$poly; # the poly being tested against
 my $n = @xy / 2;                    # Number of points in polygon.
 my @i = map {2 * $_} 0 .. (@xy/2);  # The even indices of @xy.
 my @x = map {$xy[$_ ]} @i;          # Even indices: x-coordinates.
 my @y = map {$xy[$_ + 1]} @i;       # Odd indices: y-coordinates.
 my ($i, $j);                        # Indices.
 my $side = 0;                       # 0 = outside, 1 = inside.
 for ($i = 0, $j = $n - 1 ; $i < $n; $j = $i++) {
   if (
       (
         # If the y is between the (y-) borders ...
         (($y[$i] <= $y) && ($y < $y[$j])) ||
         (($y[$j] <= $y) && ($y < $y[$i]))
       )
       and
         # ...the (x,y) to infinity line crosses the edge
         # from the ith point to the jth point...
         ($x
           <
           ($x[$j] - $x[$i] ) *
           ($y - $y[$i]) / ($y[$j] - $y[$i]) + $x[$i] )) {
        $side = not $side; # Jump the fence.
   }
 }
 return $side ? 1 : 0;
}


>>> P Kishor <[EMAIL PROTECTED]> 7/20/2006 7:19:06 AM >>>
Folks,

I have seen this question asked in the past, and have a similar
problem. I want to convert many lon-lat to a shapefile of points
quickly. The lon-lat are in a csv text file (values.txt) like so

foo,bar,baz,long,lat,qux
1, bar1, baz1, -87.796341, 41.907504, s10
...
...

So, I installed FW_Tools 1.05, created a DSN out of the text file
(this is and has to be on a Windows box), and ogrinfo happily gave me
information that I wanted, but with a twist...

OGRFeature(values.txt):1056
  foo (Integer) = 1
  bar(String) = ba1
  baz(String) = baz1
  long (Real) = -87.796341
  lat (Real) = 41.907504
  qux(String) = 10.0000

So, what's with interpreting qux? It converted my "s10" into a string
of value "10.0000".

Anyway, then I created an OVF file (values.ovf) like so

<OGRVRTDataSource>
  <OGRVRTLayer  name="values">
  <SrcDataSource>ODBC:values</SrcDataSource>
  <SrcLayer>values.txt</SrcLayer>
  <GeometryField encoding="PointFromColumns" x="Long" y="Lat"/>
  <GeometryType>wkbPoint</GeometryType>
  </OGRVRTLayer>
</OGRVRTDataSource>

well,

ogr2ogr -f "ESRI Shapefile" . values.ovf values

but while that seemed to do something without any error, I got no
output.

So, is what I am doing the right way of going about this? Is there a
better or another way to try out and benchmark? Oh, did I mention? --
The values.txt file has about 5.25 million rows, so speed is
important.