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
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
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)]) %>%
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.
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 ...
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
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 ...
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,...
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 ...
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.
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/...
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 ...
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, ...
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 ...
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 ...
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.
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 ...
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 ...
areas = 
for line_feature in line_layer.getFeatures():
cands = area_layer.getFeatures(QgsFeatureRequest().setFilterRect(line_feature.geometry().boundingBox()))
for area_feature in cands:
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 ...
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 ...