Hot answers tagged

6

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


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

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


2

I believe your second solution is way slower that the first one, since Clip has an inner Selection By Location plus other geometrical operations. So I suggest some improvements to your first algorithm: As per @artwork21 comment, merge the shapefiles if possible. e.g. they share the same schema. This procedure removes an extra loop over the shapefiles, ...


2

Solved using Programmatic raster-vector calculation I had some troubles using directly gdal.RasterizeLayer() on layers like propose here in the cookbook but it seems that using "MEM" data source and finally writing it on the disk is maybe better. My solution (which is one of many possible solutions) def rasterizer(shapePath, rasterPath, attribute, ...


2

You should be able to just use the sql argument in ogr2ogr. For instance, with the following polygon Shapefile with two features and two attributes: $>ogrinfo -so -al polygon.shp INFO: Open of `polygon.shp' using driver `ESRI Shapefile' successful. Layer name: polygon Geometry: Polygon Feature Count: 2 Extent: (-1.206294, -0.828671) - (0.727273, ...


2

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


2

Interpolate shapes using tool from 3d analyst. Run field calculator expression on numeric field type float, parser - Python def getMedian(shp): aList=[] part=shp.getPart(0) n=len(part) for i in xrange(n): p=part.getObject(i) z=p.Z aList.append(z) aList.sort() return aList[n/2] To call it use: ...


1

Some things to think about. If your example is any indicator, you are going to have over 5 million records when you are finished. Neither shapefiles or personal geodatabases will handle that very well, so you should use either a file geodatabase or SDE. Processing 48 million records won't be quick, especially when it involves a spatial selection. Perhaps ...


1

I thought I'd add my experience here. I've been working with situations where I might need to merge hundreds or thousands of polygons and the Union solution didn't quite cut it - the time to get the intersection of 1 polygon in 1 layer against a couple of hundred polygons in 1 other layer took ~15 secs. I managed to cut it down to 0.08ms by adding each ...


1

Tidiest approach that I have seen is to use python to trigger the pgsql2shp command line utility which is part of the postgis package. Not very pythonic, but very simple and robust. e.g. See How can I get a shapefile from a postgis query? Once you have a shapefile it is straight forward to convert to most common formats using ogr2ogr, also over the ...



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