__geo_interface__ for shapes

[Christian Ledermann]

Hello,

It woulld be nice if the shape would implement the __geo_interface__ convention:

https://gist.github.com/2217756 to be able to easily exchange records

from pyshp to shapely, pysal, geojson, ...

I wrote a library http://pypi.python.org/pypi/pygeoif/ that supports this 

protocol so you might want to take a look at it

[Christian Ledermann]

it looks to me as if the implementation would be trivial
(code not yet tested just of the top of my head something
like this should work).

add to the class _Shape :


@property
def __geo_interface__(self):
if self.shapeType in [POINT, POINTM, POINTZ]:
return {
'type': 'Point',
'coordinates': tuple(self.points)
}
elif self.shapeType in [MULTIPOINT, MULTIPOINTM, MULTIPOINTZ):
return {
'type': 'MultiPoint',
#'coordinates': tuple(self.points)
}
elif self.shapeType in [POLYLINE, POLYLINEM, POLYLINEZ]:
return {
'type': 'LineString',
#'coordinates': tuple(self.points)
}
elif self.shapeType in [POLYGON, POLYGONM, POLYGONZ]:
return {
'type': 'Polygon',
#'coordinates': tuple(self.points)
}



I am not acquainted with the shapefile format enough to actually implement this
right now, but if you think this is a valuable addition to pyshp I
fork it and give it
a shot.

Q: does POLYLINE represent a LineString or a MultiLineString?

--
Best Regards,

Christian Ledermann

Nairobi - Kenya
Mobile : +254 702978914

<*)))>{

If you save the living environment, the biodiversity that we have left,
you will also automatically save the physical environment, too. But If
you only save the physical environment, you will ultimately lose both.

1) Don’t drive species to extinction

2) Don’t destroy a habitat that species rely on.

3) Don’t change the climate in ways that will result in the above.

}<(((*>

[Joel Lawhead]

Christian,

I think this is brilliant. Ultimately pyshp will be useful to more people as a pure python shapefile driver in a larger API than as a simple data access library.

POLYLINES can be both. "Poly" and "Multi"-type shapefiles are a collection shapes made of parts. Parts are just a collection of points. Each shape is associated with a dbf attribute record. Parts within a shape may or may not be connected (i.e. mulitlinestring?). In pyshp all of the shape points are stored in the ".points" list. The ".parts" attribute is a list of indexes in the points list where each part starts.

Is a MultiLineSting equivalent to a shapefile shape with disconnected parts?

- Joel

--
Joel Lawhead, PMP
NVision Solutions Inc.
Phone: 228-242-0014
Fax: 228-242-0013
Cell: 228-342-1891
Email: jlaw...@nvs-inc.com
Web: http://www.nvs-inc.com

[Christian Ledermann]

Joel,

attached a first patch, this works for simple shapefiles with Points and
(Multi)Linestrings.

to test it I wrote a simple script that converts shapefiles to kml:
https://github.com/cleder/geo_file_conv


Q:
How do i find out what the exteriors and holes of a multipolygon are?
parts just gives me start and end of the LinearRings but no indication
what these linear rings are.

Does anybody have a shapefile with MultiPoints so i can test this?

I came across some issues:
The encoding of the dbf is not processed, i ran into decode issues
Did you have a look at http://pypi.python.org/pypi/dbf/ this library
handles the
encoding pretty well

some shapefiles could not be processed, I cannot quite put my finger on it
what is wrong, I will send you the shapefiles when i confirmed the errors.

On Thu, Aug 23, 2012 at 5:20 AM, Joel Lawhead
<jlaw...@geospatialpython.com> wrote:
> Christian,
>
> I think this is brilliant. Ultimately pyshp will be useful to more people as a pure python shapefile driver in a larger API than as a simple data access library.
>
> POLYLINES can be both. "Poly" and "Multi"-type shapefiles are a collection shapes made of parts. Parts are just a collection of points. Each shape is associated with a dbf attribute record. Parts within a shape may or may not be connected (i.e. mulitlinestring?). In pyshp all of the shape points are stored in the ".points" list. The ".parts" attribute is a list of indexes in the points list where each part starts.
>
> Is a MultiLineSting equivalent to a shapefile shape with disconnected parts?

yes it is

[Christian Ledermann]

BTW,
have a look at:

https://github.com/Vizzuality/cartodb/blob/develop/lib/importer/misc/shp_normalizer.py

nice way to determine the EPSG from a prj file ;)

On Thu, Aug 23, 2012 at 7:11 PM, Christian Ledermann

[Christian Ledermann]

Hi Joel,

> https://github.com/cleder/geo_file_conv
>
>
> Q:
> How do i find out what the exteriors and holes of a multipolygon are?
> parts just gives me start and end of the LinearRings but no indication
> what these linear rings are.
>

I found something in the blog:

-----------------
The points is just a list containing lists of the point coordinates.
Two things to note: If the geometry has multiple parts, such as
multiple polygons, the points for all parts are just lumped together.
You must separate them by referencing the parts index list. Some
shape types allow for z and m values which may appear in addition to
the x,y pairs.


Joel LawheadNovember 6, 2011 8:07 AM
Yes - point order matters in a polygon. Points should be clockwise and
polygon holes are counterclockwise according to the spec.

----------------------------------------------------


so the pnpoly_orig from
http://geospatialpython.com/2011/01/point-in-polygon.html
should give me the orientation of a LinearRing ?

unluckily the leading whitespace is missing, do you have the correct
formatted version?

[Christian Ledermann]

Ahh OK, it was a long Monday, the above was to determine if a point is inside
a polygon, not if it is clockwise or counterclockwise.

this should do the trick (from http://toblerity.github.com/fiona/manual.html ):

def signed_area(coords):
"""Return the signed area enclosed by a ring using the linear time
algorithm at http://www.cgafaq.info/wiki/Polygon_Area. A value >= 0
indicates a counter-clockwise oriented ring.
"""
xs, ys = map(list, zip(*coords))
xs.append(xs[1])
ys.append(ys[1])
return sum(xs[i]*(ys[i+1]-ys[i-1]) for i in range(1, len(coords)))/2.0


I'll continue with the investigation in a few days

[Christian Ledermann]

OK it took more than a few days, but now it works:

https://github.com/cleder/pyshp

you can test it with

https://github.com/cleder/geo_file_conv

shp_to_kml.py