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
Go to Settings > Options > Processing and under General change to Ignore features with invalid geometries.
Alternatively, it's also worth checking the answer by A.Oikonomidis as well as other tools available in the processing toolbox to fix invalid geometries in the original dataset.
Here's an alternate approach using the new sf package, which is meant to replace sp. Everything is much cleaner, and pipe friendly:
# 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)]) %>%
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.
# Example data from raster package
p1 <- shapefile(system.file("external/lux.shp", ...
If you group, you should get only unique points.
CREATE TABLE test_points as
roads as a,
roads as b
AND a.gid != b.gid
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.
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
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 ...
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 ...
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 ...
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 ...
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 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/...
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 ...
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://...
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, ...
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 ...
areas = 
for line_feature in line_layer.getFeatures():
cands = area_layer.getFeatures(QgsFeatureRequest().setFilterRect(line_feature.geometry().boundingBox()))
for area_feature in cands:
The simplest approach would be
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.
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 ...
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 ...
I ran the Check validity from the dropdown
Used the default settings
The result was an Invalid Output
I then copied and pasted the co-ordinated into the QGIS project window screen (Centre bottom) and zoomed in until I found an intersection...
I deleted three nodes to remove the spike intersecting and ran the Intersection to get the result.
This is a great application for a user-defined aggregate function, and I'm a bit wondering why this particular aggregate doesn't already exist in PostGIS. At its core, an aggregate function needs to do nothing more than iterate over a set of rows, maintaining a state (of type stype), and repeatedly calling a function (sfunc) that transforms a state and a row ...
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 priori (...
If you've got a polygon you want to use as a declared variable and intersect it with a table containing existing geometry, your query (including your polygon variable declaration) would look something like this:
(MSSQL Server syntax)
declare @polygon geometry = 'POLYGON((-9486683.581 4810152.256, -9282073.762 4821688.121, -9262037.786 4625578.413, -9477576....
If your polygons are slivers the eliminate command works well to merge them into either the larger area polygon or the longest edge.
If the polygons are overlaps then there may be an easier way, but I would select out the overlaps to a separate layer, then union them back in, creating the slivers and using the eliminate command.