fastest way to convert lon-lat to points

[P Kishor]

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


-- 
Puneet Kishor http://punkish.eidesis.org/
Nelson Inst. for Env. Studies, UW-Madison http://www.ies.wisc.edu/
Open Source Geospatial Foundation https://edu.osgeo.org/