Hot answers tagged

15

Following the advice at the postgis-users maillist it was possible to identify the existence of outdated ligbeos libraries as the source of the problem: $ locate libgeos_c.so /usr/lib/libgeos_c.so /usr/lib/libgeos_c.so.1 /usr/lib/libgeos_c.so.1.8.2 /usr/local/lib/libgeos_c.so /usr/local/lib/libgeos_c.so.1 /usr/local/lib/libgeos_c.so.1.7.1 /usr/local/src/...


14

In this case, it looks like you need to help the dissolve tool along by fixing the topology first. Here's how you can do this in GRASS. When you load the shapefile, this is how GRASS sees it: Each cross marks an area. You can see that there are some crosses on the boundary lines: These are the sliver polygons caused by the slight offset between the ...


11

If you're comfortable with C/C++, GEOS: http://trac.osgeo.org/geos If you're comfortable with C#, NTS: http://code.google.com/p/nettopologysuite/ If you're comfortable with Java, JTS: http://tsusiatsoftware.net/jts/main.html If you're comfortable with Python, shapely: https://github.com/Toblerity/Shapely If you're comfortable with Ruby, ffi-geos: https://...


10

This happens often with ST_Intersection, irrespective of whether you use ST_SnapToGrid (which is more useful for ensuring a certain precision than for fixing geometry errors) and ST_MakeValid. The problem is to do with the fact that when you intersect polygons, often they will meet at a single point or along a line, as well as producing a (Multi)Polygon ...


7

Use GRASS GIS command v.dissolve (in Sextante plugin -> GRASS Commands -> Vector) If I quickly reproduce you example: The result with v.dissolve: Explanation The problem is easy to understand. A shapefile has no topology, if two areas shared a common border that border would be digitized two times and also stored in duplicate. You see this clearly when ...


7

Lazy geometries are not another kind of geometry. This only describes how Geodjango manages to load, instantiate and use geometry objects, and is pretty much like lazy loading in an O/R mapping framework. If you access a whole bunch of geometry data (lets say from a table), geodjango loads them in text based "Well known text" (WKT) format. If you then pick ...


7

As of Shapely version 1.2.14, coordinates are slicable. This looks very similar to GEOSExtractLine, where a subset of the LineString can be extracted. Here are some examples how you can slice coordinates to extract a new line object: from shapely.geometry import LineString, Point # Original LineString used for examples line = LineString([(30, 50), (60, ...


6

Possible error of your geometric operation depends on: overall size of the objects - bigger size increase errors, projection that you use, datum that you use (each datum suits some parts of the Earth more than the others) quality of your data. Generally you don't want to work with unprojected data at all unless there is some specific reasons like finding ...


5

Ok, got this working - and feel like putting the long answer here, as it has a lot of useful GEOS example bits in it. Here we go. Warnings I haven't compiled this - I stripped out a load of project specific stuff and replaced it with a simple Point class which probably needs a copy/assignment operator. But it worked before I did that. I'm not sure this ...


5

Haven't tried, but found a tweet that looks promising https://mobile.twitter.com/JCSanford/status/281540051203141632 Seems like the concept of "buildpack" is a way to get things loaded on Heroku, and someone has made GEOS iavailable as a buildpack. By setting a custom geos-path you can install Shapely via pip and it will find GEOS.


5

