geopandas
Python tools for geographic data
I've got two GeoDataFrames that clearly intersect (see below), but when I run the intersects method (or intersection, or touches) I get zero counts/empty geometries.
Any guidance about what I am doing wrong (or Geopandas is doing wrong) would be very helpful.
Source: (StackOverflow)
Is it possible to get counts of intersections between two geometries using GeoPandas objects? That is, I want to count up the number of polygons or line strings in one GeoDataFrame that intersect with each polygon in another GeoDataFrame. I did not see an easy way of doing this while browsing the GeoPandas docs, but wanted to check before moving on to lower-level tools.
Source: (StackOverflow)
I've been trying to use the "intersects" feature on a geodataframe, looking to see which points lie inside a polygon. However, only the first feature in the frame will return as true. What am I doing wrong?
from geopandas.geoseries import *
p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)
g1 = GeoSeries([p1,p2,p3])
g2 = GeoSeries([p2,p3])
g = GeoSeries([Polygon([(0,0), (0,2), (2,2), (2,0)])])
g1.intersects(g) # Flags the first point as inside, even though all are.
g2.intersects(g) # The second point gets picked up as inside (but not 3rd)
Source: (StackOverflow)
I have tried to install geopandas two different ways: pip install geopandas or by cloning
git clone https://github.com/kjordahl/geopandas
In both cases, the installation file setup.py runs for a while and then returns this error message:
src/fiona/ogrinit.c:300:23: fatal error: cpl_error.h: No such file or directory
compilation terminated.
error: command 'gcc' failed with exit status 1
fiona is interface to OGR so Python can read geospatial data. cpl_error.h seems to be missing. What can I do?
Source: (StackOverflow)
I suspect this is a trivial question to answer but how do I load mapinfo mid/mif files into geopandas?
While it's trivial to load and manipulate .shp files, I can't work out what code to use:
Here's a few examples of the data I'm trying to play with: http://www.stoke.gov.uk/ccm/content/business/general/open-geospatial-consortium-data-catalogue.en
cylpaths = geopandas.GeoDataFrame.from_file(path)
I tried to test the geometry, which shows:
cylpaths.crs
Out[15]:
{u'ellps': u'airy',
u'k': 0.9996012717,
u'lat_0': 49,
u'lon_0': -2,
u'no_defs': True,
u'proj': u'tmerc',
u'towgs84': u'375,-111,431,-0,-0,-0,0',
u'units': u'm',
u'x_0': 400000,
u'y_0': -100000}
But when I try to preview the data, I get an empty dataframe:
cylpaths
Out[16]:
Int64Index([], dtype='int64') Empty GeoDataFrame
I'm a bit lost in working out what I need to do next.
Source: (StackOverflow)
I am at a loss to get the following code to work. For whatever reason, the GeoPandas *.plot() doesn't work, but I want to use both Pandas and GeoPandas for some simple plots.
I have been trying to take the Shapely objects from GeoPandas and plot them on a Basemap. The problem is that the polygons won't plot. I iterate through them from the GeoPandas.geometry, add them to the axes collection, and then use plot() - to no avail. Basemap seems to work fine, the code doesn't give any errors, but the polygons - counties - don't appear...
Thank you for the help!
import geopandas as gpd
from descartes import PolygonPatch
import matplotlib as mpl
import mpl_toolkits.basemap as base
import matplotlib.pyplot as plt
counties_file = r'C:\Users\...\UScounties\UScounties.shp'
counties = gpd.read_file(counties_file)
#new plot
fig = plt.figure(figsize=(5,5),dpi=300)
#ax = fig.add_subplot(111)
ax = ax = plt.gca()
minx, miny, maxx, maxy = counties.total_bounds
#map
m = base.Basemap(llcrnrlon=minx, llcrnrlat=miny,
urcrnrlon=maxx, urcrnrlat=maxy,
resolution='h', area_thresh=100000,
projection='merc')
patches = []
#add polygons
for poly in counties.geometry:
#deal with single polygons and multipolygons
if poly.geom_type == 'Polygon':
p = PolygonPatch(poly, facecolor='blue', alpha=1)
#plt.gca().add_patch(p)
#ax.add_patch(p)
patches.append(p)
elif poly.geom_type == 'MultiPolygon':
for single in poly:
q = PolygonPatch(single,facecolor='red', alpha=1)
#ax.add_patch(p)
patches.append(q)
m.drawcoastlines(linewidth=.1)
m.fillcontinents()
m.drawcountries(linewidth=.25,linestyle='solid')
m.drawstates(linewidth=.25,linestyle='dotted')
m.drawmapboundary(fill_color='white')
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))
ax.plot()
plt.show()
Source: (StackOverflow)
Building on data from this question on faceting through looping, I was wondering if it was possible to call a ax = df.plot(kind='bar') and assign the thus generated AxesSubplot object to a specific axis position/coordinate? (like facet row 1, col 1, 2, 3, etc ... )?
The reason I am asking is not really for faceting bar plots as such, but for faceting map production using the geopandas library. If it worked with bar plots, it might also work with geopandas geodataframe.plot() calls. I can't plot a map from axes themselves, therefore I seem to need to go the other way around--get axes as a byproduct of a plot call, and then place that as appropriate in a grid.
Non-working example--the loop is really pseudo here; I don't move the axis index to plot a different panel each time (in fact, I overwrite the axes object from the subplots call). That is, though, what I would like to do--map the axis object generated from the plot call to the axes (coordinate space) from the subplots call).
N = 100
industry = ['a','b','c']
city = ['x','y','z']
ind = np.random.choice(industry, N)
cty = np.random.choice(city, N)
jobs = np.random.randint(low=1,high=250,size=N)
df_city =pd.DataFrame({'industry':ind,'city':cty,'jobs':jobs})
## how many panels do we need?
cols =df_city.city.value_counts().shape[0]
fig, axes = plt.subplots(1, cols, figsize=(8, 8))
for x, city in enumerate(df_city.city.value_counts().index.values):
data = df_city[(df_city['city'] == city)]
data = data.groupby(['industry']).jobs.sum()
axes = data.plot(kind='bar')
print type(axes)
fig.suptitle('Employment By Industry By City', fontsize=20)
<class 'matplotlib.axes._subplots.AxesSubplot'>
<class 'matplotlib.axes._subplots.AxesSubplot'>
<class 'matplotlib.axes._subplots.AxesSubplot'>
Source: (StackOverflow)
I am a newbie to Python, and I have an issue that I believe is due to circular dependencies, but I have been unable to resolve it.
How can I figure out where the circular dependency is occurring, and what can I do to resolve this error (please find the terminal output below): (I use OsX Yosemite, if that helps)
python l3_l0.py
Traceback (most recent call last):
File "l3_l0.py", line 18, in <module>
import geopandas as gpd
File "build/bdist.macosx-10.10-intel/egg/geopandas/__init__.py", line 6, in <module>
File "build/bdist.macosx-10.10-intel/egg/geopandas/geoseries.py", line 12, in <module>
Thank you!
Source: (StackOverflow)
I'm very new to projections and struggling to convert a map dataset from OSGB36 (Eastings and Northings) to WGS84 (lat/longs).
I've used geopandas to import a shapefile. Now I want to export it as GeoJSON so I can then use Folium to map it on leaflet.js.
If I check the crs variable with .crs I get:
{u'datum': u'OSGB36',
u'k': 0.999601272,
u'lat_0': 49,
u'lon_0': -2,
u'no_defs': True,
u'proj': u'tmerc',
u'units': u'm',
u'x_0': 400000,
u'y_0': -100000}
Trouble is if I'm a bit lost on the attributes to set to convert it to WGS84:
staffs.to_crs({'datum':"wgs84"})
RuntimeError: projection not named
Can you help?
Source: (StackOverflow)
I'm trying to make a choropleth map from polygons in a Geopandas GeoDataFrame. I want to symbolize the polygons by quantiles of a value in one of the GeoDataFrame columns. I'm trying to figure out different options and see what best fits my needs. Any advice on this would be greatly appreciated.
It appears that Geopandas does have some ability to do this already: http://nbviewer.ipython.org/github/geopandas/geopandas/blob/master/examples/choropleths.ipynb
tracts.plot(column='CRIME', scheme='QUANTILES', k=3, colormap='OrRd')
This works, although I can't find much documentation. I'd like to be able to add a legend that shows the quantile cut-off values, but it appears that the Geopandas plot currently only allows for legends based of categorical data. Does anyone have a work around for this?
Additionally, I'd like to be able to adjust the polygon outline width. Is this possible?
As an alternative option I've been playing around with is using polygon patches in matplotlib. This appears much more involved but does seem to offer more options to customize. If necessary to go down this route in order to build a legend I can follow-up with another question and will include my code so far.
Thanks for the help.
Source: (StackOverflow)
I am using GeoPandas and Pandas.
I have a, say, 300,000 rows Dataframe, df, with 4 columns + the index column.
id lat lon geometry
0 2009 40.711174 -73.99682 0
1 536 40.741444 -73.97536 0
2 228 40.754601 -73.97187 0
however the unique ids are only a handful (~200)
I want to generate a shapely.geometry.point.Point object for each (lat,lon) combination, similarly to what shown here: http://nbviewer.ipython.org/gist/kjordahl/7129098
(see cell#5),
where it loops through all rows of the dataframe; but for such a big dataset, I wanted to limit the loop to the much smaller number of unique ids.
Therefore, for a given id value, idvalue (i.e., 2009 from the first row) create the GeoSeries, and assign it directly to ALL rows that have id==idvalue
My code looks like:
for count, iunique in enumerate(df.if.unique()):
sc_start = GeoSeries([Point(np.array(df[df.if==iunique].lon)[0],np.array(df[df.if==iunique].lat)[0])])
df.loc[iunique,['geometry']] = sc_start
however things don't work - the geometry field does not change - and I think is because the indexes of sc_start don't match with the indexes of df.
how can I solve this? should I just stick to the loop through the whole df?
Source: (StackOverflow)
I'm trying to fix an error I have when I try to use the geocode function in geopandas.
from geopandas.geocode import geocode
df['latlong'] = geocode(df.Location, provider="mapquest")
This is what I see:
/Users/.../lib/python2.7/site-packages/geopandas-0.1.0.dev_- py2.7.egg/geopandas/geocode.pyc in geocode(strings, provider, **kwargs)
72 'mapquest': geopy.geocoders.MapQuest,
73 'openmapquest': geopy.geocoders.OpenMapQuest,
---> 74 'nominatim' : geopy.geocoders.Nominatim}
75
76 if provider not in coders:
AttributeError: 'module' object has no attribute 'Nominatim'
I did try to install the plugin using the instructions here but with no joy: https://github.com/rdeguzman/python-nominatim
Source: (StackOverflow)
I have an ESRI shapefile that I have loaded as a geopandas dataframe, and it has a geometry column containing polygons and multipolygons.
I have no idea how to operate on this column, though. For example, for each entry, I would like to count the number of pixels.
I played around with the 'apply' function, but to no avail....
mygeodf['geometry'].apply('sum')
Source: (StackOverflow)
I'm looking to do the equivalent of the ArcPy Generate Near Table using Geopandas / Shapely. I'm very new to Geopandas and Shapely and have developed a methodology that works, but I'm wondering if there is a more efficient way of doing it.
I have two point file datasets - Census Block Centroids and restaurants. I'm looking to find, for each Census Block centroid, the distance to it's closest restaurant. There are no restrictions in terms of same restaurant being the closest restaurant for multiple blocks.
The reason this becomes a bit more complicated for me is because the Geopandas Distance function calculates elementwise, matching based on index. Therefore, my general methodology is to turn the Restaurants file into a multipoint file and then set the index of the blocks file to all be the same value. Then all of the block centroids and the restaurants have the same index value.
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon, Point, MultiPoint
Now read in the Block Centroid and Restaurant Shapefiles:
Blocks=gpd.read_file(BlockShp)
Restaurants=gpd.read_file(RestaurantShp)
Since the Geopandas distance function calculates distance elementwise, I convert the Restaurant GeoSeries to a MultiPoint GeoSeries:
RestMulti=gpd.GeoSeries(Restaurants.unary_union)
RestMulti.crs=Restaurants.crs
RestMulti.reset_index(drop=True)
Then I set the index for the Blocks equal to 0 (the same value as the Restaurants multipoint) as a work around for the elementwise calculation.
Blocks.index=[0]*len(Blocks)
Lastly, I use the Geopandas distance function to calculate the distance to the nearest restaurant for each Block centroid.
Blocks['Distance']=Blocks.distance(RestMulti)
Please offer any suggestions on how any aspect of this could be improved. I'm not tied to using Geopandas or Shapely, but I am looking to learn an alternative to ArcPy.
Thanks for the help!
Source: (StackOverflow)
Good afternoon,
I'm trying to use a combination of geopandas, pandas and folium to create a polygon map that I can embed incorporate into a web page.
For some reason, it's not displaying and wonder if anyone can help.
The steps I've taken:
1) Grabbed a .shp from the UK's OS for Parliamentary boundaries.
2) I've then used geopandas to change the projection to epsg=4326 and then exported as GeoJSON which takes the following format:
{ "type": "Feature", "properties": { "PCON13CD": "E14000532", "PCON13CDO": "A03", "PCON13NM": "Altrincham and Sale West" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ -2.313999519326579, 53.357408280545918 ], [ -2.313941776174758, 53.358341455420039 ], [ -2.31519699483377, 53.359035664493433 ], [ -2.317953152796459, 53.359102954309151 ], [ -2.319855973429864, 53.358581917200119 ],... ] ] ] } },...
Then what I'd like to do is mesh this with a dataframe of constituencies in the following format, dty:
constituency count
0 Burton 667
1 Cannock Chase 595
2 Cheltenham 22
3 Cheshire East 2
4 Congleton 1
5 Derbyshire Dales 1
6 East Staffordshire 4
import folium
mapf = folium.Map(width=700, height=370, tiles = "Stamen Toner", zoom_start=8, location= ["53.0219392","-2.1597434"])
mapf.geo_json(geo_path="geo_json_shape2.json",
data_out="data.json",
data=dty,
columns=["constituency","count"],
key_on="feature.properties.PCON13NM.geometry.type.Polygon",
fill_color='PuRd',
fill_opacity=0.7,
line_opacity=0.2,
reset="True")
The output from mapf looks like:
mapf.json_data
{'../../Crime_data/staffs_data92.json': [{'Burton': 667,
'Cannock Chase': 595,
'Cheltenham': 22,
'Cheshire East': 2,
'Congleton': 1,
'Derbyshire Dales': 1,
'East Staffordshire': 4,
'Lichfield': 438,
'Newcastle-under-Lyme': 543,
'North Warwickshire': 1,
'Shropshire': 17,
'South Staffordshire': 358,
'Stafford': 623,
'Staffordshire Moorlands': 359,
'Stoke-on-Trent Central': 1053,
'Stoke-on-Trent North': 921,
'Stoke-on-Trent South': 766,
'Stone': 270,
'Tamworth': 600,
'Walsall': 1}]}
Although the mapf.create_map() function successfully creates a map, the polygons don't render.
Can anyone suggest any debugging steps?
I'm never clear on how to add the full data files if anyone needs them, so please let me know.
Thank you for your time in advance.
Source: (StackOverflow)