37
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 ...
32
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 ...
28
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', '...
26
SELECT grid.gid, count(kioskdhd3.geom) AS totale
FROM grid
LEFT JOIN kioskdhd3 ON st_contains(grid.geom,kioskdhd3.geom)
GROUP BY grid.gid;
18
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.
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....
17
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.
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 ...
13
$sql = "SELECT points.name FROM polygons, points WHERE ST_CONTAINS(polygons.geom, Point(points.longitude, points.latitude)) AND polygons.name = 'California'";
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
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
...
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
People using a wgs84 geometry will probably have trouble with this function since
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1))::int,$2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1))::int,$2) as y
only return integers. Except for very big geometries such as countries (that are laying on multiple lat, lng degrees), ...
10
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 ...
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.
10
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.
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 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(...
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)))...
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
Your guess is correct - use 'Count Unique Points in Polygon' in the processing toolbox, using "car_id" as the "Class field".
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
The problem in your first example is in this loop:
...
shpfilePoints = []
for shape in polygon:
shpfilePoints = shape.points
...
It only appends the last feature points. I tried out my approach with this shapefile:
I modified your code to:
from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile
path = '/home/zeito/pyqgis_data/...
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.
7
Use this script:
from PyQt5.QtCore import QVariant
# layers
polygon_lyr = QgsProject.instance().mapLayersByName("Polygons")[0]
point_lyr = QgsProject.instance().mapLayersByName("Points")[0]
# index number of "class" field
field_index = point_lyr.fields().indexFromName('class')
# unique class names
unique_classes = point_lyr.uniqueValues(field_index)
...
6
Found this function on the PostGIS mailinglist. I guess it is what you require:
CREATE OR REPLACE FUNCTION point_inside_geometry(param_geom geometry)
RETURNS geometry AS
$$
DECLARE
var_cent geometry := ST_Centroid(param_geom);
var_result geometry := var_cent;
BEGIN
-- If the centroid is outside the geometry then
-- calculate a box ...
6
First create a polygon grid using the Vector Grid Tool (Vector\Research Tools)
You can specify the polygon dimensions in the settings
Second run a spatial query to intersect your point dataset with the grid cells
Save selection as a new layer
Only top voted, non community-wiki answers of a minimum length are eligible
