New answers tagged

3

Modified version of createBuffer function using a column as the distance value is here: import ogr, os def createBuffer(inputfn, outputBufferfn, bufferField): inputds = ogr.Open(inputfn) inputlyr = inputds.GetLayer() shpdriver = ogr.GetDriverByName('ESRI Shapefile') if os.path.exists(outputBufferfn): shpdriver.DeleteDataSource(...


3

You can use the GeoPandas library to read in your Shapefile as a GeoDataFrame. Then you can make use of the function written in this answer. Here's a quick mock-up of what it would look like. import geopandas as gpd # Copy the definition of the make_variable_buffer function # from here: https://gis.stackexchange.com/a/398549/123828 # Reading in the input ...


0

You can merge geometries from all features in your feature collection into one multiPolygon and than use a .dissolve() function. var filterAll = ee.Filter.inList('ADM2_NAME', ['Misamis Oriental', 'Davao del Norte', 'Bukidnon']); var all = ee.FeatureCollection(ADM2.filter(filterAll)).geometry().dissolve() Link to the code: https://code.earthengine.google....


1

Use the "Join attributes by location (summary)" tool with: 'Input' : name_of_grid_layer 'Join layer' : name_of_polyline_layer 'Geometric predicate' : intersects 'Fields to summarize (leave empty to use all fields)' : x 'Summaries to calculate (leave empty to use all fields)' : sum.


8

To merge all the features of layer 'polygon2' inside the same same polygon from layer 'polygon1' automatically (without manual selection), use Menu Processing / Toolbox / Aggregate (see documentation) with 'polygon2' as input layer and this expression as Group by expression: array_first(overlay_within('polygon1', $id)) 2 purple polygons = layer 'polygon1'; ...


6

If you want to dissolve in-place: Start edit mode for the layer. Select the bounding and interior polygons Click "Merge selected features" button () in the Advanced digitizing Toolbar. Repeat 2-3 for all bounding and interior polygons separately.


6

To get as a different layer from the original one, use the following steps. Select the polygons you want to dissolve Run "Vector Geometry" -> "Dissolve tool". (Check "Selected features only") Then run "Vector geometry" -> "Multipart to singleparts" tool for the dissolved layer.


2

Here is an example of how this can be done. Replace the values of the first four variables to match your use case: # -*- coding: utf8 -*- from qgis.PyQt.QtCore import * from qgis.PyQt.QtGui import * from qgis.PyQt.QtWidgets import * from qgis.utils import iface from qgis.core import * import processing layers = [layer for ...


0

You are telling GDAL to burn all 0 values with burn_values=[0]. Additionally, setting xRes and yRes to 100 won't work as your data uses lon, lat coords (1 degree is 111,000m at the equator so 100 degrees is huuuuge). By the way, the newer gdal.Rasterize is easier to use than the old gdal.RasterizeLayer. Particularly there's no need to pre-create the output ...


1

There are more than one way. I would try Spatial Join with match option HAVE_THEIR_CENTER_IN —The features in the join features will be matched if a target feature's center falls within them. The center of the feature is calculated as follows: for polygon and multipoint the geometry's centroid is used, and for line input the geometry's midpoint is used To ...


2

https://github.com/corteva/geocube from geocube.api.core import make_geocube vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.shp' out_grid = make_geocube( vector_data=vector_fn, measurements=["column_name"], resolution=(-25, 25), fill=-9999, ) out_grid["column_name&...


0

The parameters you have entered into the Create() function are incorrect. You can find them using the python help function: from osgeo import gdal driver = gdal.GetDriverByName('GTiff') help(driver.Create) Create(Driver self, char const * utf8_path, int xsize, int ysize, int bands=1, GDALDataType eType, char ** options=None) In your case, you have entered ...


0

first, select all yellow polygons and do union. second, select all the blue ones and do union.and lastly for the gray ones. Now you should have three features. Then give a unique attribute to each. e.g. unionYellow = 0 , unionBlue = 1, unionGray = 2 then use merge shapefiles tool or any other merger tool. You can use the numbers to change their order. Now ...


1

Another approach is to use the "Extract by expression" geoalgorithm together with the following expression: overlay_intersects(@layer) it uses the overlay_intersects() function. Input: Algorithm's window: Output:


3

You can use "Field Calculator" with an expression like this: array_to_string( array:=overlay_intersects( layer:='counties', expression:=kom_kod ), delimiter:=', ' )


1

If you want to divide polygon before it's published to the web, you can use GIS tools, such as for example ArcGIS Pro (see https://pro.arcgis.com/en/pro-app/latest/help/editing/divide-a-polygon-by-a-value.htm) or QGIS with the Polygon Divider plugin (see https://github.com/jonnyhuck/RFCL-PolygonDivider). I'm not aware of any JS library than could do this on ...


4

Union the layer with itself Extract by expression "fid"<>"fid_2", replace fid with whatever id column you have.


1

Go to Vector < Data Management Tool < Split Vector Layer. There is one Advanced Parameters Select shp there.


1

You can create a "Virtual Layer" through Layer > Add Layer > Add/Edit Virtual Layer.... Replace polylayer with the name of your polygon layer and idcolumn with the name of your id column (create one if you don't have one): SELECT ST_INTERSECTION(a.geometry, b.geometry) FROM polylayer a, polylayer b WHERE a.idcolumn < b....


2

Let's assume there area two layers, a polygon layer 'layer' and a polyline layer 'lines_test', see the image below In the polygon layer create a new field "NewId" using Field calculator with the following expression: with_variable( 'lines', 'lines_test', -- name of the polyline layer array_find( array_agg($id, ...


1

As was written by @g12345kk in his question: I did find a way how to transform from MultiPolygons to Polygons without error: import geopandas as gpd small = gpd.read_file(r'R:\...\small.shp') large = gpd.read_file(r'R:\...\large.shp') smallClip = gpd.clip(small, large) smallClip_explode = smallClip['geometry'].explode() gdf = gpd.GeoDataFrame(geometry=...


0

I think the best way to do this is to add vertices to your line (for example every 2 meter) and then with python you can loop over the vertices and check for intersection with the parcellayer. If there is intersection, then you write the value of the counter (starting with 1) to a field of the intersected parcel. Next you check for the following vertex and ...


1

The reason why it happens can be the following: Input: Output: The resulting pink feature is the MULTIPOLYGON(((1)),((2))) To tackle this circumstances you can either: use the GeoDataFrame.explode(), see this thread for more details: Explanding GeoPandas Multipolygon Dataframe To One Poly Per Line use the GeoDataFrame.overlay() that contains a keyword ...


2

Thin can easily be done using QGIS expressions. To get a layer with actual geometries, use Menu Processing / Toolbox / Geomtry by expression (see here for details) and an expression like this (see first screenshot): make_triangle( make_point ( "attr_x" , "attr_y" ), make_point ( "attr_x"+2, "attr_y" ), ...


2

Here is a sample python code you can use. Starting scenario is that you have two layers points and polygons. After you create a new PointFeature you need to run the python script which will first drop all existing PolygonFeatures and then create a PolygonFeature at every point location using theadd_x and add_y constants. I am using EPSG:3857 so these values ...


2

Another solution is the "Collect geometries" geoalgorithm, can be found under Vector > Geometry Tools > Collect geometries. Takes a vector layer and collects its geometries into new multipart geometries. One or more attributes can be specified to collect only geometries belonging to the same class (having the same value for the specified ...


0

Your Option 4 is probably the closest to what I would suggest. What I would recommend is: Create linear features that represent the polygon outline and use tools that work with the linear features. eg: Parcel edge and road edge. You could create the linear features by using Polygon to line translation and then use Polygon Neighbors tool to identify the ...


0

I had a similar issue and after quite some investigation I found a solution! You need to use UNIQUE INDEX, which encapsulates hashed normalized spatial feature. CREATE TABLE GeomTable ( ID SERIAL PRIMARY KEY, Geom geometry ); CREATE UNIQUE INDEX UniqueGeometry ON GeomTable(SHA512(ST_Normalize(Geom))); Let's discuss the functions that we are using: ...


1

Yes, the returned values will be the median of the grid cells inside the buffer distance, and each point is independent so if two points are closer than the buffer distance then their returned value will be calculated including some set of common grid points. Its not like one point is "using up" some of the grid such that there's fewer left for the ...


0

I was finally able to do it by following the steps below (this works with larger and more complex polygon layers): Add an attribute with value 1 to the attribute table of all vector layers. This will be used later as a value for the burned areas that will be added with the overlapping values. (if you have the overlapping polygons in a single layer you need ...


3

Let's assume there is a polygon layer 'ne_10m_admin_0_countries_lakes' with features that have some "problematic" geometries (data was taken from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries_lakes.zip), see image below. Step 1. Use the "Geometry by expression" geoalgorithm ...


8

See: https://stackoverflow.com/questions/7595635/how-to-convert-polygon-data-into-line-segments-using-postgis : SELECT st_length(ST_MakeLine(sp,ep)) FROM -- extract the endpoints for every 2-point line segment for each linestring (SELECT ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp, ST_PointN(geom, generate_series(2, ...


1

I like to use PyQGIS (much easier than ArcPy) for problems like this. Zoom to your polygon and the raster overlapping current extent will be copied to a folder: import os from shutil import copyfile rasterfolder = r'/home/bera/GIS/Data/grid50m/Hillshades/' copy_to = r'/home/bera/GIS/Data/testdata/out/' mapextent = iface.mapCanvas().extent() for root, ...


6

Start editing Select the part Delete part Draw a polygon on top of the hole Select and merge


2

Another approach utilizes the .query(), that was also mentioned in the Use a list of values to select rows from a Pandas dataframe. import geopandas as gpd list_of_names = ['Type 1', 'Type 3'] _layer = "C:/Documents/Python Scripts/geopandas/layer.shp" layer = gpd.read_file(_layer) print(layer) fid Type ...


2

Use isin, see Use a list of values to select rows from a Pandas dataframe: import geopandas as gpd list_of_names = ['2262', '2361', '2161', '2039', '2132'] df = gpd.read_file(r"C:\GIS\data\testdata\ak_riks_lines.shp") #df.shape #Out: (294, 16) fieldname = 'KOM_KOD' df2 = df[df[fieldname].isin(list_of_names)] #df2.shape #(5, 16)


1

Probably long after the OP needs this, but here is an answer I think should be more straightforward than the existing ones. This uses QGIS 3.16. Dissolve the coral reef layer (Vector > Geoprocessing Tools > Dissolve) Run the Intersection tool (Vector > Geoprocessing Tools > Intersection) Input layer = dissolved coral reef layer from step 1 ...


3

I think you make things a bit complicated, because there is a very useful module called geopandas.overlay(). Let's assume there are two polygon layers 'grid' (comparable with your 'Small_polygon') and 'layer' (comparable with your 'Polygons'), see the image below. Using the following code: import geopandas as gpd _layer = "C:/Documents/Python Scripts/...


7

Use GeoPandas Overlay polygons = gpd.read_file("Polygons.shp") small_polygon = gpd.read_file("Small_polygon.shp") Intersection of the two GeoDataFrames: result = gpd.overlay(polygons,small_polygon, how='intersection') Result Areas of the intersection polygons result['area'] =result.apply(lambda row: row.geometry.area,axis=1)


1

library(raster) library(sf) poly <- st_read("directory.shp") # Creating a raster ##write the rows column of polygonn and xmn, xmx, ymn, ymx according to your country data r <- raster(nrows=, ncols=, xmn=, xmx=, ymn= , ymx=, crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,...


4

You're looking for a "Set Overlay" operation. Here's the way to do it in geopandas: poly_intersections = Polygons.overlay(Small_polygon, how='intersection') Here's some more documentation on how to work with it: https://geopandas.org/en/stable/docs/user_guide/set_operations.html Then you can get the area of the new cut/split polygons by doing this:...


0

I solved this problem using this code: var drawingTools = Map.drawingTools(); // Initialize empty geomtry layers as placeholder and fill them with drawn shapes while (drawingTools.layers().length() === 0) { var dummyLayer1 = ui.Map.GeometryLayer({geometries: null, name: 'geometries1', color: 'grey'}); var dummyLayer2 = ui.Map.GeometryLayer({geometries: ...


5

I don't use geopandas but in any GIS system what you are describing is a Union or Intersect style operation. This you can apply at the dataset level and I suspect would be far more efficient than looping over individual geometries as you are doing now. It appears that geopandas offers up such geoprocessing using overlay().


4

One way to do this would be to create a second polygon (you could use the Envelope Tool) in ArcMap. You may need to edit the new polygon around the area you have highlighted (where everything flows into the sea). Run the "Erase" tool in ArcTool box using the North American polygon and new envelope polygon you created. The output will be a polygon ...


1

Okay, I figured it out. Here is the sequence of actions in order to obtain the proper result : first, buffer your lines next, create a symetric difference finally, use 'multiple to individual' This will get you the individual street blocks, each as a single polygon


0

Method in ui.Map.DrawingTools for drawing objects is 'addLayer'; as it can be corroborated in following code. var polygon = ee.Geometry.Polygon( [[[51.50311161691894, 25.26703875176188], [51.50427033121338, 25.26502063795976], [51.50637318308105, 25.266146128641104],]], null, false); var line = ee.Geometry.LineString( [[...


0

The Plugins you use in QGIS do not get transfered into QField. There is no tracing in QField.


4

Another option is using the "Display" through RMC > Properties > Display. Paste the following expression format_number($area,2,'de') into the 'Display Name' and press Apply. Note: I made it in very intuitive way, more sophisticated output can be achieved via the 'HTML Map Tip'. After hover over the feature and get the output


7

I personally find the identify tool sometimes too clunky. Here another approach: To show the area for a polygon without using an additional field, you can just use a formula for a label with the $area function, e.g.: round($area/10000,1) will give you the area of the features in hectares. As you just want to see the area when pointing/clicking, you could ...


3

Also possible by means of PyQGIS. def feature_area(selected, _, __): n = len(selected) if n == 1: g = iface.activeLayer().selectedFeatures()[0].geometry() print("The area of selected feature is " + '{:0,.2f}'.format(g.area()) + " m²") elif n > 1: print("Select only one feature"...


Top 50 recent answers are included