Python 第七周作业

本文章使用工具从 Jupyter Notebook 转换而来。
本文将在课程结束后撤下。

Data sources

https://mapshaper.org/ for shapefile
https://tj.jiangsu.gov.cn/2024/index.htm for population data
Thanks to Google's Gemini 2.5 Pro for converting HTML to CSV

Do note that assets files are .gitignored and will not be included.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
import shapefile
import csv
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
import numpy as np
import matplotlib.colors as mcolors

plt.rcParams['font.sans-serif'] = ['HarmonyOS Sans SC']  # Huawei make good fonts lol

r = shapefile.Reader('assets/shapefiles/js.shp')

geom = r.shape(0)
print(geom)
print(geom.shapeTypeName)
print(geom.shapeType)
print(geom.bbox)
print(geom.points)
print(r.fields)
print(r.records())

city_to_data = {}
with open('assets/js-data.csv', 'r', encoding='utf-8') as f:
    reader = csv.DictReader(f)
    for row in reader:
        city_to_data[row['city']] = {
            'population': float(row['population']),
            'gdp': float(row['gdp']),
            'lat': float(row['lat']),
            'lon': float(row['lon'])
        }

city_field_name = 'name'

w = shapefile.Writer('assets/shapefiles/js-edited.shp', shapeType=r.shapeType, encoding='gbk')
w.fields = r.fields[:]

w.field('lat', 'F', 10, 6)
w.field('lon', 'F', 10, 6)
w.field('population', 'F', 15, 0)
w.field('gdp', 'F', 15, 2)
w.field('per_gdp', 'F', 15, 2)

for i in range(len(r)):
    w.shape(r.shape(i))
    old_rec = r.record(i)

    city_index = [f[0] for f in r.fields[1:]].index(city_field_name)
    city = old_rec[city_index]

    lat = 0.0
    lon = 0.0
    people = 0.0
    gdp = 0.0
    per_gdp = 0.0

    if city in city_to_data:
        data = city_to_data[city]
        lat = data['lat']
        lon = data['lon']
        people = data['population']
        gdp = data['gdp']
        if people > 0:
            per_gdp = gdp / people
        else:
            per_gdp = 0.0

    new_rec = old_rec + [lat, lon, people, gdp, per_gdp]
    w.record(*new_rec)

w.close()

