fastest way to convert lon-lat to points

Manfred Meier m.meier at SPIEKERMANN.DE
Thu Jul 20 21:33:30 PDT 2006


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 <Steve.Lime at dnr.state.mn.us> 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 <punkish at EIDESIS.ORG> 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.
>>
> 
> 



More information about the MapServer-users mailing list