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

To process large files, you need to use a generator which only read/write one line at a time and Python has a command that does that: with. The with statement handles opening and closing the file, including if an exception is raised in the inner block. The for line in f treats the file object f as an iterable, which automatically uses buffered IO and ...


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

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


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

Python does not contain a separate library for GDAL but rather bindings to access the GDAL libraries. As GDAL is updated so are the affected bindings, ensuring that you'll have full access to GDALs functionality from Python. If you upgrade GDAL you will also upgrade the Python bindings (if the version you upgrade to supports Python bindings). If you build ...


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


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

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


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


3

A python solution is fairly simple using ogr.GetDriverByName(in_format).Open(in_file). import sys from osgeo import ogr def main(in_file, in_format, out_file, out_format): in_ds = ogr.GetDriverByName(in_format).Open(in_file) out_ds = ogr.GetDriverByName(out_format).CopyDataSource(in_ds, out_file) if __name__ == '__main__': main(*sys.argv[1:]) ...


3

Since OGR version 1.10.0 the sqlite SQL dialect has been able to be applied to any spatial datset. Which is great, as it means that you can apply it to your GeoJSON files. Looking at the OGR GeoJSON documentation you can see that the layer name for a GeoJSON file is OGRGeoJSON which means that the SQL that selects from the GeoJSON file will translate from ...


3

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


3

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


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

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

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


3

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


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


2

I just ran into the same problem. I ended using the ExportToJson function in ogr and then reading the Json string into a dictionary. Using my data and the notation from the original question, this looks like: import json ... ring_dict = json.loads(ring.ExportToJson()) ring_dict {'coordinates': [[-4.94237, 55.725449], [-4.941922, 55.725585], ...


2

Just an update on this question using the new GDAL OpenFileGDB Driver this is possible. As mentioned in this answer you'll first have to unzip the .lpk file which is actually a 7zip archive. With Homebrew $ brew install p7zip And to unzip $ 7z x yourLayerPackage.lpk Now that you have the uncompressed .gdb file you can use install GDAL/OGR to convert ...


2

With the OGR's function UnionCascaded() from osgeo import ogr shp= "yourfile" shapef = ogr.Open(shpor) lyr = shapef.GetLayer() unionc = ogr.Geometry(ogr.wkbMultiPolygon) for feat in xrange(lyr.GetFeatureCount()): fit= lyr.GetFeature(feat) geom= fit.GetGeometryRef() unionc.AddGeometry(geom) union= unionc.UnionCascaded()


2

In addition to @Ali's answer... My guess is the Clone() method calls for a deep copy of the object. Where the other method provides a shallow copy. The differences come as one only copies the memory pointer while the other copies the entirety and stores it in another memory location. Since OGR is a C++ library, you can check out this tutorial on C++ ...


2

From my experience the difference is subtle in that GetGeometryRef() will store a reference to the geometry which will be removed when the underlying feature is destroyed. If you use Clone() when storing the geometry you can access it and its methods even after destroying the feature. To highlight this, suppose you want to open a shapefile, read a feature ...


2

Just figured out a way to do this in R, inspired by the link I posted in the comments (which uses outdated functionalities but the right packages); as an example, I'll use the "Southwesternmost" counties in Michigan (shapefile here--Allegan, Berrien, Cass, Kalamazoo, St. Joseph, and Van Buren counties) So, ordering counties alphabetically, the (length>0, ...



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