36
Here is an example using a SpatialGrid object:
### read shapefile
library("rgdal")
shp <- readOGR("nybb_13a", "nybb")
proj4string(shp) # units us-ft
# [1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333
# +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"
### define ...
27
over() from package sp can be a little confusing but works well. I'm assuming you've already made "A" spatial with coordinates(A) <- ~longitude+latitude:
# Overlay points and extract just the code column:
a.data <- over(A, B[,"code"])
Instead of a point spatial object, this simply gives you a data frame, with the same no. rows as A, and a single ...
26
ST_DWithin was faster in my test than ST_Intersects. That is surprising, especially since the prepared geometry algorithm is supposed to kick in on cases like this. I think there is a chance that this will be quite a lot faster than I showed here.
I did some more tests and two things almost 10-doubled the speed.
First, I tried on a newer computer, but ...
24
Shapefiles have no type MultiPolygon (type = Polygon), but they support them anyway (all rings are stored in one feature = list of polygons, look at Converting huge multipolygon to polygons)
The problem
If I open a MultiPolygon shapefile, the geometry is 'Polygon'
multipolys = fiona.open("multipol.shp")
multipolys.schema
{'geometry': 'Polygon', '...
20
SELECT grid.gid, count(kioskdhd3.geom) AS totale
FROM grid LEFT JOIN kioskdhd3
ON st_contains(grid.geom,kioskdhd3.geom)
GROUP BY grid.gid;
17
This should do what you need:
A select query:
SELECT polygons.id, Count(*)
FROM points
JOIN polygons
ON polygons.ogr_geometry.STContains(points.ogr_geometry) = 1
GROUP BY polygons.id
With an update:
UPDATE polygons
SET [countcolumn] = counts.pointcount
FROM polygons
JOIN
(
SELECT polygons.id, Count(*)
FROM points
JOIN polygons
ON polygons....
16
If "fastest" includes the amount of your time that is spent, the solution will depend on what software you are comfortable with and can use expeditiously. The following remarks consequently focus on ideas for achieving the fastest possible computing times.
If you use a canned program, almost surely the best you can do is pre-process the polygons to set up ...
15
In ArcGIS you can use Hawth's Tools to generate your points or in QGIS you can use Vector->Research Tools->Random Points or Regular Points.
However, you say the points will represent windmill bases and you want to maximise production. So, personally, I would not use either of these methods. Maximising windfarm production is dependent on much more than ...
14
Short answer: No. With this type of UPDATE query, we are updating each row in locations ("Seq Scan"), and the GiST index on the_geom in regions is sufficient in helping limit rows for the ST_Within condition to pair-up the right row from regions.
Longer answer: The magic to figuring this out is to compare what you get from explain query. From pgAdmin III, ...
13
The tool you're looking for is now called Count points in polygons, and it can be found in Processing Toolbox -> QGIS Geoalgorithms -> Vector analysis tools.
12
As with almost all such questions, the optimal approach depends on the "use cases" and how the features are represented. The use cases are typically distinguished by (a) whether there are many or few objects in each layer and (b) whether either (or both) layers allow for precomputing some data structures; that is, whether one or both of them is sufficiently ...
12
Yet another option, this is more of a theory and programmatic one, using arcpy.
A polygon can consist not only of a single outer ring with a single inner donut hole -- they can be nested to an arbitrary number of levels.
Consider the following:
Difference between outer and inner rings http://edndoc.esri.com/arcobjects/8.3/componenthelp/esricore/.%...
12
$sql = "SELECT points.name FROM polygons, points WHERE ST_CONTAINS(polygons.geom, Point(points.longitude, points.latitude)) AND polygons.name = 'California'";
12
Convex hull - as mentioned by Kazuhito - is one option, but - depending on the cluster shape - you will get more appropriate polygons using concave hulls, for example implemented in ConcaveHull plugin.
11
If we examine your polygon:
polygon = shapefile_record['geometry']
print polygon.bounds
(77.84476181915733, 30.711096140487314, 78.59476181915738, 31.28199614048725)
From Shapely manual, object.bounds:
Returns a (minx, miny, maxx, maxy) tuple (float values) that bounds the object.
Here minx = 77.84476181915733, miny = 30.711096140487314
= here, min ...
11
In QGIS, many of the really good tools are in the processing toolbox; you need 'concave hull':
Try it with different threshold values for different levels of detail:
Finally, add a 10% buffer around the outside to make it resemble the sketch you provided:
10
Try Concave hull, What are Definition, Algorithms and Practical Solutions for Concave Hull?
Concave hull has a smaller area, and most of implementations allows you to tune how small and precise resulting polygon should be.
9
Create 2 new columns (x and y) type real. In field calculator use option update existing field and execute following expressions: $x (for x column) and $y (for y column). Actually you can create new columns in field calculator directly using create new column option. Now you can export your layer as CSV (save layer as) and work with it in Excell or whatever ...
9
Another option, without needing the function
update points set country = t1.country from
(
select points.oid, countries.name as country from
countries INNER JOIN points on st_contains(countries.wkb_geometry,points.wkb_geometry)
) t1
where t1.oid = points.oid
I suspect (although I haven't tested) that this will be faster than using a nested ...
9
The Google maps API does not already provide a method for checking points in polygons.
After researching a bit I stumbled across the Ray-casting algorithm which will determine if an X-Y coordinate is inside a plotted shape. This will translate to latitude and longitude. The following extends the google.maps.polygon.prototype to use this algorithm. Simply ...
9
The New York dataset provided in the question is no longer available for download. I use the nc dataset from sf package to demonstrate a solution using sf package:
library(sf)
library(ggplot2)
# read nc polygon data and transform to UTM
nc <- st_read(system.file('shape/nc.shp', package = 'sf')) %>%
st_transform(32617)
# random sample of 5 points
...
8
If getting PostGIS set up right now is more than you want to get involved with, you can get by with probably much less effort in the program you have chosen.
You will want to assign to each of your points the name of the lake so you can sum the catch by the lake variable. This is what ArcGIS folks call a spatial join. In qgis parlance, you can do a couple ...
8
You first need to understand the different ways to create geometries and the geometrical relations in GeoDjango (GeoDjango: GEOS API), without a database:
1) create valid geometries:
# with Point, Polygon objects of GeoDjango
from django.contrib.gis.geos import Point, Polygon,
poly = Polygon(((0.0, 0.0), (0.0, 50.0), (50.0, 50.0), (50.0, 0.0), (0.0, 0.0)))...
8
The point.in.poly function in the spatialEco package returns a SpatialPointsDataFrame object of the points that intersect an sp polygon object and optionally adds the polygon attributes.
First lets add the require packages and create some example data.
require(spatialEco)
require(sp)
data(meuse)
coordinates(meuse) = ~x+y
sr1=Polygons(list(Polygon(cbind(c(...
7
If you had the polygon bounding boxes stored in something like a quad tree then you could use that to quickly determine which polygons to check. At the very minimum you could just see if the point is inside each polygon bounding box as opposed to doing a full point in polygon for each polygon. Personally I'd setup a web service which would cache the ...
7
1) Select your polygon of interest - you can do this manually or with one of the selection tools.
2) Select your 'marginal' points using the Select By Attributes tool with the Method drop-down set to 'Add to current selection'.
You will now have a selected polygon in one layer and selected points in another layer.
3) Use the Select By Location tool to ...
7
Here's a way in R:
Make a test raster, 20x30 cells, make 1/10 of the cells set to 1, plot:
> require(raster)
> m = raster(nrow=20, ncol=30)
> m[] = as.numeric(runif(20*30)>.9)
> plot(m)
For an existing raster in a file, for example a geoTIFF, you can just do:
> m = raster("mydata.tif")
Now get a matrix of the xy coordinates of the 1 ...
7
Spatial Join your points to your polygons, use INTERSECT or WITHIN, no need to keep all the attributes just the OID of the polygon is what's needed on the joined points.
Using summary statistics you can count the points.. use a summary field of FID or OBJECTID depending on what sort of data you have (shape or GDB), summary type of count and a case field of ...
7
You may be interested in Convex Hull which is in Processing | QGIS geoalgorithms | Vector geometry tools. There is Field option which can be used with Method Create convex hulls based on field.
Or from the menu Vector | Geoprocessing Tools | Convex Hull(s).
Many Thanks, Techie_Gus and underdark for information.
7
If you want a really simple algorithm how about this:
Take your point and draw a straight line to the bounding box of your polygon.
Count how many times it crosses the polygon boundary.
If number is odd it must be inside, if even it must be outside.
Only top voted, non community-wiki answers of a minimum length are eligible
