Hot answers tagged

36

Fiona returns Python dictionaries and you can not use poly['properties']['score']) += point['properties']['score']) with a dictionary. Example of summing attributes using the references given by Mike T: # read the shapefiles import fiona from shapely.geometry import shape polygons = [pol for pol in fiona.open('poly.shp')] points = [pt for pt in fiona.open(...


33

If you've got an GDAL/OGR dev environment (headers, libs), you could radically simplify your code by using Fiona. To read features from a shapefile, add new attributes, and write them out as GeoJSON is just a handful of lines: import fiona import json features = [] crs = None with fiona.collection("docs/data/test_uk.shp", "r") as source: for feat in ...


30

I've designed Fiona (an OGR wrapper) to make this kind of processing simple. from fiona import collection import logging log = logging.getLogger() # A few functions to shift coords. They call eachother semi-recursively. def shiftCoords_Point(coords, delta): # delta is a (delta_x, delta_y [, delta_y]) tuple return tuple(c + d for c, d in zip(coords,...


21

Using JEQL This can be done with three lines: ShapefileReader t file: "shapefile.shp"; out = select * except (GEOMETRY), Geom.translate(GEOMETRY,100,100) from t; ShapefileWriter out file: "ahapefile_shift.shp";


21

The question is about Fiona and Shapely and the other answer using GeoPandas requires to also know Pandas. Moreover GeoPandas uses Fiona to read/write shapefiles. I do not question here the utility of GeoPandas, but you can do it directly with Fiona using the standard module itertools, specially with the command groupby ("In a nutshell, groupby takes an ...


18

You can use the shape function of Shapely: from shapely.geometry import shape c = fiona.open('data/boroughs/boroughs_n.shp') pol = c.next() geom = shape(pol['geometry']) and a MultiPolygon is a list of Polygons,so Multi = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('data/boroughs/boroughs_n.shp')]) Example with one of my data: # the ...


18

Why would you want install Fiona if you can use PyQGIS to read the shapefile attributes ? Fiona is for reading geometries and attributes of a shapefile file (as PyShp), therefore you don't need dbfpy (look at Python Script examples for geoprocessing shapefiles without using arcpy). Fiona is a Python module so you must install it as any Python module in the ...


16

I have reproduced your example with shapefiles. You can use Shapely and Fiona to solve your problem. 1) Your problem (with a shapely Point): 2) starting with an arbitrary line (with an adequate length): from shapely.geometry import Point, LineString line = LineString([(point.x,point.y),(final_pt.x,final_pt.y)]) 3) using shapely.affinity.rotate to ...


15

