Permafacture / Py-Visvalingam-Whyatt

Python implementation of Visvalingam and Whyatt's algorithm for reducing the complexity of poly-lines
44 stars 12 forks source link

IndexError: index out of bounds #5

Closed DemersM closed 9 years ago

DemersM commented 9 years ago

I tried your script on your zoning shapefiles and encounter the following error on the feature with ZONING_ID : 225943.

ERROR:django.request:Internal Server Error: /basqui/layer/shapefile/simplify
Traceback (most recent call last):
File "C:\Python27\lib\site-packages\django\core\handlers\base.py", line 112, in get_response
    response = wrapped_callback(request, *callback_args, **callback_kwargs)
File "C:\Python27\lib\site-packages\django\contrib\auth\decorators.py", line 22, in _wrapped_view
    return view_func(request, *args, **kwargs)
File "C:\basqui\geodjango\shapefile\views.py", line 331, in simplifyVectorLayer
    simplifier = VWSimplifier(pts)
File "C:\basqui\geodjango\shapefile\simplify.py", line 68, in __init__
    self.thresholds = self.build_thresholds()
File "C:\basqui\geodjango\shapefile\simplify.py", line 79, in build_thresholds
    real_areas = array([triangle_area(pts[n-1],pts[n],pts[n+1]) if n not in [0,n max-1] else np.inf for n in range(nmax)])
File "C:\basqui\geodjango\shapefile\simplify.py", line 45, in triangle_area
    return abs(p1[0]*(p2[1]-p3[1])+p2[0]*(p3[1]-p1[1])+p3[0]*(p1[1]-p2[1]))/2.
IndexError: index out of bounds

Note that prior to running the script, the shapefile was transform to EPSG:3857

Permafacture commented 9 years ago

Usage would be like this. Does this work?

simplifiers = []
for geom in layer.get_geoms():
    simplifiers.append(GDALSimplifier(geom))
DemersM commented 9 years ago

Im using a django queryset as input like this:

features = Feature.objects.filter(shapefile__pk=layer_id)
simplifier = []
    for feat in features:
        simplifier.append(VWSimplifier(feat.geom_multipolygon)) #geom_multipolygon is:  models.MultiPolygonField(srid=3857)

the geometries are in epsg:3857. The feature with ZONING_ID : 225943 returns the error above, the other features seems ok and are simplified successfuly

Permafacture commented 9 years ago

VWsimplifier just takes a polyline as a list of points. I added GDALSimplifier to take an OGRGeometry object. I guess I misunderstood your request for gdal support. I'll at the multipolygonfield later today.

DemersM commented 9 years ago

Allright, I understand. I have tested with OGRGeomtry input and GDALSimplifier as you suggested. There was no error. However, I just wander how to recreated OGRGeometry from the numpy array output. Have a suggestion for that?

Permafacture commented 9 years ago

It would be very slow to create those objects just to appease the simplifier. It just turns them into arrays and (optionally) back again anyway. I'm interested in coming up with a solution that fits common workflows.

Quick hack for you would be to make simplifiers for each polyline in the multipolygon, like GDALSimplifier does. I'll checkout the geodjango geometry field objects today and hopefully push a solution.

Permafacture commented 9 years ago

Can you describe your workflow? Are you just simplifying once? Or depending on the request? You are loading into a database?

DemersM commented 9 years ago

The workflow is very basic I guess. All data are in a database.

1- Create a queryset of features 2- Pass all features's geometry field into the simplify script 3- Update the geometry field of each features with the simplified geometry

On Mon, Aug 4, 2014 at 2:00 PM, Elliot Hallmark notifications@github.com wrote:

Can you describe your workflow? Are you just simplifying once? Or depending on the request? You are loading into a database?

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-51094597 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Permafacture commented 9 years ago

I don't have a geodjango database, so no queryset to look at. But, I'm missing something here. A queryset is giving you a multipolygonfield as a result? Shouldn't it give you the value that field contains? Like if I had an IntegerField, the queryset would tell me the value of the integer there, not give me an IntegerField object to work with. The simplifier can't simplify something without having the points to work with.

How do you get the value of the MultiPolygonField when it is a result in a queryset? What is the type of the value you are getting? Or, how do you get the points to do something with them?

DemersM commented 9 years ago

