38

You could try 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 of ...


25

This method uses the intersect() function from the raster package. The example data I've used aren't ideal (for one thing they're in unprojected coordinates), but I think it gets the idea across. library(sp) library(raster) library(rgdal) library(rgeos) library(maptools) # Example data from raster package p1 <- shapefile(system.file("external/lux.shp", ...


22

Go to Settings > Options > Processing and under General change to Ignore features with invalid geometries.


21

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) ;


21

Why don't you: Select the two polygons you want to intesect Enable Editing Edit-> Merge Selected Features Save edits


20

Try Intersect tool in the analyses toolbox if you are on ArcGIS. Set output to point.


20

Here's an alternate approach using the new sf package, which is meant to replace sp. Everything is much cleaner, and pipe friendly: library(sf) library(tidyverse) # example data from raster package soil <- st_read(system.file("external/lux.shp", package="raster")) %>% # add in some fake soil type data mutate(soil = LETTERS[c(1:6,1:6)]) %>% ...


18

The main difference will be in the attributes of the results. When using Clip only the input feature’s attributes will be in the output (none from the clip feature), where if you used Intersect the attributes form all features used will be in the output.


16

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


16

The nodes: You want two things, the end points of the polylines (without intermediate nodes) and the intersection points. There are an additional problem, some polylines end points are also intersection points: A solution is to use Python and the modules Shapely and Fiona 1) Read the shapefile: from shapely.geometry import Point, shape import fiona lines ...


15

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


15

While I'm a big user of both shapely and fiona, I wouldn't go this approach. This is a task of writing an effective SQL statement. Using ogr2ogr with an SQLITE dialect, you can process this from a command line. Change directory to one before the shapefiles, so that all of the shapefiles are in one directory called data. OGR treats directories of shapefiles ...


13

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): (source: arcgis.com) Select line using line INTERSECT A, C, D, E, F, G, H,...


13

Look at Martin Davis (creator of the JTS Topology Suite), Lin.ear th.inking: Quirks of the "Contains" Spatial Predicate Geometry A contains Geometry B if no points of B lie in the exterior of A, and at least one point of the interior of B lies in the interior of A Geometry A covers Geometry B if no points of B lie in the exterior of A All that is ...


12

You can do it with the Vector tools built into the newer versions of QGIS. I am using 1.8 at the moment, and it has a Vector menu. Got to: Vector > Geoprocessing Tools > Intersect. Alternatively, you can probably devise query on the attribute table that will select just the features you want. Once selected you can export that to another layer.


12

I encountered similar issues as well with polygons. Maybe you have a similar problem. Error Message by ESRI: "Invalid Topology (Incomplete Void Poly)" Actual Error: "Invalid Geometry" Fix: Run "Repair Geometry" (changes data in-place, be careful, there is no undo) What happens is that the error reported is not using the ESRI terminology of Topology/...


12

There's a couple of ways of going about this but you probably want to dissolve the features (Vector->Geoprocessing Tools->Dissolve). With dissolve you don't need to select anything first as it is all done from the attributes. So, let's say you have a field called 'Type' (for example). Then in your example your polygons would all be type 'A' (and you ...


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

The behaviour of gIntersection is not to pass any intersected data by design: Since there are no general matches between intersected spatial objects, any arbitrary operations on attributes require assumptions about unknown user intentions. This is why no data slots should be passed through ... ... The design of gIntesection() is inentional, ...


10

You're right, using ST_Intersection slows down your query noticeable. Instead of using ST_Intersection it is better to clip (ST_Clip) your raster with the polygons (your fields) and dump the result as polygons (ST_DumpAsPolygons). So every raster cell will be converted into a little polygon rectangle with distinct values. For receiving min, max or mean ...


10

The accepted method did not work for me because my basemap layer wouldn't show up in the Input Features dropdown. I solved this by doing the following: At the View menu, choose Data Frame Options. At the Data Frame Tab look for "Clip Options" Choose Clip to shape Then click the Specify Shape button Then select the boundary layer as input Apply and the ...


9

The simplest approach would be library(raster) x <- spdf1 - spdf2 # or, more formally y <- erase(spdf1, spdf2) See ?'raster-package' (section XIV) for more functions that deal with polygon overlay. These functions use the base-functions of rgeos under the hood, in 'user-level' (as opposed to 'developer-level') functions.


9

Here's a query that should get you what you want: SELECT a.name, b.name, ST_LENGTH(ST_Intersection(a.geom, b.geom)) FROM polygons a, lines b WHERE ST_Intersects(a.geom, b.geom);


9

Super easy manual process. You use the tool Select by Location. Export all points in B to a new layer C Select points in A that match B. "Switch" the selection in A. You now have selected all the points in A that are not in B. Append those selected A points to your C layer. The C layer now contains all points that are in both A & B, uniquely in A, and ...


8

My guess would be that it fails if the intersection returns NULL. So you should add a where clause checking if there actually is an intersection before you try to create the WKT.


8

I would create intersection points, then buffer them by the distance of 'space' you want. Intersect the buffer with the lines, and delete those that intersect. Done.


8

There are some errors in your script but it is not the most important problem: You cannot create a valid shapefile without specifying the geometry of the layer: driver = ogr.GetDriverByName('ESRI Shapefile') dstshp = driver.CreateDataSource('SomeFilename.shp') dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon) And you don't know a priory ...


8

areas = [] for line_feature in line_layer.getFeatures(): cands = area_layer.getFeatures(QgsFeatureRequest().setFilterRect(line_feature.geometry().boundingBox())) for area_feature in cands: if line_feature.geometry().intersects(area_feature.geometry()): areas.append(area_feature.id()) area_layer.select(areas)


8

To get the appropriate cell coordinates from your latitude and longitude you need to know the coverage of the netCDF file (for example, the coordinate of the upper left cell and the size of the cells in the x and y directions). @Luke 's solution below works for a srtaight x/y raster with no cell rotation. If you don't have that, you could look at the affine ...


8

In theory the queries that you have done should return the polygons you said haven't been returned. That makes me suspect that you might be encountering floating point error issues that SQL Server has with it spatial data types. Hence my comment about buffering the bounding polygon with a minimal amount. So something like the following should get the ...


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