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