1
2
3
4
5
6
7
    Polygon #0
    POLYGON
    5
    (118.364371, 31.229373, 119.241664, 32.608232)
    [(119.067452, 32.462398), (119.047022, 32.462882), (119.035909, 32.452428), (119.02373, 32.456073), (119.022495, 32.440303), (119.029455, 32.428836), (119.016434, 32.413193), (119.017556, 32.399261), (119.02749, 32.388316), (119.037425, 32.388844), (119.038435, 32.38159), (119.022551, 32.380447), (119.02171, 32.375436), (119.043823, 32.365191), (119.027322, 32.35134), (119.028445, 32.345095), (119.041971, 32.333264), (119.035853, 32.325435), (119.035572, 32.308718), (119.039838, 32.303746), (119.041073, 32.271709), (119.037088, 32.259868), (119.056563, 32.25168), (119.074748, 32.248951), (119.086029, 32.241158), (119.132332, 32.245164), (119.153772, 32.243668), (119.177232, 32.238912), (119.229036, 32.222179), (119.241664, 32.216278), (119.221403, 32.201127), (119.216239, 32.190686), (119.191095, 32.1854), (119.184529, 32.189629), (119.154446, 32.186193), (119.13671, 32.193374), (119.121332, 32.187603), (119.086871, 32.178527), (119.078677, 32.178747), (119.076039, 32.161695), (119.047976, 32.161254), (119.058472, 32.15645), (119.039333, 32.157023), (119.031756, 32.14865), (119.022383, 32.130886), (119.00925, 32.12467), (119.008913, 32.115941), (119.026649, 32.115941), (119.049772, 32.109768), (119.057686, 32.102228), (119.069248, 32.107916), (119.079799, 32.107343), (119.089733, 32.089704), (119.098882, 32.091027), (119.098489, 32.07162), (119.08659, 32.053224), (119.092652, 32.037604), (119.097928, 32.011256), (119.092371, 32.004105), (119.10281, 32.004899), (119.112464, 32.000441), (119.114204, 31.988212), (119.121612, 31.984812), (119.121107, 31.977571), (119.113923, 31.979823), (119.111061, 31.96949), (119.102979, 31.965603), (119.090126, 31.974126), (119.073738, 31.969622), (119.06487, 31.97342), (119.046854, 31.968562), (119.02923, 31.95708), (119.029455, 31.950013), (119.037537, 31.93897), (119.051512, 31.935922), (119.06762, 31.940207), (119.104157, 31.933624), (119.108423, 31.922977), (119.111173, 31.901634), (119.116449, 31.890497), (119.107806, 31.885149), (119.100341, 31.865787), (119.092427, 31.860437), (119.075421, 31.864195), (119.069753, 31.868793), (119.043374, 31.854424), (119.032822, 31.853761), (119.025919, 31.846376), (119.003357, 31.845757), (118.979672, 31.842352), (118.969457, 31.834613), (118.970411, 31.826077), (118.986407, 31.81785), (118.998137, 31.80166), (119.006387, 31.798607), (119.002908, 31.783431), (118.982534, 31.782988), (118.980458, 31.764667), (118.998867, 31.76719), (119.00504, 31.777766), (119.021541, 31.778607), (119.029679, 31.784714), (119.055665, 31.788962), (119.065768, 31.78197), (119.077105, 31.783962), (119.077835, 31.773916), (119.089902, 31.769624), (119.093999, 31.756036), (119.105392, 31.751344), (119.114428, 31.740233), (119.12498, 31.735806), (119.129245, 31.723586), (119.138674, 31.722125), (119.157701, 31.699362), (119.186437, 31.694136), (119.190422, 31.687049), (119.186605, 31.677658), (119.192274, 31.659847), (119.187335, 31.649168), (119.204734, 31.646642), (119.212423, 31.627718), (119.232516, 31.632283), (119.234087, 31.627053), (119.221852, 31.59478), (119.216239, 31.585423), (119.21579, 31.570389), (119.197662, 31.551849), (119.189075, 31.55176), (119.181105, 31.538096), (119.180656, 31.526959), (119.198616, 31.530021), (119.199458, 31.524519), (119.183518, 31.500199), (119.1664, 31.493453), (119.152986, 31.493186), (119.149114, 31.483998), (119.156017, 31.480135), (119.152481, 31.464374), (119.141481, 31.451674), (119.140919, 31.444969), (119.155568, 31.438618), (119.164155, 31.441105), (119.169207, 31.421739), (119.169207, 31.394015), (119.174089, 31.380906), (119.190085, 31.379084), (119.2003, 31.368817), (119.204509, 31.357838), (119.214107, 31.355883), (119.218035, 31.348859), (119.198336, 31.330674), (119.193172, 31.318935), (119.195361, 31.309373), (119.181105, 31.300478), (119.176222, 31.293939), (119.15815, 31.294784), (119.149787, 31.276945), (119.140695, 31.276011), (119.13222, 31.262752), (119.120041, 31.261684), (119.10702, 31.250603), (119.10528, 31.235293), (119.083728, 31.239744), (119.073794, 31.23289), (119.062344, 31.238498), (119.048818, 31.233557), (119.014526, 31.241524), (119.0033, 31.233735), (118.986463, 31.231154), (118.984218, 31.237474), (118.92495, 31.2397), (118.915184, 31.244195), (118.891219, 31.237964), (118.869666, 31.242237), (118.840481, 31.236361), (118.831277, 31.229462), (118.819491, 31.236005), (118.806357, 31.229507), (118.79575, 31.229373), (118.781382, 31.239833), (118.778575, 31.261684), (118.761401, 31.277702), (118.726603, 31.28224), (118.707914, 31.296786), (118.704434, 31.302791), (118.719588, 31.294784), (118.721103, 31.32267), (118.735078, 31.344769), (118.745742, 31.372684), (118.75461, 31.385084), (118.765386, 31.387928), (118.770886, 31.376951), (118.767575, 31.363973), (118.784525, 31.368595), (118.80024, 31.368995), (118.829481, 31.377529), (118.847385, 31.393927), (118.851145, 31.392549), (118.86961, 31.421561), (118.866074, 31.44488), (118.870508, 31.463042), (118.883193, 31.496471), (118.886279, 31.527669), (118.873539, 31.532062), (118.875111, 31.551982), (118.881958, 31.565156), (118.873315, 31.578417), (118.877299, 31.589902), (118.862258, 31.595001), (118.870789, 31.599346), (118.868768, 31.611139), (118.858554, 31.62395), (118.832848, 31.630377), (118.819547, 31.628516), (118.814495, 31.619252), (118.802765, 31.619473), (118.79603, 31.628294), (118.789576, 31.646554), (118.782953, 31.656258), (118.799005, 31.666361), (118.774141, 31.682619), (118.768304, 31.677481), (118.760503, 31.680271), (118.748099, 31.675753), (118.745125, 31.656524), (118.7394, 31.651074), (118.735976, 31.633436), (118.72868, 31.634233), (118.719419, 31.627452), (118.71695, 31.639463), (118.694949, 31.639729), (118.685127, 31.636139), (118.673846, 31.640527), (118.65785, 31.641103), (118.643594, 31.649656), (118.647972, 31.65936), (118.643257, 31.671367), (118.646962, 31.679873), (118.654931, 31.67562), (118.671432, 31.688201), (118.674407, 31.696483), (118.680749, 31.696616), (118.697699, 31.709416), (118.687372, 31.716147), (118.684902, 31.724914), (118.677269, 31.728589), (118.653977, 31.730139), (118.652125, 31.73979), (118.641686, 31.743996), (118.648365, 31.748821), (118.641686, 31.758249), (118.627037, 31.759444), (118.619068, 31.75006), (118.609358, 31.757187), (118.600995, 31.754133), (118.593082, 31.744438), (118.57529, 31.740941), (118.571698, 31.746253), (118.55312, 31.729032), (118.529772, 31.742358), (118.529379, 31.748555), (118.538696, 31.749352), (118.53965, 31.761569), (118.533869, 31.767234), (118.523599, 31.762011), (118.510746, 31.765641), (118.502608, 31.778342), (118.482178, 31.778342), (118.486107, 31.795421), (118.504741, 31.841733), (118.490597, 31.843148), (118.475443, 31.850312), (118.466744, 31.858005), (118.467081, 31.868881), (118.4723, 31.879623), (118.452263, 31.884663), (118.44014, 31.892707), (118.41595, 31.901457), (118.402986, 31.914007), (118.386036, 31.914803), (118.375316, 31.925319), (118.364371, 31.930178), (118.369198, 31.94299), (118.368974, 31.956196), (118.379974, 31.968341), (118.376719, 31.981721), (118.387158, 31.983267), (118.390975, 31.989536), (118.389628, 32.01205), (118.401526, 32.018671), (118.38598, 32.033147), (118.394342, 32.048414), (118.385025, 32.060768), (118.393276, 32.062665), (118.405455, 32.080794), (118.433518, 32.086661), (118.441095, 32.092041), (118.462871, 32.10064), (118.473647, 32.11065), (118.472974, 32.114839), (118.489474, 32.120746), (118.500812, 32.121584), (118.494694, 32.130974), (118.500587, 32.14173), (118.507435, 32.137895), (118.495648, 32.165176), (118.495592, 32.177073), (118.50648, 32.182096), (118.509848, 32.194299), (118.521859, 32.188219), (118.532354, 32.195797), (118.539594, 32.192449), (118.551717, 32.199629), (118.561708, 32.196898), (118.575627, 32.202536), (118.580902, 32.19866), (118.598357, 32.199497), (118.619629, 32.205443), (118.626588, 32.202536), (118.644099, 32.20998), (118.667504, 32.235918), (118.675137, 32.255862), (118.667447, 32.257535), (118.659253, 32.266119), (118.670759, 32.285617), (118.667335, 32.296354), (118.659983, 32.296354), (118.660263, 32.305418), (118.678448, 32.316549), (118.696352, 32.317736), (118.683948, 32.324071), (118.687652, 32.330009), (118.703199, 32.328954), (118.707914, 32.33841), (118.697531, 32.342632), (118.699832, 32.349757), (118.692928, 32.351692), (118.702021, 32.357936), (118.696071, 32.367961), (118.698878, 32.383085), (118.690178, 32.378293), (118.678672, 32.393943), (118.685744, 32.404975), (118.688887, 32.418202), (118.687091, 32.431297), (118.692648, 32.465122), (118.68928, 32.47351), (118.672892, 32.469865), (118.647691, 32.47654), (118.643426, 32.47048), (118.628104, 32.467669), (118.608572, 32.470567), (118.608067, 32.475837), (118.59673, 32.476628), (118.592127, 32.482029), (118.592464, 32.500207), (118.608011, 32.501744), (118.617496, 32.512411), (118.614465, 32.5204), (118.604812, 32.52299), (118.608853, 32.537034), (118.598975, 32.54467), (118.584663, 32.549014), (118.56805, 32.558755), (118.563728, 32.564371), (118.569509, 32.586435), (118.598077, 32.601347), (118.601781, 32.597356), (118.623838, 32.592444), (118.633043, 32.578979), (118.641237, 32.587093), (118.658411, 32.594812), (118.668794, 32.5949), (118.689561, 32.588541), (118.700842, 32.588541), (118.69798, 32.60639), (118.725312, 32.608232), (118.725986, 32.596523), (118.740691, 32.589374), (118.743273, 32.597751), (118.757528, 32.603978), (118.775769, 32.597751), (118.774815, 32.587663), (118.784749, 32.582444), (118.801811, 32.58433), (118.805291, 32.593453), (118.820501, 32.604241), (118.824542, 32.602531), (118.824654, 32.589637), (118.820501, 32.582532), (118.827573, 32.581172), (118.827011, 32.572662), (118.842502, 32.572969), (118.848507, 32.568187), (118.855972, 32.574373), (118.872641, 32.576084), (118.874998, 32.58297), (118.884203, 32.580426), (118.885101, 32.586348), (118.901995, 32.589023), (118.909347, 32.586874), (118.896382, 32.577268), (118.898627, 32.560466), (118.887907, 32.55485), (118.914005, 32.548839), (118.919786, 32.557263), (118.936568, 32.557526), (118.940609, 32.547785), (118.964799, 32.526721), (118.967942, 32.514562), (118.97984, 32.503851), (118.997463, 32.503895), (119.00751, 32.514079), (119.01924, 32.517986), (119.041409, 32.514782), (119.043318, 32.502841), (119.050502, 32.492919), (119.062793, 32.484576), (119.067452, 32.462398), (118.86585, 31.518972), (118.855691, 31.504726), (118.846711, 31.50304), (118.846262, 31.50912), (118.856421, 31.522256), (118.86585, 31.518972)]
    [Field(name="DeletionFlag", field_type=FieldType.C, size=1, decimal=0), Field(name="adcode", field_type=FieldType.N, size=6, decimal=0), Field(name="name", field_type=FieldType.C, size=12, decimal=0), Field(name="center", field_type=FieldType.N, size=0, decimal=0), Field(name="centroid", field_type=FieldType.N, size=0, decimal=0), Field(name="childrenNu", field_type=FieldType.N, size=2, decimal=0), Field(name="level", field_type=FieldType.C, size=4, decimal=0), Field(name="parent", field_type=FieldType.N, size=0, decimal=0), Field(name="subFeature", field_type=FieldType.N, size=2, decimal=0), Field(name="acroutes", field_type=FieldType.N, size=0, decimal=0)]
    [Record #0: [320100, '南京市', None, None, 11, 'city', None, 0, None], Record #1: [320200, '无锡市', None, None, 7, 'city', None, 1, None], Record #2: [320300, '徐州市', None, None, 10, 'city', None, 2, None], Record #3: [320400, '常州市', None, None, 6, 'city', None, 3, None], Record #4: [320500, '苏州市', None, None, 9, 'city', None, 4, None], Record #5: [320600, '南通市', None, None, 7, 'city', None, 5, None], Record #6: [320700, '连云港市', None, None, 6, 'city', None, 6, None], Record #7: [320800, '淮安市', None, None, 7, 'city', None, 7, None], Record #8: [320900, '盐城市', None, None, 9, 'city', None, 8, None], Record #9: [321000, '扬州市', None, None, 6, 'city', None, 9, None], Record #10: [321100, '镇江市', None, None, 6, 'city', None, 10, None], Record #11: [321200, '泰州市', None, None, 6, 'city', None, 11, None], Record #12: [321300, '宿迁市', None, None, 5, 'city', None, 12, None]]

1
2
3
4
5
new_r = shapefile.Reader('assets/shapefiles/js-edited.shp', encoding='gbk')
print("\n新的属性表字段:")
print(new_r.fields)
print("\n新的属性表记录:")
print(new_r.records())

1
2
3
4
5
    新的属性表字段:
    [Field(name="DeletionFlag", field_type=FieldType.C, size=1, decimal=0), Field(name="adcode", field_type=FieldType.N, size=6, decimal=0), Field(name="name", field_type=FieldType.C, size=12, decimal=0), Field(name="center", field_type=FieldType.N, size=0, decimal=0), Field(name="centroid", field_type=FieldType.N, size=0, decimal=0), Field(name="childrenNu", field_type=FieldType.N, size=2, decimal=0), Field(name="level", field_type=FieldType.C, size=4, decimal=0), Field(name="parent", field_type=FieldType.N, size=0, decimal=0), Field(name="subFeature", field_type=FieldType.N, size=2, decimal=0), Field(name="acroutes", field_type=FieldType.N, size=0, decimal=0), Field(name="lat", field_type=FieldType.F, size=10, decimal=6), Field(name="lon", field_type=FieldType.F, size=10, decimal=6), Field(name="population", field_type=FieldType.F, size=15, decimal=0), Field(name="gdp", field_type=FieldType.F, size=15, decimal=2), Field(name="per_gdp", field_type=FieldType.F, size=15, decimal=2)]

    新的属性表记录:
    [Record #0: [320100, '南京市', None, None, 11, 'city', None, 0, None, 31.85, 118.82, 954, 17421.4, 18.25], Record #1: [320200, '无锡市', None, None, 7, 'city', None, 1, None, 31.57, 120.31, 749, 15456.19, 20.62], Record #2: [320300, '徐州市', None, None, 10, 'city', None, 2, None, 34.27, 117.23, 902, 8900.44, 9.87], Record #3: [320400, '常州市', None, None, 6, 'city', None, 3, None, 31.8, 119.97, 537, 10116.36, 18.82], Record #4: [320500, '苏州市', None, None, 9, 'city', None, 4, None, 31.3, 120.63, 1295, 24653.37, 19.03], Record #5: [320600, '南通市', None, None, 7, 'city', None, 5, None, 32.07, 120.87, 774, 11813.27, 15.25], Record #6: [320700, '连云港市', None, None, 6, 'city', None, 6, None, 34.6, 119.22, 459, 4363.61, 9.5], Record #7: [320800, '淮安市', None, None, 7, 'city', None, 7, None, 33.62, 119.15, 453, 5015.06, 11.05], Record #8: [320900, '盐城市', None, None, 9, 'city', None, 8, None, 33.35, 120.15, 668, 7403.87, 11.07], Record #9: [321000, '扬州市', None, None, 6, 'city', None, 9, None, 32.45, 119.42, 458, 7423.26, 16.19], Record #10: [321100, '镇江市', None, None, 6, 'city', None, 10, None, 119.45, 32.2, 322, 5264.07, 16.32], Record #11: [321200, '泰州市', None, None, 6, 'city', None, 11, None, 32.45, 119.92, 450, 6731.66, 14.94], Record #12: [321300, '宿迁市', None, None, 5, 'city', None, 12, None, 33.97, 118.3, 497, 4398.07, 8.84]]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
fig_pop, ax_pop = plt.subplots(figsize=(10, 10))
patches_pop = []
populations = []

people_index = [f[0] for f in new_r.fields[1:]].index('population')

for i in range(len(new_r)):
    shape = new_r.shape(i)
    rec = new_r.record(i)
    pop = rec[people_index]
    populations.append(pop)

    parts = list(shape.parts) + [len(shape.points)]

    for j in range(len(parts) - 1):
        start = parts[j]
        end = parts[j + 1]
        points = shape.points[start:end]
        if len(points) > 2:
            poly = Polygon(points, closed=True)
            patches_pop.append(poly)

if not patches_pop:
    print("No valid polygons to plot.")
else:
    # Normalize population values for color mapping
    pop_array = np.array(populations)
    if np.max(pop_array) > np.min(pop_array):
        norm_pop = plt.Normalize(vmin=np.min(pop_array), vmax=np.max(pop_array))
    else:
        norm_pop = plt.Normalize(vmin=0, vmax=1)

    pcollection_pop = PatchCollection(patches_pop, cmap='Reds', alpha=0.8)
    pcollection_pop.set_array(pop_array)
    pcollection_pop.set_norm(norm_pop)

    ax_pop.add_collection(pcollection_pop)

    bbox = new_r.bbox
    ax_pop.set_xlim(bbox[0], bbox[2])
    ax_pop.set_ylim(bbox[1], bbox[3])

    fig_pop.colorbar(pcollection_pop, ax=ax_pop, label='人口数')

    ax_pop.set_title('江苏省人口分布图')
    ax_pop.set_xlabel('经度')
    ax_pop.set_ylabel('纬度')
    ax_pop.set_aspect('equal')

    plt.show()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
fig_per, ax_per = plt.subplots(figsize=(10, 10))
patches_per = []
per_gdps = []

per_gdp_index = [f[0] for f in new_r.fields[1:]].index('per_gdp')

for i in range(len(new_r)):
    shape = new_r.shape(i)
    rec = new_r.record(i)
    per = rec[per_gdp_index]
    per_gdps.append(per)

    parts = list(shape.parts) + [len(shape.points)]

    for j in range(len(parts) - 1):
        start = parts[j]
        end = parts[j + 1]
        points = shape.points[start:end]
        if len(points) > 2:
            poly = Polygon(points, closed=True)
            patches_per.append(poly)

if not patches_per:
    print("No valid polygons to plot.")
else:
    # Convert to array
    per_array = np.array(per_gdps)

    # Determine classification: assume 5 levels, based on min and max linear classification
    if np.max(per_array) > np.min(per_array):
        levels = 5
        boundaries = np.linspace(np.min(per_array), np.max(per_array), levels + 1)
        norm_per = mcolors.BoundaryNorm(boundaries, ncolors=plt.get_cmap('Blues').N, clip=True)
    else:
        norm_per = plt.Normalize(vmin=0, vmax=1)
        boundaries = np.array([np.min(per_array), np.max(per_array) + 1])  # 虚拟边界

    pcollection_per = PatchCollection(patches_per, cmap='Blues', alpha=0.8)
    pcollection_per.set_array(per_array)
    pcollection_per.set_norm(norm_per)
    ax_per.add_collection(pcollection_per)

    bbox = new_r.bbox
    ax_per.set_xlim(bbox[0], bbox[2])
    ax_per.set_ylim(bbox[1], bbox[3])

    fig_per.colorbar(pcollection_per, ax=ax_per, label='人均GDP', boundaries=boundaries, ticks=boundaries)

    ax_per.set_title('江苏省人均GDP分布图')
    ax_per.set_xlabel('经度')
    ax_per.set_ylabel('纬度')
    ax_per.set_aspect('equal')

    plt.show()