Shapely is cool and elegant, but why not using still ogr, with its spatial operators (in OGRGeometry class)?
sample code:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
polyshp = driver.Open('/home/pcorti/data/shapefile/multipoly.shp')
polylyr = polyshp.GetLayer(0)
pointshp = driver.Open('/home/pcorti/data/shapefile/point.shp')
pointlyr = pointshp.GetLayer(0)
point = pointlyr.GetNextFeature()
polygon = polylyr.GetNextFeature()
point_geom = point.GetGeometryRef()
polygon_geom = polygon.GetGeometryRef()
print point_geom.Within(polygon_geom)
Note that you need to compile GDAL with GEOS support.
Sorry for reviving an old thread, the same question came up in the context of the geopackage specification so I thought I might as well answer it here too.
The reason you need all these empty geometry types is because the ISO SQL/MM part 3 specification requires them (not explicitly, but indirectly).
As an example ST_Intersection
needs to be able to return a 'there is no intersection' value. The type of the returned value (what ST_GeometryType
will return when called on the value) depends on the argument types.
ST_Intersection(ST_GeomFromText('POINT (1 1)'), ST_GeomFromText('POINT (0 0)')
for instance is specified in section 5.1.22 as returning 'an empty set of type point'. So this needs to return a valid ST_Geometry
value, ST_IsEmpty
should return true
for it and ST_GeometryType
should return ST_Point
.
So what should you get when you execute ST_AsText
and ST_AsBinary
on that returned value? For WKT the only correct answer I think is POINT EMPTY
.
For WKB this is poorly defined. WKB does not provide an equivalent of WKT's EMPTY
and points don't have a natural empty representation. This is why in GeoPackage there's an empty header bit and it's required to use NaN
as the coordinates of the point. The intention is to avoid people accidentally using empty point sets as a regular point value. NaN works as an absorbing element, so it's pretty obvious when you've made a mistake since the outcome of your calculations is very likely to be NaN as well.
Note that in many implementations, like PostGIS for instance, the above example will return GEOMETRYCOLLECTION EMPTY
instead. That's strictly speaking not compliant with the ISO spec since the type of the return value is not correct.
So to answer the question 'why do you have all these different empty types in WKT and WKB'? Because WKT and WKB are, as far as I know, defined in the SQL/MM part 3 specification, that specification defines distinct typed empty sets for each geometry type and as a consequence WKT and WKB need to have a representation for these typed empty sets.
Best Answer
The structure of WKT was originally defined in the standard OpenGIS® Implementation Standard for Geographic information - Simple feature access - Part 1: Common architecture https://portal.ogc.org/files/?artifact_id=13227. The current version of the standard is https://portal.ogc.org/files/?artifact_id=25355 and that is used as reference below.
EDIT
Excerpt from 7.2.1 BNF Introduction (bolding added):
One can also try to interpret the BNF definition. In WKT geometry names point, linestring etc. are "names". "Name" is made of "letters", "letters" are made of "simple Latin Letters", which are made from a set of small case and upper case letters.
Thus
Point
andPOINT
are both OK, and nothing wrong withpoinT
either.