If you inverse the coordinates, it does not work (geopy uses (latitude,longitude) in the WGS84 crs) dublin = (53.33306,-6.24889) liverpool = ( 53.41058,-2.97794) print distance(dublin, liverpool).km 217.863019038 print(vincenty(dublin, liverpool).kilometers) 217.863019038 print(great_circle(dublin, liverpool).kilometers) 217.211596704 GEOS (...


4

I had the same problem, and here is how I solved it: The problem lies in PolygonBuilder.cpp within the libgeos library. The line causing the problem has actually been fixed already, but the fix is not in the current ubuntu/debian repository that is installed via apt-get. Details of the change can be found here. I downloaded the libgeos from github and ...


4

No need to clone any geometries. Also, iterate over geometries[1:] to save one Union(). If osgeo.ogr has something like Shapely's cascaded_union use that instead of accumulating unions in a loop.


4

I used the 2to3 tool (see geographika's answer) to convert Shapely-1.2.13 and then manually had to change 3 lines of code to handle Python 3.x's explicit handling of string and byte encoding. I have put this as a separate answer just for clarity and in case anybody else needs to do the same here are the lines I manually changed: \geos.py (line 97) old: v =...


4

It's not a bug, more like a strategic laziness. 1.X versions of PostGIS only supported GEOMETRYCOLLECTION EMPTY, not other forms of empty. For 2.X I (perhaps foolishly) embraced the full variety of nothings. The result is a not-entirely-complete support of varieties of nothing, made slightly ill-formed by the fact that support libraries like GEOS have their ...


4

You can use Select By Location and then switch the selection by right clicking the layer that has the selection and choosing "Selection" and then "Switch Selection". The image shows something similar to how you might select the points that are identical. The final step would be to export the selection that remains after the switch by right-clicking the layer ...


4

Yes, if you want ST_Area() on a bare geometry to return a planar area, you need to use an area-preserving CRS. GEOS is not magic, it just works in whatever units are handed to it, assuming cartesian math, so the expectation in a planar CRS if you want useful areas/lengths etc.


4

It is a GEOS well-known problem (with Shapely, RGeos-RGdal example or PotGIS example and ...). It simply means there is a line crossing another line and no intermediate coordinate records the intersection or by intersecting line segments which are nearly parallel What are the solutions ? 1) rounding the numerical precision of the input geometries (...


4

It is not clear to me whether you are looking specifically for information on loading point clouds, or geometries in general. GEOS may not be the ideal tool for manipulating point clouds, but you can certainly transport a subset as point geometries. My example below uses a polygon, but the general strategy is geometry-agnostic. (I hadn't tried this before,...


3

Christian, you can't have your circle center and point in long,lat space and have R with mile units. Computational geometry algorithms require you to choose one coordinate system. Latitude and longitude is not a good one to work in, by the way unless your objects are very small in extent. After you transform your points to a local coordinate system and scale ...


3

You can use the distance_lte function in geodjango to check if your point is in a circle or not. It is the same as ST_Distance method in PostGIS. The method which I have mentioned is drawing a 5 meter radius of a circle to find points from Zipcode models. Returns models where the distance to the geometry field from the lookup geometry is less than or ...


3

In regards to Shapely and Python 3.2: We're not yet testing or deploying Shapely on Python 3 + Windows. There's probably some non-conforming code in the shapely.geometry.geo module that prevents its import. June, 2011 - http://old.nabble.com/Shapely-package-installation-td31879542.html So you could try getting the Shapely source from GitHub and ...


3

Python osgeo inlcudes a UnionCascaded method that implements the underlying GEOS CascadedUnion method. It is a little strange in that while Union takes two geometries (in form geom1.Union(geom2)), UnionCascaded operates on a geometry collection. According to GEOS documentation, this method is more efficient than performing a union piecewise on individual ...


3

In my experience using Ubuntu with GIS products will not be problematic. But you should try to keep to Long Term Support (LTS) versions of Ubuntu. The reasons are among others: Cutting edge versions of Ubuntu tend to have fewer compiled versions of GIS-software (both desktop and server) LTS releases are designed to be stable platforms that you can stick ...


3

The points being returned to you by Google Maps are in lat/lng, which almost always these days means SRID=4326. So change your geometry SRID in your PostGIS table to that: location = geos.Point(data['lng'], data['lat'], srid=4326) Then, when displaying points on a web map the projection of the map is commonly set to the Spherical Mercator coordinate system,...


3

The only way I could get POSTGis on "Trusty" was to build from Sources, here is what I used .. sudo mkdir postgis-build cd postgis-build ------------------------------------------------------ then put the script below in a shell script file sudo nano build.sh make it executable and run it sudo chmod a+x build.sh -----------------------------------------...


3

The build is broken. After a long debugging, I found out that if PROJ4 isn't compiled before gdal, the resulting library will NOT work correctly with different reference spaces.


3

Why with Geos? I've used pygc (based on Spheroid) for another Python solution (excerpt below) from pygc import great_circle # New point from initial point, distance, and azimuth great_circle(distance=111000, azimuth=65, latitude=30, longitude=-74) {'latitude': 30.41900364921926, 'longitude': -72.952930949727573, 'reverse_azimuth': 245.52686122611451} ...


3

another method I found import pyproj endLon,endLat,backAzimuth = (pyproj.Geod(ellps='WGS84') .fwd(startLon,startLat,forwardAzimuth,distance)) See http://jswhit.github.io/pyproj/pyproj.Geod-class.html#fwd


3

ST_ConvexHull It is often used to determine an affected area based on a set of point observations. Also worth nothing, It is usually used with MULTI and Geometry Collections. Although it is not an aggregate - you can use it in conjunction with ST_Collect to get the convex hull of a set of points. ST_ConvexHull(ST_Collect(somepointfield)). So run ...


Only top voted, non community-wiki answers of a minimum length are eligible