Hot answers tagged

7

I think the speed problem you are having is looping around each feature and using geoprocessing tools inside the loop. They are not designed for that. They expect to process the whole dataset at once. So restructure your script to avoid the loop and it should complete "within the time to drink a cup of coffee". The data is in my opinion so small that it ...


6

Ryan, You should be able to with the SQLite dialect: $ ogrinfo -ro -geom=no -where "TZID='America/Boise' OR TZID='America/Denver'" tz_world.shp tz_world INFO: Open of `tz_world.shp' using driver `ESRI Shapefile' successful. Layer name: tz_world Metadata: DBF_DATE_LAST_UPDATE=2012-06-25 Geometry: Polygon Feature Count: 3 Extent: (-179.999900, -89.999900) ...


5

This sounds like a projection issue. You can try reprojecting to EPSG:3857 -- WGS84 Web Mercator (Auxiliary Sphere) and see what happens when you upload your next shapefile to CartoDB.


5

Problems: 1: the outcome of UnionSpatialPolygons is a spatial polygon 2: converting the result back into a spatial polygon data frame is a real pain -a. you need a very exact data frame to attach to a spatial polygon -b. data you used for UnionSpatialPolygons has more rows than the output and is not formatted in the way that is needed. My (ugly) ...


5

Assuming you downloaded and installed QGIS through the OSGeo4W Network Installer you will have access to the QGIS Processing Framework. Among other algorithm providers this will provide you access to the GDAL/OGR Provider which includes (under OGR miscellaneous) the Information tool (which is ogrinfo). The results will be output to a temporary HTML file. If ...


5

If pol.shp contains your polygons and pnt.shp - points: Create input.vrt: <OGRVRTDataSource> <OGRVRTLayer name="pol"> <SrcDataSource>pol.shp</SrcDataSource> <SrcLayer>pol</SrcLayer> </OGRVRTLayer> <OGRVRTLayer name="pnt"> <SrcDataSource>pnt.shp</SrcDataSource> <SrcLayer>pnt&...


5

Keep it simple, you don't need if f1 is Noneand others (assuming that you know the shapefiles to handle) in a simple personal script. Here is a simple solution with your data (with one polygon in each shapefile) from osgeo import ogr poly1 = ogr.Open('poly1.shp') poly2 = ogr.Open('poly2.shp') layer1 = poly1.GetLayer() layer1.GetFeatureCount() 1 # first ...


5

I do not understand why beginners try to start with the GDAL/OGR Python bindings (not very "Pythonic" and difficult) when other easier alternatives are available. With your script, you need to know osgeo.ogr and the SQL dialect of SQLite. The solution proposed by Mike T is powerful but not "classic" and performs only the intersection of shapefiles. What ...


4

This gave me endless hours of frustration, and the above answers (even the ones in the comments) didn't work for me either, but I found my own solution that did the trick. It turns out that to do an intersection or a clip between two layers, they must be saved on disk with the same CRS. So even if they look like they overlap in QGIS, it may be doing it's ...


4

Since you are using python you could also have a look at Shapely. It is a library for manipulation and analysis of geometric objects in the Cartesian plane. For your problem you would proceed like Load both geometries into Shapely objects on the obeject you want to substract the other from call object.difference(other) save the result into a new variable ...


4

