Tag Info

Hot answers tagged

11

If I understand correctly, You want to take (A is the left geometry, B is the right): And extract: A∖AB AB and B∖AB That is - three different geometries for every intersecting pair. First, let's create a view of all intersecting geometries. Assuming your table name is polygons_table, we will use: CREATE OR REPLACE VIEW p_intersections AS -- ...


11

If you group, you should get only unique points. CREATE TABLE test_points as SELECT ST_Intersection(a.geom, b.geom), Count(Distinct a.gid) FROM roads as a, roads as b WHERE ST_Touches(a.geom, b.geom) AND a.gid != b.gid GROUP BY ST_Intersection(a.geom, b.geom) ;


10

Since you said you get a group of intersecting polygons for each polygon you're interested in, you may want to create what is referred to as a "polygon overlay". This isn't exactly what Adam's solution is doing. To see the difference, take a look at this picture of an ABC intersection: I believe Adam's solution will create an "AB" polygon that covers ...


10

the reason it doesn't work with "not intersects" is that you only compare geometries in pair. and there will be the same problem with disjoint. every housepoint will disjoint some parcels even if it intersects one parcell. underdarks suggestion doesn.t have that problem. there is also another trick to use that probably will make the use of indexes more ...


9

GIS (as a field) hasn't done too hot when it comes to really grappling with the surface of the globe. For example, your problem isn't fully defined. Unlike in 2D, where we know the edges of a polygon are composed of straight lines, what are they on the globe? Arcs of great circles, minimizing distance between vertices, are a good choice but not the only ...


9

For a single feature at a time, you can do this pretty easily interactively using the normal Select By Location dialog, using the following key as a guide to the spatial relationship types for line on line overlays (from Select by Location: graphic examples): Select line using line INTERSECT A, C, D, E, F, G, H, I, J CONTAINS ...


8

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


8

Have you tried with shapely ? They describe spatial relationships and it work on windows The spatial data model is accompanied by a group of natural language relationships between geometric objects – contains, intersects, overlaps, touches, etc. – and a theoretical framework for understanding them using the 3x3 matrix of the mutual intersections ...


7

In case there is no specialized function: CREATE table t_intersect AS SELECT hp.gid, hp.st_address, hp.city, hp.st_num, hp.the_geom FROM public.housepoints as hp WHERE hp.gid NOT IN ( SELECT h.gid FROM public.parcel as p, public.housepoints as h WHERE ST_Intersects(h.the_geom,p.the_geom) ) AS foo ...


7

The intersects method in OpenLayers will only return a boolean value, not the intersection of the features. if(polygon1.intersects(polygon2)){ //polygons intersect Edit: I see you mean using the OpenLayers polygon type rather than the OpenLayers intersect method. The JSTS library has only been out a few days, and the intersects function you ...


6

For a Python solution, you may want to look at Shapely http://gispython.org/shapely/docs/1.2/ and RTree http://pypi.python.org/pypi/Rtree/ Rtree will help you create spatial indexes.


6

I have reproduced your example with shapefiles. You can use Shapely and Fiona to solve your problem. 1) Your problem (with a shapely Point): 2) starting with an arbitrary line (with an adequate length): from shapely.geometry import Point, LineString line = LineString([(point.x,point.y),(final_pt.x,final_pt.y)]) 3) using shapely.affinity.rotate to ...


6

You can use the Haversine formula in conjunction with basic trig to iterate a series of vertices describing your circle. Alternatively, if you have access to a GIS (e.g. ArcGIS, QGIS, PostGIS etc) or a GIS API (e.g. OGR, Shapely, GeoTools) you could simply buffer the point by one mile. E.g. for PostGIS you could use ST_Distance_sphere or ...


6

You can use the GDAL/OGR Python bindings for that. from osgeo import ogr wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))" wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 ...


5