you can get a list of the drivers with >>> import fiona >>> fiona.supported_drivers {'ESRI Shapefile': 'raw', 'ARCGEN': 'r', 'PCIDSK': 'r', 'SUA': 'r', 'DGN': 'raw', 'SEGY': 'r', 'MapInfo File': 'raw', 'GeoJSON': 'rw', 'PDS': 'r', 'FileGDB': 'raw', 'GPX': 'raw', 'DXF': 'raw', 'GMT': 'raw', 'Idrisi': 'r', 'GPKG': 'rw', 'OpenFileGDB': 'r',...


14

Happily OGR can do this for you as both ogr.Feature and ogr.Geometry objects have ExportToJson() methods. In your code; fe.ExportToJson() And since geojson FeatureCollection objects are simply dictionaries with a type of FeatureCollection and a features object containing a list of Feature objects. feature_collection = {"type": "FeatureCollection", ...


14

I'd do it like this: def explode(coords): """Explode a GeoJSON geometry's coordinates object and yield coordinate tuples. As long as the input is conforming, the type of the geometry doesn't matter.""" for e in coords: if isinstance(e, (float, int, long)): yield coords break else: for f in ...


14

Your solution is a little out of date (look at Fiona - Preffered method for defining a schema). For a better solution look at that proposed by Sean Gillies in gistfile1.py to parse a delimited text file data and create a new shapefile 1) How do I define the spatial reference: Fiona crs Module from fiona.crs import from_epsg from_epsg(3857) {'init': 'epsg:...


14

You are on the right track and the geopandas GeoDataFrame is a good choice for rasterization over Fiona. Fiona is a great toolset, but I think that the DataFrame is better suited to shapefiles and geometries than nested dictionaries. import geopandas as gpd import rasterio from rasterio import features Set up your filenames shp_fn = '...


13

And though the post is tagged with python, since JEQL has already been mentioned, here's an example with JavaScript (using GeoScript). /** * Shift all coords in all features for all layers in some directory */ var Directory = require("geoscript/workspace").Directory; var Layer = require("geoscript/layer").Layer; // offset for all geometry coords var dx =...


13

Using GDAL >= 1.10.0 compiled with SQLite and SpatiaLite: ogr2ogr data_shifted.shp data.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,1,10) FROM data" where shiftX = 1 and shiftY = 10.


13

Suppose we have two polygons (green and blue): They are not equal (as Fetzer suggest): green.equals(blue) False and blue.equals(green) False And we can can determine the difference (in red): blue.difference(green) and green.difference(blue) gives an empty geometry Thus, you can use a supplementary condition: if not poly1.difference(poly2).is_empty:...


13

Additionally - geopandas now optionally includes rtree as a dependency, see the github repo So, instead of following all the (very nice) code above, you could simply do something like: import geopandas from geopandas.tools import sjoin point = geopandas.GeoDataFrame.from_file('point.shp') # or geojson etc poly = geopandas.GeoDataFrame.from_file('poly.shp')...


13

To change projections with Fiona, use the pyproj module. Example with a point shapefile (you can simplify the algorithm): from pyproj import Proj, transform import fiona from fiona.crs import from_epsg shape = fiona.open('sample.shp') original = Proj(shape.crs) # EPSG:4326 in your case destination = Proj(init='EPSG:...') # your new EPSG with fiona.open('...


13

The question is about Shapely and Fiona in pure Python without QGIS ("using command line and/or shapely/fiona"). A solution is from shapely import shape, mapping import fiona # schema of the new shapefile schema = {'geometry': 'Polygon','properties': {'area': 'float:13.3','id_populat': 'int','id_crime': 'int'}} # creation of the new shapefile with the ...


12

That's the gist of it. The R-tree allows you to make a very fast first pass and gives you a set of results that will have "false positives" (bounding boxes may intersect when the geometries precisely do not). Then you go over the set of candidates (fetching them from the shapefile by their index) and do a mathematically precise intersection test using, e.g., ...


12

You can pass the json directly to the GeoDataFrame constructor: import geopandas as gpd import requests data = requests.get("https://data.cityofnewyork.us/api/geospatial/arq3-7z49?method=export&format=GeoJSON") gdf = gpd.GeoDataFrame(data.json()) gdf.head() Outputs: features type 0 {'type': '...


11

fc.next() is a simple iterator: fc = fiona.open("my.shp") first_feature = fc.next() second_feature = fc.next() ... Or more simply: for feat in fiona.open("my.shp") print feat The result is a Python dictionary. For example with one result (feat=) {'geometry': {'type': 'Point', 'coordinates': (180627.0, 330190.0)}, 'type': 'Feature', 'id': '154', '...


11

I highly recommend GeoPandas for dealing with large assortments of features and performing bulk operations. It extends Pandas dataframes, and uses shapely under the hood. from geopandas import GeoSeries, GeoDataFrame # define your directories and file names dir_input = '/path/to/your/file/' name_in = 'cb_2013_us_county_20m.shp' dir_output = '/path/to/...


10

Use Fiona of Sean Gillies , a very simple wrapper of the OGR library (The Fiona User Manual) All the elements of a shapefile (schema, records) are processed using Python dictionaries: schema of one of my shapefiles as example: {'geometry': 'LineString', 'properties': {u'faille': 'str:20', u'type': 'str:20', u'id': 'int'}} one record in the shapefile: ...


10

You need to use the function shape of Shapely: from shapely.geometry import shape c = fiona.open("ne_10m_admin_0_countries.shp") # first record country = c.next() print "country name :",country['properties']['NAME'] country name : Aruba # shape(country['geometry']) -> shapely geometry print "bounds:", shape(country['geometry']).bounds bounds: (-70....


10

At the first question, 'epsg:32054' code has feet units. For this reason, it is necessary to use 'preserve_units=True' as parameter in inProj = Proj(init='epsg:32054') line. Now, next code works well: from pyproj import Proj, transform # wisconsing EPSG:32054 # epsg:4326 is for the entire world, wgs 84...not obvious inProj = Proj(init='epsg:32054', ...


9

Use Rtree as an index to perform the much faster joins, then Shapely to do the spatial predicates to determine if a point is actually within a polygon. If done properly, this can be faster than most other GISes. See examples here or here. The second part of your question concerning 'SUM', use a dict object to accumulate populations using a polygon id as ...


9

A MultiLineString is a list of lines: from shapely.geometry import MultiLineString, mapping, shape coords = [((0, 0), (1, 1)), ((-1, 0), (1, 0))] lines = MultiLineString(coords) print lines MULTILINESTRING ((0 0, 1 1), (-1 0, 1 0)) for line in lines: print line LINESTRING (0 0, 1 1) LINESTRING (-1 0, 1 0) # convert to GeoJSON format: print mapping(...


9

1) With Fiona, you don't need shapely to count the number of points in a polygon/multipolygon. Simply use the resulting GeoJSON format (= a Python dictionary). Polygon simple: >>> import fiona >>> shape = fiona.open("simplePoly.shp") >>> # first feature >>> feature = shape.next() >>> geom = feature['geometry'] &...


9

You've almost got it, but you've made a small error. You need to use the intersection method on the spatial index, rather than passing the index to the intersection method on the buffered point. Once you've found a list of features where the bounding boxes overlap, then you need to check if your buffered point actually intersects the geometries. import ...


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