For a shapefile in CRS 3857 (Pseudo Mercator) QGIS writes a prj projection file with the following WKT (formatted): PROJCS["WGS_84_Pseudo_Mercator", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["...


4

Here are the one-liners for you, please check them for field/filenames etc. Add a new field: ogrinfo input.shp -sql "ALTER TABLE input ADD COLUMN field3 integer(5)" Do the math (you will need your SQLite driver to be compiled with SQLite SQL dialect): ogrinfo input.shp -dialect SQLite -sql "UPDATE input SET field3 = field1/field2" I don't see where ...


4

I ran your script (slightly modified) at the Python Console of QGIS: from osgeo import ogr vlayer = iface.activeLayer() provider = vlayer.dataProvider() path = provider.dataSourceUri() tmp = path.split("|") path_to_shp_data = tmp[0] driver = ogr.GetDriverByName("ESRI Shapefile") dataSource = driver.Open(path_to_shp_data, 1) layer = dataSource.GetLayer(...


4

I don't understand your problem import osgeo.ogr as ogr import osgeo.osr as osr driver = ogr.GetDriverByName("ESRI Shapefile") srs = osr.SpatialReference() srs.ImportFromEPSG(4326) # create the data source data_source = driver.CreateDataSource("test_ogr.shp") # create the layer layer = data_source.CreateLayer("test_ogr", srs, ogr.wkbPoint) # Add the field ...


4

It is understandable to read only the manual page of ogr2ogr http://www.gdal.org/ogr2ogr.html and miss the page about options which are general to all OGR utilities http://www.gdal.org/ogr_utilities.html because the first one does not mention that the latter exists. However, you can save you -sql "<query>" terms into text files and reuse them by ...


4

You don't need Gdal/Python here. It is easier to use Shapely 1) Transform the list of points to a shapely geometry list = [(0.0, 0.0), (0.0, 1.0), (0.0, 2.0), (0.0, 3.0), (0.0, 4.0), (0.0, 5.0), (0.0, 6.0), (0.0, 7.0), (0.0, 8.0), (0.0, 9.0), (1.0, 0.0), (1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0), (1.0, 5.0), (1.0, 6.0), (1.0, 7.0), (1.0, 8.0), (1.0, 9....


4

The Gohlke GDAL/OGR wheel includes the FileGDB driver compiled as a plugin. To get the FileGDB driver working: Copy the Esri bin64\FileGDB.dll to [python install/virtualenv dir]\Lib\site-packages\osgeo (use bin\FileGDB.dll if using 32bit python). Do not copy the FileGDB.dll to the gdalplugins directory. Set GDAL_DRIVER_PATH environment variable, either:...


3

As asked by the OP, I put this in answer. The problem wasn't with the method used, as ogr.Geometry.Union() is the appropriate one but with the second polygon which was not closed. It can be done with : tri2.CloseRings() or by repeating the first point : tri2.AddPoint(1,1) If there is many features to dissolve (like the entire layer) it can be done with ...


3

I would utilize the GDAL SQLite dialect. Here comes an example: The two polygons in the image are as WKT: POLYGON (( -113 39, -113 42, -110 42, -110 39, -113 39 )) POLYGON (( -112 40, -112 41, -111 41, -111 40, -112 40 )) You can make a hole into polygon 1 with polygon 2 by using Spatialite function ST_Difference https://www.gaia-gis.it/gaia-...


3

IMHO homebrew will be best/easiest, once their build issues are fixed. They are on the case, so if you have a GitHub account you may want to subscribe to notifications on the relevant issue - that way you'll know the moment it is available: https://github.com/Homebrew/homebrew/pull/42845 I just tried the Python easy_install route and it blew up my meagre ...


3

Assuming you have an older version of GDAL/OGR, you can use OGR SQL to cast the geom field to a geometry. For example: ogr2ogr -f "ESRI Shapefile" "sample.shp" "geo.csv" -sql "SELECT *, CAST(geom as geometry) FROM geo" This will create a new shapefile using the WKT data as the geometry. You can use the same query with ogrinfo as well: ogrinfo -ro -sql "...


3

You can use PyQGIS code (I know that you use Ubuntu) at the Python Console of QGIS. To test my sugestion, I used the next code where a QgsExpression object has your class 1 values: layer = iface.activeLayer() expression = QgsExpression( u'"values" >= 30 AND "values" <= 50' ) idx = [ feat.attributes()[0] for feat in layer.getFeatures() if ...


3

When you tell GDAL to use SQLite SQL dialect the most part of the work is done by Spatialite library. Couple of weeks ago another user reported that ST_Collect in Spatialite was pretty slow for any bigger tasks https://groups.google.com/forum/#!topic/spatialite-users/t3QrvA1qQyg. Week after that Alessandro Furieri announced about a code optimization which ...


3

You should exchange ymin and ymax in -te (target extent). -te 102.3375307079206 10.350077240783799 107.6323189079206 14.6874037407838


3

I've written my own geojson map tiler driver in ruby. it's a quick and dirty ruby script to create geojson tiles from a postgis database using gdal/ogr. I could not find out how to cut the geojson directly, so I imported the json into a postgis database and exported them using ogr. ogr is supposed to read geojson but I could not get the driver to work ...


3

QGIS contains the standard Processing plugin which contains hundreds of in-built algorithms from a number of open source softwares such as GRASS, GDAL/OGR, SAGA etc. Using the Python Console in QGIS allows users to not only write pure python but also call the available algorithms from the Processing plugin and specify the required parameters without the ...


3

You weren't requesting the first feature, you were requesting a feature with an FID of "0" and your WFS doesn't have a feature with an FID of "0". Why: The ogr.Layer.__getitem__ method which allows element and slice access, calls the ogr.Layer.GetFeature method. GetFeature fetches a feature by its identifier, not by its position. This is intentional, so ...


3

You can do everything that is possible to do with SQL also with ogrinfo by utilizing the SQLite dialect http://www.gdal.org/ogr_sql_sqlite.html. Once you know that the problem turns into a question about how to make such a query and the answer can be found from Stackoverflow http://stackoverflow.com/questions/17194145/sql-count-based-on-column-value The ...


3

What is strange if the intersection result is a point ? The intersect predicate is Returns True if the boundary and interior of the object intersect in any way with those of the other. With a common point between the geometries, the intersects predicate returns TRUE because the boundary of the first geometry intersects the boundary of the second ...


3

The key information is here: I am also able to get OGR formats: ogrinfo --formats prints a bunch of them in the Windows cmd (with no FileGDB there though). This tells me that your GDAL_DRIVER_PATH environment variable is not set. See the instructions at https://trac.osgeo.org/gdal/wiki/FileGDB#Testingthedriver. Once "FileGDB" shows in ogrinfo --formats ...



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