在地理 pandas 中移动阿拉斯加和夏威夷以获取脉络膜
问题描述
在R中，我可以像这样移动阿拉斯加和夏威夷： https://www.storybench.org/how-to-shift-alaska-and-hawaii-below-the-lower-48-for-your-interactive-choropleth-map/
我可以在Cartopy中这样做： Showing Alaska and Hawaii in Cartopy map
如何使用Geopandas执行类似操作？
这是我的……
import requests
import geopandas
import numpy as np
np.random.seed(5)
URL = 'https://github.com/kjhealy/us-county/raw/master/data/geojson/gz_2010_us_040_00_500k.json'
us_json = requests.get(URL).json()
# save to disk as fname then load from file
usa = geopandas.read_file(fname)
usa["value"] = np.random.uniform(low=0,high=100, size=52)
usa.sort_values("NAME", inplace=True)
usa
现在绘制它：
usa.plot(column='value', cmap='Reds', scheme='quantiles', k=7, legend=True)
且仅限内地：
not_mainland = ['Alaska', 'Hawaii', "Puerto Rico"]
mainland_usa = usa.query('NAME not in @not_mainland')
mainland_usa.sort_values("NAME")
mainland_usa.plot(column='value', cmap='Reds', scheme='quantiles', k=7, legend=True)
如何添加阿拉斯加和夏威夷的次要情节？
解决方案
对于连续状态，可以照常使用geoDataFrame.plot()。但阿拉斯加和夏威夷必须作为插图地图单独绘制。在下面完整且可运行的代码中，有许多插入的注释，它们解释了帮助读者的重要步骤。关于专题数据，我更喜欢使用更真实的人口数据(而不是问题中的化妆数据)来展示地块上的专题地图。包mapclassify用于对数据进行分类，以便进行适当的主题数据处理。
import matplotlib.pyplot as plt
import matplotlib
import cartopy
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import shapely.geometry as sgeom #for box drawing
import geopandas as gpd
import numpy as np
import mapclassify as mc
import pandas as pd
import requests
import json
# population per sq-km
# include `Puerto Rico` but not used
popdensity = {
'Alaska': 0.8,
'District of Columbia': 4251,
'Hawaii': 86,
'Puerto Rico': 360,
'New Jersey': 438.00,
'Rhode Island': 387.35,
'Massachusetts': 312.68,
'Connecticut': 271.40,
'Maryland': 209.23,
'New York': 155.18,
'Delaware': 154.87,
'Florida': 114.43,
'Ohio': 107.05,
'Pennsylvania': 105.80,
'Illinois': 86.27,
'California': 83.85,
'Virginia': 69.03,
'Michigan': 67.55,
'Indiana': 65.46,
'North Carolina': 63.80,
'Georgia': 54.59,
'Tennessee': 53.29,
'New Hampshire': 53.20,
'South Carolina': 51.45,
'Louisiana': 39.61,
'Kentucky': 39.28,
'Wisconsin': 38.13,
'Washington': 34.20,
'Alabama': 33.84,
'Missouri': 31.36,
'Texas': 30.75,
'West Virginia': 29.00,
'Vermont': 25.41,
'Minnesota': 23.86,
'Mississippi': 23.42,
'Iowa': 20.22,
'Arkansas': 19.82,
'Oklahoma': 19.40,
'Arizona': 17.43,
'Colorado': 16.01,
'Maine': 15.95,
'Oregon': 13.76,
'Kansas': 12.69,
'Utah': 10.50,
'Nebraska': 8.60,
'Nevada': 7.03,
'Idaho': 6.04,
'New Mexico': 5.79,
'South Dakota': 3.84,
'North Dakota': 3.59,
'Montana': 2.39,
'Wyoming': 1.96}
# use this simple colormap
my_colormap = matplotlib.cm.Reds
# some settings
edgecolor = "gray"
# use this column for thematic mapping
theme_value = "pop_per_sqkm"
# A function that draws inset map
# ===============================
def add_insetmap(axes_extent, map_extent, state_name, facecolor, edgecolor, geometry):
# create new axes, set its projection
use_projection = ccrs.Mercator() # preserves shape
#use_projection = ccrs.PlateCarree() # large distortion in E-W, bad for for Alaska
geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
sub_ax = plt.axes(axes_extent, projection=use_projection) # normal units
sub_ax.set_extent(map_extent, geodetic) # map extents
# option to add basic land, coastlines of the map
# can comment out if you don't need them
sub_ax.add_feature(cartopy.feature.LAND)
sub_ax.coastlines()
sub_ax.set_title(state_name)
# add map `geometry`
sub_ax.add_geometries([geometry], ccrs.PlateCarree(),
facecolor=facecolor, edgecolor=edgecolor, lw=0.3)
# +++ more features can be added here +++
# plot box around the map
extent_box = sgeom.box(map_extent[0], map_extent[2], map_extent[1], map_extent[3])
sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), color='none')
# access USA shapefile
# use data from internet
URL = 'https://github.com/kjhealy/us-county/raw/master/data/geojson/gz_2010_us_040_00_500k.json'
# more hard-code for settings
state_name = "NAME" #defined column header
fname = "usa_json.json"
# Only on the first-run, set True in the next statement
if True:
us_json = requests.get(URL).json()
with open(fname, 'w') as file:
file.write(json.dumps(us_json))
# saved it to disk as fname
# then load from file
usa = gpd.read_file(fname)
usa.sort_values(state_name, inplace=True)
# make dataframe from population data `popdensity` dict object
usa_popden = pd.DataFrame.from_dict(popdensity, orient='index',
columns=['pop_per_sqkm'])
# merge `usa_popden` to main dataframe, `usa`
newusa = usa.merge(usa_popden, how='left', left_on='NAME', right_index=True)
# take only some columns in `newusa` for our operation
newusa = newusa[['NAME', 'pop_per_sqkm', 'geometry']]
# Data classification for thematic mapping
# choose number of classes of population density
# classes --> assigned colors in thematic mapping
num_classes = 7 #will get class-values: 0,1,2,3,4,5,6
sclass = mc.Quantiles(newusa["pop_per_sqkm"].values, k=num_classes)
#print(sclass)
# add new column, "sclass", for raw class values, sclass.yb
# its values will be used to assign color for polygon's facecolor
newusa["sclass"] = sclass.yb
# extract parts of the whole 'newusa' geodataframe for separate plotting/manipulation
# 'usa_main': excluding non-conterminous states
usa_main = newusa[~newusa[state_name].isin(["Alaska", "Hawaii", "Puerto Rico"])] # exclude these
# re-project usa_main to equal-area conic projection "EPSG:2163"
usa_main.crs = {'init': 'epsg:4326'}
usa_main = usa_main.to_crs(epsg=2163)
# 'usa_more': non-conterminous states, namely, Alaska and Hawaii
usa_more = newusa[newusa[state_name].isin(["Alaska", "Hawaii"])] # include these
# ------------ Plot --------------
# plot 1st part, using usa_main and grab its axis as 'ax2'
ax2 = usa_main.plot(column="sclass", legend=False,
cmap=matplotlib.cm.Reds, ec=edgecolor, lw=0.4)
# manipulate colorbar/legend
fig = ax2.get_figure()
cax = fig.add_axes([0.9, .25, 0.02, 0.5]) #[left,bottom,width,height]
sm = plt.cm.ScalarMappable(cmap=my_colormap,
norm=plt.Normalize(vmin=min(newusa["sclass"]),vmax=max(newusa["sclass"])))
# clear the array of the scalar mappable
sm._A = []
cb = fig.colorbar(sm, cax=cax)
cb.set_label("Pop-density class")
# manipulate the axis seetings
ax2.set_frame_on(False)
ax2.set_xticks([])
ax2.set_yticks([])
ax2.set_xticklabels([])
ax2.set_yticklabels([])
ax2.set_title("Population Density Plot")
# add more features on ax2
# plot Alaska, Hawaii as inset maps
for index,state in usa_more.iterrows():
if state[state_name] in ("Alaska", "Hawaii"):
st_name = state[state_name]
# set fill color, using normalized `sclass` on `my_colormap`
facecolor = my_colormap( state["sclass"] / max(newusa["sclass"] ))
if st_name == "Alaska":
# (1) Alaska
# Custom extent, relative size
map_extent = (-178, -135, 46, 73) # degrees: (lonmin,lonmax,latmin,latmax)
axes_extent = (0.04, 0.06, 0.29, 0.275) # axes units: 0 to 1, (LLx,LLy,width,height)
if st_name == "Hawaii":
# (2) Hawaii
# Custom extent, relative size
map_extent = (-162, -152, 15, 25)
axes_extent = (0.27, 0.06, 0.15, 0.15)
# add inset maps
add_insetmap(axes_extent, map_extent, st_name,
facecolor,
edgecolor,
state["geometry"])
plt.show()
数据分类：
在上面的代码中，sclass携带了人口密度的分类详细信息。要获取类列表，请运行以下代码：
low_val = 0
max_cval = len(sclass.bins)-1
print("Population density (per square_km)")
print("Class Value_ranges")
for ix,val in enumerate(sclass.bins):
#print(low_val, "< x <=", val)
print("{ix:} {low_val:.2f} < density <= {val:.2f}".format(ix=ix, low_val=low_val, val=val))
#print(my_colormap(ix/max_cval))
low_val = val
并且您应该会得到下面的结果，该结果应该与上面的地图一起出现。
Population density (per square_km)
Class Value_ranges
0 0.00 < density <= 7.48
1 7.48 < density <= 18.56
2 18.56 < density <= 30.50
3 30.50 < density <= 51.70
4 51.70 < density <= 75.38
5 75.38 < density <= 155.09
6 155.09 < density <= 4251.00
额外的色条操作
可以操纵色条将标签文本从简单的类号0、1、2、.6更改为分类值范围。只需在上述代码的plt.show()行之前插入这些代码行，然后重新运行。
class_txts = []
low_val = 0
for ix,val in enumerate(sclass.bins):
class_txts.append("{low_val:.1f}, {val:.1f}".format(ix=ix, low_val=low_val, val=val))
low_val = val
cb.ax.set_yticklabels( class_txts )
cb.set_label("Pop-density range")
新绘图将如下所示：
相关文章