88

Different approach. Knowing that the pain is in ST_Intersection, and that true/false tests are fast, trying to minimize the amount of geometry passing through the intersection might speed things up. For example, parcels that are totally contained in a jurisdiction don't need to be clipped, but ST_Intersection will still probably go to the trouble of building ...


47

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


33

UTM uses a transverse Mercator projection with a scale factor of 0.9996 at the central meridian. In the Mercator, the distance scale factor is the secant of the latitude (one source: http://en.wikipedia.org/wiki/Mercator_projection), whence the area scale factor is the square of this scale factor (because it applies in all directions, the Mercator being ...


20

If you were in a plane, then the point that is r meters away at a bearing of a degrees east of north is displaced by r * cos(a) in the north direction and r * sin(a) in the east direction. (These statements more or less define the sine and cosine.) Although you are not in a plane--you're working on the surface of a curved ellipsoid that models the Earth'...


18

To save you some time here is @MerseyViking answer in javascript: function radians(n) { return n * (Math.PI / 180); } function degrees(n) { return n * (180 / Math.PI); } function getBearing(startLat,startLong,endLat,endLong){ startLat = radians(startLat); startLong = radians(startLong); endLat = radians(endLat); endLong = radians(endLong); ...


18

The compactness of an object can be measured using the Polsby-Popper test by determining the Polsby-Popper (PP) score. The PP score is determined by: multiplying the polygon's area by 4pi and dividing by the perimeter squared. Using this, a circle will have a score of 1 and any other geometric shape has a smaller ratio. disc :(4*PI)* PI*R² / 4PI²R²= 1 ...


17

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


17

The "minimum bounding geometry" and "clip polygon" algorithms in QGIS are implemented in /python/plugins/processing/algs/qgis/MinimumBoundingGeometry.py and /src/analysis/processing/qgsalgorithmclip.cpp. If you follow through the source of these, you'll find that they rely on geometry-related functions from a C++ class called QgsGeometry, specifically ...


16

Introduction Because this issue (of discrepancies in standard deviations, variances, or other statistical summaries) comes up periodically, especially when a thoughtful and careful GIS analyst checks their work, I thought it would be good to share the "forensic analysis" of the discrepancy so that readers can carry out similar checks in their own ...


16

This likely requires some scripting in any GIS platform. The most efficient method (asymptotically) is a vertical line sweep: it requires sorting the edges by their minimum y-coordinates and then processing the edges from bottom (minimum y) to top (maximum y), for a O(e * log(e)) algorithm when e edges are involved. The procedure, although simple, is ...


15

When the curve is comprised of line segments, then all interior points of those segments are inflection points, which is not interesting. Instead, the curve should be thought of as being approximated by the vertices of those segments. By splining a piecewise twice-differentiable curve through those segments, we can then compute the curvature. An ...


15

Use a spatial projection library to do the hard work. Adapting from a previous answer, use a dynamic azimuthal equidistant projection to do a geodesic buffer. from functools import partial import pyproj from shapely.ops import transform from shapely.geometry import Point proj_wgs84 = pyproj.Proj('+proj=longlat +datum=WGS84') def geodesic_point_buffer(lat, ...


14

It is a consequence of a theorem of Archimedes (c. 287-212 BCE) that for a spherical model of the earth, the area of a cell spanning longitudes l0 to l1 (l1 > l0) and latitudes f0 to f1 (f1 > f0) equals (sin(f1) - sin(f0)) * (l1 - l0) * R^2 where l0 and l1 are expressed in radians (not degrees or whatever). l1 - l0 is calculated modulo 2*pi (e.g., -179 - ...


14

This is a graph coloring problem. Recall that a graph coloring is an assignment of a color to the vertices of a graph in such a way that no two vertices which share an edge will also have the same color. Specifically, the (abstract) vertices of the graph are the polygons. Two vertices are connected with an (undirected) edge whenever they intersect (as ...


14

I've figured out an algorithm for the grid approach using several Python tools. Rasterising and polygonising is done with rasterio, which is based on GDAL/OGR. Here are most of the imports: import rasterio import numpy as np from rasterio import Affine, features from shapely.geometry import mapping, shape from shapely.ops import cascaded_union from math ...


14

I couldn't stop thinking about this... I was able to come up with a Stored Procedure to do the loop counting. The example path contains 109 loops! Here are the flight points shown with the loop centroids in red: Basically, it runs through the points in the order they were captured and builds a line as it iterates through the points. When the line we are ...


13

It looks to me like you need to perform the trigonometry in radians not degrees. You use a function toDeg() so presumably you have one called toRad() (or possibly fromDeg() if you're odd). Call that function with your latitude and longitude values before the calculations, and you should be set. Edit I just tried this in Python (the syntax isn't dissimilar ...


13

Not sure if this is helpful, let's go down the rabbit hole. First we need to know what's happening when you call the function in SQL. To do this we reference the output of \dS+. \dS+ shows the override table with one entry for every function and prototype, as well as the function that it is dispatching to. You're calling ST_Distance($1::geog,$2::geog). ...


12

If you are interested in an implementation look at jsts a Javascript implementation of the much used Java Topology Suite library -- depending on whether you prefer reading Javascript or Java, I suppose. A general idea of how the algorithm works. For points, it is trivial, you simply buffer them by a given radius. If you have multiple points, you will have ...


11

Contraction Hierarchy is a very fast algorithm: http://algo2.iti.kit.edu/1087.php This algorithm is RAM friendly while executing a query (to hold a contracted graph some more RAM is necessary as well as massive preprocessing) There are some other algorithms - including the ones that solve public transit routing: http://i11www.iti.uni-karlsruhe.de/members/...


11

After reading your Question, and the various Answers, I got interested in this problem. After doing a bit of reading on Map-matching algorithms, I have understood the following: To Match the gps Location to road, you need the actual road data in vector format It will help if you have different weights for different roads. So the chances of a point matching ...


11

Here is an R function that implements the Alpha Hull model. The output is an sp polygon object. Please see the example in the header. It requires the sp, alphahull and maptools packages. **Update (01-15-2018) There have been numerous changes to the resulting objects produced by the alphahull package. As such, I needed to rewrite the function. I added a ...


11

you can check out k-means clustering algorithm here. In data mining, k-means clustering is a method of cluster analysis which aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean. This results into a partitioning of the data space into Voronoi cells. kmeans-postgresql implementation ...


11

What if you used a grid approach? Convert each of the polygons to a raster, assign each cell a value of 1, then add all of the rasters together. The resulting polygon would be formed by the highest-value cells whose area equals the average area of the input polygons. The cell size would have to be small enough to make the difference between the different ...


11

You could have a look at the following method : skeletonize your polygons and rather work on line type features related to your original polygon with a unique source polygon ID. I guess there's some guesses to do (for example, when to consider a polyline as a real centerline : minimal length for a polyline to be eligible to centerline status). When the ...


10

An elegant principle provides a simple answer: All points on a smooth curved surface are flat at a sufficiently large scale. This means that after affine change of coordinates (usually involving just a rescaling of one of them), we can use formulas of Euclidean geometry, such as the Pythagorean Theorem for computing distances and the negative-reciprocal-...


10

This requires a kind of "field calculation" in which the value computed (based on a latitude, longitude, central azimuth, uncertainty, and distance) is the bowtie shape rather than a number. Because such field calculation capabilities were made much more difficult in the transition from ArcView 3.x to ArcGIS 8.x and have never been fully restored, nowadays ...


10

The Flow Direction tool in Spatial Analyst only supports the D8 flow direction model, as stated in the How flow direction works page: There are eight valid output directions relating to the eight adjacent cells into which flow could travel. This approach is commonly referred to as an eight-direction (D8) flow model and follows an approach presented in ...


10

The following relies on the Wikipedia article on seven-parameter Helmert transformations. The data ("double points") consist of ordered pairs ((x,y,z), (x',y',z')) where (x,y,z) are earth-centered Cartesian coordinates in the source datum and (x',y',z') are the corresponding points in the target datum, all measured in meters. The latter are presumed ...


10

If I understand you right you want to cluster lines that is about the same without respect to direction. Here is an idea that I think could work. split the lines in start point and end point Cluster the points and get cluster id Find lines with the same combination of cluster id. Those are a cluster This should be possible in PostGIS (of course :-) ) ...


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