Geospatial Analysis 地理空间分析

Your First Map

读取数据

import geopandas as gpd
world_loans = gpd.read_file(loans_filepath)
world = gpd.read_file(world_filepath)

展示数据

ax = world.plot(figsize=(20,20), color='whitesmoke', linestyle=':', edgecolor='black')#画出一张地图
world_loans.plot(ax=ax, markersize=2)#在地图上画数据

Interactive Maps

画普通地图

# Create a map
m_1 = folium.Map(location=[42.32,-71.0589], tiles='openstreetmap', zoom_start=10)

# Display the map
m_1

Plotting points

# Create a map
m_2 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=13)

# Add points to the map
for idx, row in daytime_robberies.iterrows():
    Marker([row['Lat'], row['Long']]).add_to(m_2)

# Display the map
m_2

folium.plugins.MarkerCluster

# Create the map
m_3 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=13)



# Add points to the map
mc = MarkerCluster()
for idx, row in daytime_robberies.iterrows():
    if not math.isnan(row['Long']) and not math.isnan(row['Lat']):
        mc.add_child(Marker([row['Lat'], row['Long']]))
m_3.add_child(mc)

# Display the map
m_3

Bubble maps

# Create a base map
m_4 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=13)

def color_producer(val):
    if val <= 12:
        return 'forestgreen'
    else:
        return 'darkred'

# Add a bubble map to the base map
for i in range(0,len(daytime_robberies)):
    Circle(
        location=[daytime_robberies.iloc[i]['Lat'], daytime_robberies.iloc[i]['Long']],
        radius=20,
        color=color_producer(daytime_robberies.iloc[i]['HOUR'])).add_to(m_4)
# Display the map
m_4

Heatmaps

# Create a base map
m_5 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=12)

# Add a heatmap to the base map
HeatMap(data=crimes[['Lat', 'Long']], radius=10).add_to(m_5)

# Display the map
m_5

Choropleth maps

# GeoDataFrame with geographical boundaries of Boston police districts
districts_full = gpd.read_file('../input/geospatial-learn-course-data/Police_Districts/Police_Districts/Police_Districts.shp')
districts = districts_full[["DISTRICT", "geometry"]].set_index("DISTRICT")
districts.head()
# Number of crimes in each police district
plot_dict = crimes.DISTRICT.value_counts()
plot_dict.head()
# Create a base map
m_6 = folium.Map(location=[42.32,-71.0589], tiles='cartodbpositron', zoom_start=12)

# Add a choropleth map to the base map
Choropleth(geo_data=districts.__geo_interface__, 
           data=plot_dict, 
           key_on="feature.id", 
           fill_color='YlGnBu', 
           legend_name='Major criminal incidents (Jan-Aug 2018)'
          ).add_to(m_6)

# Display the map
m_6

Manipulating Geospatial Data

from geopy.geocoders import Nominatim
#从名字反查地理位置
geolocator = Nominatim(user_agent="kaggle_learn")
location = geolocator.geocode("Pyramid of Khufu")

print(location.point)
print(location.address)

Attribute join

# Use an attribute join to merge data about countries in Europe 按照 name 匹配合并
europe = europe_boundaries.merge(europe_stats, on="name")
europe.head()

Spatial join

#按照地理位置信息匹配,就是地图上的位置有交集
european_universities = gpd.sjoin(universities, europe)

Proximity Analysis 附近位置分析

计算距离

# Select one release incident in particular
recent_release = releases.iloc[360]

# Measure distance from release to each station
# 这里返回的应该是一个Series
distances = stations.geometry.distance(recent_release.geometry)
distances
# 获取距离这个点最近的station
print(stations.iloc[distances.idxmin()][["ADDRESS", "LATITUDE", "LONGITUDE"]])

创建一个缓冲区

# 以各个点为中心创建类似足球的12边形
two_mile_buffer = stations.geometry.buffer(2*5280)
two_mile_buffer.head()


# Plot each polygon on the map 打印出这些区域
GeoJson(two_mile_buffer.to_crs(epsg=4326)).add_to(m)
# 汇集成一个区域 Type: <class 'shapely.geometry.multipolygon.MultiPolygon'>
my_union = two_mile_buffer.geometry.unary_union
# 是否包含一个位置
my_union.contains(releases.iloc[360].geometry)