Nearly perfect answer from @ThingumaBob
Only thing to extend ist that in his script the geometry type will only by geometry instead of Multipolygon. So i extended your example with the ST_Multi() function:
CREATE TABLE export.osm_adminareas_a_clip5 AS
ST_Multi(ST_Intersection(a.geom, ST_MakeEnvelope(12, 47, 13, 48, 4326)))::...
@Vince gave an educational comment, @Cyril named it directly, and I like to avoid further confusion; the simple answer is:
clip those grometries!
All queries present in this post do one thing, one way or another, that is selecting every row whose geometries minimum bounding rectangle, the bbox, intersects the created envelope, another bbox.
And that's ...
Try the following with both the overlap condition and the intersects
CREATE TABLE export.osm_adminareas_a AS
WHERE geom && ST_MakeEnvelope(12, 47, 13, 48, 4326)
AND ST_INTERSECTS(geom, ST_MakeEnvelope(12,47,13,48, 4326))
The reason why I include the overlap condition is that the ST_INTERSECTS ...
I noticed later that, even being only one feature, there were tiny islands around the boundary, and maybe the issue was in those. But i wasn't able to delete all of them, because i was in a hurry.
So, to solve my issue, i used Simplify, with "Distance (Douglas-Peucker)" method, using a 1km tolerance. Any type of polygon based directly on the original (...
If you create a clipping layer with only the circular feature to use, your code can be simplified in this way (I use my own paths):
from osgeo import gdal, ogr
OutTile = gdal.Warp("/home/zeito/pyqgis_data/cut.tif",
You can also use raster calculator for this:
Select 'Raster' > 'Raster Calculator...' use the expression ("Mask_raster@1"=1)*"Raster1@1" + ("Mask_raster@1"=0)*-9999 (replace Mask_raster and Raster1 with your layer names).
Select 'Raster' > 'Conversion' > 'Translate (Convert Format)...'. Set the input layer as the layer you just created in raster calculator ...
You can clip a raster by a vector easily, so first step is to convert your raster to vector.
GDAL > Raster Conversion > Polygonize.
You will have polygons with attribute called DN with values 0 and 1. Select all polygons with value 1 with expression like DN = 1
Use GDAL > Raster Extraction > Clip raster by mask layer. Select the polygonized layer as the ...
It is not clear for me what do you want to do with above code but, when you have overlapped features in the same layer you can use itertools python module for avoiding this kind of lines:
if outerFeature.id() != innerFeature.id()
for repeated features. So, your code can be written as follow for producing a memory layer of all intersection of primary ...
ogr2ogr -f "ESRI Shapefile" out.shp road.shp -dialect sqlite -sql "select ST_Diffrence(road.geometry, seg_buff.geometry) from road, 'seg_buff.shp'.seg_buff seg_buff
Put an extra seg_buff at the end
The -spat parameter in ogr2ogr https://gdal.org/programs/ogr2ogr.html#ogr2ogr affects all available layers. This selects data withing the given rectangle from all layers if you write only the name of the datasource into the command but not the layername. So just "...input.gpkg" instead of "...input.gpkg selected_layer".
ogr2ogr -f gpkg -spat 100000 6000000 ...