Sorry, I was not clear enough because I thought you were working with geodjango already.

The queryset is a list of feature objects where the geometry are stored in a field (in my case the name of that field is geom_multipolygon). The geometry is stored as a GEOS geometry object (django.contrib.gis.geos https://github.com/django/django/blob/master/django/contrib/gis/geos/geometry.py), The GEOS object has many function to retrieve coordinates like wkt or to transform it to OGRGeometry object.

I can use the GDALSimplifier() with OGR Geometry as input, but probably it would be more straightfoward for django application to allow GEOS Geometry as input and output

I hope that answer your questions.

Thank you for your time and support

-Max Demars

On Mon, Aug 4, 2014 at 4:04 PM, Elliot Hallmark notifications@github.com wrote:

I don't have a geodjango database, so no queryset to look at. But, I'm missing something here. A queryset is giving you a multipolygonfield as a result? Shouldn't it give you the value that field contains? Like if I had an IntegerField, the queryset would tell me the value of the integer there, not give me an IntegerField object to work with. The simplifier can't simplify something without having the points to work with.

How do you get the value of the MultiPolygonField when it is a result in a queryset? What is the type of the value you are getting? Or, how do you get the points to do something with them?

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-51109993 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

DemersM commented 9 years ago

Im thinking about it... It would be easier to skip the GEOS geometry creation by doing a raw sql query on the db to return the WKT.

So if your script can have WKT in input and output that would be perfect for my needs. Is this something not too complicate to do?

On Tue, Aug 5, 2014 at 8:14 AM, Max Demars burton449geo@gmail.com wrote:

Sorry, I was not clear enough because I thought you were working with geodjango already.

The queryset is a list of feature objects where the geometry are stored in a field (in my case the name of that field is geom_multipolygon). The geometry is stored as a GEOS geometry object (django.contrib.gis.geos https://github.com/django/django/blob/master/django/contrib/gis/geos/geometry.py), The GEOS object has many function to retrieve coordinates lie wkt or to transform it to OGRGeometry object.

I can use the GDALSimplifier() with OGR Geometry as input, but probably it would be more straightfoward for django application to allow GEOS Geometry as input and output

I hope that answer your questions.

Thank you for your time and support

-Max Demars

On Mon, Aug 4, 2014 at 4:04 PM, Elliot Hallmark notifications@github.com wrote:

I don't have a geodjango database, so no queryset to look at. But, I'm missing something here. A queryset is giving you a multipolygonfield as a result? Shouldn't it give you the value that field contains? Like if I had an IntegerField, the queryset would tell me the value of the integer there, not give me an IntegerField object to work with. The simplifier can't simplify something without having the points to work with.

How do you get the value of the MultiPolygonField when it is a result in a queryset? What is the type of the value you are getting? Or, how do you get the points to do something with them?

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-51109993 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

DemersM commented 9 years ago

Shapely could be a good solution to translate wkt to numpy array and vice versa...

http://toblerity.org/shapely/manual.html#numpy-and-python-arrays

On Tue, Aug 5, 2014 at 10:04 AM, Max Demars burton449geo@gmail.com wrote:

Im thinking about it... It would be easier to skip the GEOS geometry creation by doing a raw sql query on the db to return the WKT.

So if your script can have WKT in input and output that would be perfect for my needs. Is this something not too complicate to do?

On Tue, Aug 5, 2014 at 8:14 AM, Max Demars burton449geo@gmail.com wrote:

Sorry, I was not clear enough because I thought you were working with geodjango already.

The queryset is a list of feature objects where the geometry are stored in a field (in my case the name of that field is geom_multipolygon). The geometry is stored as a GEOS geometry object (django.contrib.gis.geos https://github.com/django/django/blob/master/django/contrib/gis/geos/geometry.py), The GEOS object has many function to retrieve coordinates lie wkt or to transform it to OGRGeometry object.

I can use the GDALSimplifier() with OGR Geometry as input, but probably it would be more straightfoward for django application to allow GEOS Geometry as input and output

I hope that answer your questions.

Thank you for your time and support

-Max Demars

On Mon, Aug 4, 2014 at 4:04 PM, Elliot Hallmark <notifications@github.com

wrote:

I don't have a geodjango database, so no queryset to look at. But, I'm missing something here. A queryset is giving you a multipolygonfield as a result? Shouldn't it give you the value that field contains? Like if I had an IntegerField, the queryset would tell me the value of the integer there, not give me an IntegerField object to work with. The simplifier can't simplify something without having the points to work with.

How do you get the value of the MultiPolygonField when it is a result in a queryset? What is the type of the value you are getting? Or, how do you get the points to do something with them?

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-51109993 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Permafacture commented 9 years ago

I am familiar with django and am working on a gis project but everything on that so far has been without a geospatial database. Lots of data processing though.

Anyways, I don't know how you got all the way to that index error (GEOSGeometry not having a geom_name method clearly comes up first) but I have added GEOSGeometry. It is only slightly different from OGRGeometry, so it wasn't tough. It is slower though, because the tuple method is much slower.

For WKT, all that would have to be changed is to add a short section to GDALSimplifier.init

WKT is already generated from the numpy array, so self.Geometry = lambda w: w
would be sufficient for wkt->geometry object conversion.

Converting WKT to numpy is the only trick. It is possible to do this poorly so that it is very slow (my first attempt at array to wkt was super slow). Shapely does this by converting wkt to a geometry object as an intermediary. I don't want to add any dependencies either.

Permafacture commented 9 years ago

If you want to look into wkt to numpy, I'd look at parsing modifying the string just a little and then using ast to make a tuple or list. Something like wkt.replace('(((','[[[[').replace(')))',']]]]').replace(',','],[').etc...

Permafacture commented 9 years ago

Also, there was an error because GEOSGeometry doesn't like polygons of only two points. I fixed this in from_number and from_threshold by forcing the number to be at least three.

DemersM commented 9 years ago

Hi Elliot,

Im actually trying to implement your VWSimplifier into a FME transformer. Im trying to send WKT of linestring geometry into the WKTSimplifer but I got the following error. Do you have an idea what is wrong?

wkt = "LINESTRING (-69.66260388 47.72742705,-69.66261612 47.72749302,-69.66248112 47.72760003, [....])" simplifier = WKTSimplifier(wkt)

Traceback (most recent call last): File "", line 1, in File "C:\Python27\lib\site-packages\polysimplify.py", line 184, in init VWSimplifier.init(self,args,*kwargs) File "C:\Python27\lib\site-packages\polysimplify.py", line 68, in init self.thresholds = self.build_thresholds() File "C:\Python27\lib\site-packages\polysimplify.py", line 78, in build_thresh olds nmax = len(pts) TypeError: len() of unsized object

On Tue, Aug 5, 2014 at 5:39 PM, Elliot Hallmark notifications@github.com wrote:

Also, there was an error because GEOSGeometry doesn't like polygons of only two points. I fixed this in from_number and from_threshold by forcing the number to be at least three.

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-51264101 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Permafacture commented 9 years ago

Sorry, the classes are poorly named. Wktsimplifier is just code for converting points to wkt, which gdalsimplifier uses. It shouldn't be it's own class and I'll probably put it in one of the others.

If you want to use wkt directly, without making a geometry object, you'll need to modify the GDALSimplifier as described above.

Permafacture commented 9 years ago

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re
from ast import literal_eval
q = re.compile( '([ 0123456789.]+) ([0123456789.]+)')
def str2tuple(q):
      #for taking the coords seperated by spaces and making a tuple
      return '(%s,%s)' % (q.group(1),q.group(2))
name,pts = wkt.split(' ',1)  #pts is a string
result = literal_eval(q.sub(str2tuple,pts))
DemersM commented 9 years ago

Thank you Elliot!

I suppose you are using Python3. In python 2.7 there is a problem with literal_eval()

see: http://stackoverflow.com/a/20748308

Tell me if you have a solution, I'll try to find something on my side.

-Max Demars

On Wed, Aug 13, 2014 at 7:20 PM, Elliot Hallmark notifications@github.com wrote:

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re from ast import literal_eval q = re.compile( '([ 0123456789.]+) ([0123456789.]+)') def str2tuple(q):

for taking the coords seperated by spaces and making a tuple

  return '(%s,%s)' % (q.group(1),q.group(2))

name,pts = wkt.split(' ',1) #pts is a string result = literal_eval(q.sub(str2tuple,pts))

  • FYI, the VW algorithm is really really fast if you are saving the ordered list of thresholds and reusing it. Otherwise, it is probably slower than C based implementations (like in GEOSGeometry)

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52125063 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

DemersM commented 9 years ago

Furthermore, is this is still correct?

from polysimplify import VWSimplifier import numpy as np from time import time

n = 5000 thetas = np.linspace(0,2*np.pi,n) pts = np.array([[np.sin(x),np.cos(x)] for x in thetas])

start=time() simplifier = VWSimplifier(pts) VWpts = simplifier.from_number(n/100)

I got the following error: : 'numpy.ndarray' object has no attribute 'getAttribute'

On Thu, Aug 14, 2014 at 8:38 AM, Max Demars burton449geo@gmail.com wrote:

Thank you Elliot!

I suppose you are using Python3. In python 2.7 there is a problem with literal_eval()

see: http://stackoverflow.com/a/20748308

Tell me if you have a solution, I'll try to find something on my side.

-Max Demars

On Wed, Aug 13, 2014 at 7:20 PM, Elliot Hallmark <notifications@github.com

wrote:

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re from ast import literal_eval q = re.compile( '([ 0123456789.]+) ([0123456789.]+)') def str2tuple(q):

for taking the coords seperated by spaces and making a tuple

  return '(%s,%s)' % (q.group(1),q.group(2))

name,pts = wkt.split(' ',1) #pts is a string result = literal_eval(q.sub(str2tuple,pts))

  • FYI, the VW algorithm is really really fast if you are saving the ordered list of thresholds and reusing it. Otherwise, it is probably slower than C based implementations (like in GEOSGeometry)

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52125063 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

DemersM commented 9 years ago

I am not enable to use GDALSimplifier neither:

from osgeo import ogr from polysimplify import GDALSimplifier

geom = ogr.CreateGeometryFromWkt(wktGeometry) type(geom)

<class 'osgeo.ogr.Geometry'>

simplifier = GDALSimplifier(geom)

: 'Geometry' object has no attribute 'tuple'

On Thu, Aug 14, 2014 at 10:03 AM, Max Demars burton449geo@gmail.com wrote:

Furthermore, is this is still correct?

from polysimplify import VWSimplifier import numpy as np from time import time

n = 5000 thetas = np.linspace(0,2*np.pi,n) pts = np.array([[np.sin(x),np.cos(x)] for x in thetas])

start=time() simplifier = VWSimplifier(pts) VWpts = simplifier.from_number(n/100)

I got the following error: : 'numpy.ndarray' object has no attribute 'getAttribute'

On Thu, Aug 14, 2014 at 8:38 AM, Max Demars burton449geo@gmail.com wrote:

Thank you Elliot!

I suppose you are using Python3. In python 2.7 there is a problem with literal_eval()

see: http://stackoverflow.com/a/20748308

Tell me if you have a solution, I'll try to find something on my side.

-Max Demars

On Wed, Aug 13, 2014 at 7:20 PM, Elliot Hallmark < notifications@github.com> wrote:

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re from ast import literal_eval q = re.compile( '([ 0123456789.]+) ([0123456789.]+)') def str2tuple(q):

for taking the coords seperated by spaces and making a tuple

  return '(%s,%s)' % (q.group(1),q.group(2))

name,pts = wkt.split(' ',1) #pts is a string result = literal_eval(q.sub(str2tuple,pts))

  • FYI, the VW algorithm is really really fast if you are saving the ordered list of thresholds and reusing it. Otherwise, it is probably slower than C based implementations (like in GEOSGeometry)

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52125063 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Permafacture commented 9 years ago

That is unrelated. I am using literal_eval as intended in python 2.7 On Aug 14, 2014 7:38 AM, "DemarsM" notifications@github.com wrote:

Thank you Elliot!

I suppose you are using Python3. In python 2.7 there is a problem with literal_eval()

see: http://stackoverflow.com/a/20748308

Tell me if you have a solution, I'll try to find something on my side.

-Max Demars

On Wed, Aug 13, 2014 at 7:20 PM, Elliot Hallmark notifications@github.com

wrote:

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re from ast import literal_eval q = re.compile( '([ 0123456789.]+) ([0123456789.]+)') def str2tuple(q):

for taking the coords seperated by spaces and making a tuple

return '(%s,%s)' % (q.group(1),q.group(2)) name,pts = wkt.split(' ',1) #pts is a string result = literal_eval(q.sub(str2tuple,pts))

  • FYI, the VW algorithm is really really fast if you are saving the ordered list of thresholds and reusing it. Otherwise, it is probably slower than C based implementations (like in GEOSGeometry)

— Reply to this email directly or view it on GitHub < https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52125063>

.

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52177615 .

Permafacture commented 9 years ago

from osgeo import ogr

Well, you are using yet another type of input that the code hasn't been written for (which is Geodjago's OGR and GEOS objects.) The error is very clear.

Permafacture commented 9 years ago

Please start new bugs for new issues. I can't tell what the problem is. The test code at the bottom of polysimplify.py, which looks essentially the same as yours, works for me. Does running polysimplify.py work for you? If not, would you post the full trace so I can see better? On Aug 14, 2014 9:03 AM, "DemarsM" notifications@github.com wrote:

Furthermore, is this is still correct?

from polysimplify import VWSimplifier import numpy as np from time import time

n = 5000 thetas = np.linspace(0,2*np.pi,n) pts = np.array([[np.sin(x),np.cos(x)] for x in thetas])

start=time() simplifier = VWSimplifier(pts) VWpts = simplifier.from_number(n/100)

I got the following error: : 'numpy.ndarray' object has no attribute 'getAttribute'

On Thu, Aug 14, 2014 at 8:38 AM, Max Demars burton449geo@gmail.com wrote:

Thank you Elliot!

I suppose you are using Python3. In python 2.7 there is a problem with literal_eval()

see: http://stackoverflow.com/a/20748308

Tell me if you have a solution, I'll try to find something on my side.

-Max Demars

On Wed, Aug 13, 2014 at 7:20 PM, Elliot Hallmark < notifications@github.com

wrote:

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re from ast import literal_eval q = re.compile( '([ 0123456789.]+) ([0123456789.]+)') def str2tuple(q):

for taking the coords seperated by spaces and making a tuple

return '(%s,%s)' % (q.group(1),q.group(2)) name,pts = wkt.split(' ',1) #pts is a string result = literal_eval(q.sub(str2tuple,pts))

  • FYI, the VW algorithm is really really fast if you are saving the ordered list of thresholds and reusing it. Otherwise, it is probably slower than C based implementations (like in GEOSGeometry)

— Reply to this email directly or view it on GitHub < https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52125063>

.

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52186852 .

DemersM commented 9 years ago

Hi Elliot,

Thanks for all your patience. I finally make it works correctly. I am using the VWSimplifier to simplify linestring. It works really well. I am using the from_number() and it works fast enought for my needs.

I manage to create a FME process that can simplify many multipolygon layers at once in a maneer that keeps the topology between layers.

Also, I have noticed that long linestrings with many vertex are much more simplified than small linestring with few vertex. Is this something normal for this algorithm? To keep the equality between long and small linestring, I am chopping them to a maximum number of vertex and recreating the complete linestring after the simplifier.

Thanks again for your great job!

-Max Demars

On Thu, Aug 14, 2014 at 11:55 AM, Elliot Hallmark notifications@github.com wrote:

Please start new bugs for new issues. I can't tell what the problem is. The test code at the bottom of polysimplify.py, which looks essentially the same as yours, works for me. Does running polysimplify.py work for you? If not, would you post the full trace so I can see better? On Aug 14, 2014 9:03 AM, "DemarsM" notifications@github.com wrote:

Furthermore, is this is still correct?

from polysimplify import VWSimplifier import numpy as np from time import time

n = 5000 thetas = np.linspace(0,2*np.pi,n) pts = np.array([[np.sin(x),np.cos(x)] for x in thetas])

start=time() simplifier = VWSimplifier(pts) VWpts = simplifier.from_number(n/100)

I got the following error: : 'numpy.ndarray' object has no attribute 'getAttribute'

On Thu, Aug 14, 2014 at 8:38 AM, Max Demars burton449geo@gmail.com wrote:

Thank you Elliot!

I suppose you are using Python3. In python 2.7 there is a problem with literal_eval()

see: http://stackoverflow.com/a/20748308

Tell me if you have a solution, I'll try to find something on my side.

-Max Demars

On Wed, Aug 13, 2014 at 7:20 PM, Elliot Hallmark < notifications@github.com

wrote:

Turns out regular expressions are pretty fast. I have this going pretty fast but the literal_eval takes 4x longer than all of the string manipulation (q.sub), making it no faster than just converting to OGRGeometry (not GEOS). I'll look into it more, but not necessarily soon. You could always make OGRGeometry as an intermediary also, though I don;t think that helps you.

Here is what I did. result can be a tuple, list or array. Tuple seems the most straight forward

import re from ast import literal_eval q = re.compile( '([ 0123456789.]+) ([0123456789.]+)') def str2tuple(q):

for taking the coords seperated by spaces and making a tuple

return '(%s,%s)' % (q.group(1),q.group(2)) name,pts = wkt.split(' ',1) #pts is a string result = literal_eval(q.sub(str2tuple,pts))

  • FYI, the VW algorithm is really really fast if you are saving the ordered list of thresholds and reusing it. Otherwise, it is probably slower than C based implementations (like in GEOSGeometry)

— Reply to this email directly or view it on GitHub <

https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52125063>

.

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

— Reply to this email directly or view it on GitHub < https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52186852>

.

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52202437 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Permafacture commented 9 years ago

Short linestrings probably have verticies that are all more significant than in long ones, which might be more redundant. Use from_threshold to make sure you are filtering them all to the same level of detail, as opposed to fraction of number of points which is arbitrary.

DemersM commented 9 years ago

That makes a lot of sense! Is this the correct way to use the thresholds?

pts = np.array(feature.getAllCoordinates())
simplifier = VWSimplifier(pts)
thresholds = simplifier.build_thresholds()
VWpts = simplifier.from_threshold(thresholds)

On Fri, Aug 15, 2014 at 9:26 AM, Elliot Hallmark notifications@github.com wrote:

Short linestrings probably have verticies that are all more significant than in long ones, which might be more redundant. Use from_threshold to make sure you are filtering them all to the same level of detail, as opposed to fraction of number of points which is arbitrary.

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52304965 .

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

Permafacture commented 9 years ago

No need to call build_thresholds. It is called by Init. On Aug 15, 2014 8:48 AM, "DemarsM" notifications@github.com wrote:

That makes a lot of sense! Is this the correct way to use the thresholds?

pts = np.array(feature.getAllCoordinates()) simplifier = VWSimplifier(pts) thresholds = simplifier.build_thresholds() VWpts = simplifier.from_threshold(thresholds)

On Fri, Aug 15, 2014 at 9:26 AM, Elliot Hallmark notifications@github.com

wrote:

Short linestrings probably have verticies that are all more significant than in long ones, which might be more redundant. Use from_threshold to make sure you are filtering them all to the same level of detail, as opposed to fraction of number of points which is arbitrary.

— Reply to this email directly or view it on GitHub < https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52304965>

.

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52307229 .

Permafacture commented 9 years ago

Also, that assumes feature is a simple polyline. For polygon and multipolygon, you'll need to use django's OGR or GEOM objects. On Aug 15, 2014 10:53 AM, "Elliot Hallmark" Permafacture@gmail.com wrote:

No need to call build_thresholds. It is called by Init. On Aug 15, 2014 8:48 AM, "DemarsM" notifications@github.com wrote:

That makes a lot of sense! Is this the correct way to use the thresholds?

pts = np.array(feature.getAllCoordinates()) simplifier = VWSimplifier(pts) thresholds = simplifier.build_thresholds() VWpts = simplifier.from_threshold(thresholds)

On Fri, Aug 15, 2014 at 9:26 AM, Elliot Hallmark < notifications@github.com> wrote:

Short linestrings probably have verticies that are all more significant than in long ones, which might be more redundant. Use from_threshold to make sure you are filtering them all to the same level of detail, as opposed to fraction of number of points which is arbitrary.

— Reply to this email directly or view it on GitHub < https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52304965>

.

Stack Overflow: http://stackoverflow.com/users/1914034/burton449 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449 LastFm: http://www.lastfm.fr/user/burton449

— Reply to this email directly or view it on GitHub https://github.com/Permafacture/Py-Visvalingam-Whyatt/issues/5#issuecomment-52307229 .