53

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


28

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


24

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


21

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


19

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


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


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


14

You need to separate them via semi-colons: field1;field2;field3


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


13

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


12

I just implemented this myself and posted my answer over on StackOverflow, but I figured I'd drop my version here for others to view: import numpy as np from scipy.spatial import ConvexHull def minimum_bounding_rectangle(points): """ Find the smallest bounding rectangle for a set of points. Returns a set of points representing the ...


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


12

you need add the layer to project. For add without showing it use: QgsProject.instance().addMapLayer(layer1, False) Example using only one layer: layer1 = QgsVectorLayer(r"C:\test\grassland.shp", 'layer 1', 'ogr') QgsProject.instance().addMapLayer(layer1, False) parameters = {'INPUT': QgsProcessingFeatureSourceDefinition(layer1.id(), True), ...


11

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


10

Yes - processing.run accepts a "feedback" argument, which must be an instance of a QgsProcessingFeedback subclass. If you construct your own feedback object to pass to this function, you can connect to the progressChanged signal and handle progress reports: def progress_changed(progress): print(progress) f = QgsProcessingFeedback() f.progressChanged....


9

The formula for calculating grid convergence (sometimes called meridian convergence) for spherical UTM projections was given (very incorrectly until just now) at How to Calculate North? In case that is not clear γ = arctan [tan (λ - λ0) × sin φ] where γ is grid convergence, λ0 is longitude of UTM zone's central ...


8

To add to previous answers for this post, at least as of QGIS 2.6 does have concave hull algorithm Parameters Input point layer [vector: point] put parameter description here Threshold (0-1, where 1 is equivalent with Convex Hull) [number] put parameter description here Default: 0.3 Allow holes [boolean] put ...


8

You can have a look at the following python versions of those algorithms: Visvalingam–Whyatt Reumann–Witkam Opheim simplification Lang simplification Couldn't find an example for Zhao-Saalfeld (yet).


8

According to the Wikipedia page Longest path problem, this problem ... is NP-hard, meaning that it cannot be solved in polynomial time for arbitrary graphs unless P = NP. Stronger hardness results are also known showing that it is difficult to approximate. However, it has a linear time solution for directed acyclic graphs, which has important ...


8

I think you may be looking for the Graph Diameter of your network. There are a couple of questions on stackexchange that mention this topic, e.g.: The time complexity of finding the diameter of a graph Good algorithm for finding the diameter of a (sparse) graph? What is meant by diameter of a network? Difference between diameter of a graph vs longest path ...


8

The export / save layer as dialog provides a list of fields where users can chose which fields should be included in the export. To reproduce this in a Processing model, you can use the "Drop Field(s)" tool:


7

There is a ready-made tool called ET Spatial's "GeoTools". Where there is a tool called "ET Miscellaneous". In this tool there are two types of way to split polygons viz. By Percentage and By Area ( several units e.g sq m, hectares etc.). This tool splits polygons from any of the four side i.e NESW as the pictures shows. I used "up to down" i.e North to ...


7

Yes, there are such algorithms (see for example here, there and here) but since you seem to have only one single road, would it not be easier to do it by hand? (!). Using for example QGIS, you could import your GPS traces, create a new layer, digitalise the centerline of the bundle, and then export it in whatever format.


7

Looking over your code, I don't see how it will be faster than pure Python. It's exclusively calling Python methods (rasterio methods, np.dstack is Python) and those aren't executed any faster just because the function is compiled with Cython. The key to speeding things up is in here: # Iterate over rows and columns in this block for x in xrange(...


7

There is a paper called "Curved Reconstruction from Unorganized Points" by In-Kwon Lee which looks into constructing lines/curves from a set of points without any ordering by exploiting the moving least-squares method. Although it focues on 2D applications, it mentions the possibility of extending this to higher dimensions. The following image is taken from ...


7

Using a spherical model of the earth may give adequate accuracy and leads to simple fast calculations. Convert all coordinates into earth-centered (3D) cartesian coordinates. For example, the formula (cos(lon)*cos(lat), sin(lon)*cos(lat), sin(lat)) will do. (It uses a distance measure in which the earth's radius is one unit, which is convenient.) ...


7

It is an error. Your link and presumably image source is the 10.0 help. The error is still present in the 10.1 help. However in the 10.2 help, it has been corrected to show only a 1. It's also worth noting that if you work through the math matrix as whuber has done at Arcmap 10 restrict Flow Accumulation, that error is apparent as well as the adjacent 35 ...


7

This is not the answer, I just thought I post Python solution for those who interested: # --------------------------------------------------------------------------- # PAGE MAKER # # --------------------------------------------------------------------------- # Import arcpy module import arcpy, traceback, os, sys from arcpy import env width=650 height=500 ...


7

You can check the source code for the Nearest Neighbour Analysis tool from GitHub. More specifically, the following lines of code which shows how the different parameters are calculated: do = float(sumDist) / count de = float(0.5 / math.sqrt(count / A)) d = float(do / de) SE = float(0.26136 / math.sqrt(( count ** 2) / A)) zscore = float((do - de) / SE) ...


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