Use the return value of the Intersect method instead of the TopologicalOperator. Try the following instead (I use C#, not VB.NET, so hopefully this works. The casting business is really confusing): Dim topoOp As ITopologicalOperator = TryCast(pTestPoly2, ITopologicalOperator) Dim pOutPointCol As IPointCollection = ...


5

You can use the Union tool to create a single featureclass from the two polygons and then calculate the area of each feature within the polygon.


5

Just like set intersection, geometric intersection (which is essentially an intersection on two sets with spatially defined bounds) is both commutative and associative, e.g. (A ∩ B) ∩ C = A ∩ (B ∩ C) = (C ∩ B) ∩ A = C ∩ (B ∩ A). In practice, software often uses the order of inputs to determine how attributes are associated with the output; but the spatial ...


5

Yup, it looks like that is the behaviour from JTS and GEOS. The problem is that your LINESTRING is invalid. If you have PostGIS 2.0, you can use ST_MakeValid(geometry) to fix the LINESTRING to a POINT. This query verifies your bug, and uses ST_MakeValid as a workaround. WITH data AS (SELECT 'POLYGON((150 280, 99 215, 190 210, 150 280))'::geometry AS poly, ...


5

I would recommend you use ST_HausdorffDistance() to caluclate the similarity of your GPS tracks; using ST_Buffer() or ST_Distance() will cause you issues if your tracks get near to each other at some point, but otherwise are divergent. Hausdorff distance is a bit funny to get your head around, but it's ideal for this sort of query, and after a bit of trial ...


5

A layer is composed of one or several geometries. For the intersection of layers, you must iterate through each layer geometries. With shapely it is easy, example with two shapefiles: from osgeo import ogr from shapely.wkb import loads from shapely.geometry import * # first layer, a polygon shapefile first = Polygon() # open shapefile source1 = ...


5

Yes. They are exposed in the bindings: >>> from osgeo import ogr >>> help(ogr.Layer.Intersection) Help on method Intersection in module osgeo.ogr: Intersection(self, *args, **kwargs) unbound osgeo.ogr.Layer method Intersection(self, Layer method_layer, Layer result_layer, char options = None, GDALProgressFunc callback = None, ...


5

You should use esriSpatialRelEnum.esriSpatialRelRelation instead and set the SpatialRelDescription property to pSFilter.SpatialRelDescription = "T***F****". This will only return the features that have the interiors intersecting.


5

Your example may not be the best since it only contains two groups so I added line(point(6 12) point(8 7)) which intersects with 5. This should be possible with a recursive CTE. This works on the test data: WITH RECURSIVE inter_agg AS ( SELECT r.id,array[r.id]||array_agg(r2.id) as arr FROM roads AS r JOIN roads AS r2 ON r.id<>r2.id WHERE ...


4

I'm actually looking for a general solution that will work for many large polygons. Most of the answers I have read so far are naturally focusing on finding a suitable projection for your features. This may be counter-intuitive, but here you need to consider where a projection fails to work, not where it works well. Some projections fail to work only ...


4

What you're describing is the way that the Union operator works in ArcGIS, but its a little different than either Union or Intersection in the GEOS world. The Shapely manual has examples of how sets work in GEOS. However, the PostGIS wiki does have a good example using linework which should do the trick for you. Alternatively, you could compute three ...


4

While an algebraic solution is possible if one assumes the earth is a sphere, we can still handle an ellipsoidal earth using Newton's method and a Esri's projection engine. The projection engine is a c style dll (pe.dll) and is bundled with the freely downloadable ArcGIS Explorer. I think the question could be rephrased as ... A plane is flying from ...


4

Have a look at the answers to this question - Calculating Lat/Lng X miles from point. Using one of the formulas for calculating location from bearing and distance, with a fixed distance and varying bearings, will give you the points of your small circle.


4

What you want is a pivot table. You can create that with the pivot functionality in PostgreSQL, but I think it is easier to create it like this in your case: SELECT m.id, (length_type='00')::int * sum(ST_Length(r.geom)) length_type_00, (length_type='01')::int * sum(ST_Length(r.geom)) length_type_01, (length_type='02')::int * ...


4

Do you need to keep all of your intermediate intersects (ie- temp_poly0, temp_poly1)? If not, you will find that sending results to "in_memory" workspace will not only be easier to manage (all the data is cleaned up from memory on exit of the program, or you could just delete it from memory with arcpy.Delete to clear the memory space if you need to), you ...


4

If it's just search that you're looking for you can do a WFS request with a spatial filter of type Contains or Intersects, something like: <wfs:GetFeature service="WFS" version="1.1.0" xmlns:topp="http://www.openplans.org/topp" xmlns:wfs="http://www.opengis.net/wfs" xmlns="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml" ...



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