A message from our CEO about the future of Stack Overflow and Stack Exchange. Read now.
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 ...


45

Yes, there is an analytical solution for this problem. The algorithm you are looking for is known in polygon generalisation as "smallest surrounding rectangle". The algorithm you describe is fine but in order to solve the problems you have listed, you can use the fact that the orientation of the MAR is the same as the one of one of the edges of the point ...


44

Every polygon has, at a minimum, four distinct "centers": The barycenter of its vertices. The barycenter of its edges. Its barycenter as a polygon. A GIS-specific "center" useful for labeling (usually calculated with undocumented proprietary methods). (They may accidentally coincide in special cases, but for "generic" polygons they are distinct points.) A ...


42

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


40

To supplement @julien's great solution, here is a working implementation in R, which could serve as pseudocode to guide any GIS-specific implementation (or be applied directly in R, of course). Input is an array of point coordinates. Output (the value of mbr) is an array of the vertices of the minimum bounding rectangle (with the first one repeated to ...


31

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


19

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

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

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


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

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


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


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

It's simple but messy. Because you're working in ECEF, presumably you have the ray's origin (x,y,z) and direction vector (u,v,w) in ECEF coordinates, too. For the moment let's assume that during the time of travel to the earth's surface, the earth does not appreciably move. (The fastest part of the rotating earth, the Equator, moves about 0.45 km/sec and ...


12

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


12

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


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

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

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


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


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