技術標籤：python進行時空資料處理pythonapi

這是《地圖時空大資料爬取》第6章的內容，這篇部落格主要是抓取一下行政區劃資料，最小是能抓到區縣的行政區劃資料。然後書裡是用arcpy和arcgis再加上python一起來處理的，有些麻煩，我統一用python來處理了，改寫了下程式碼，因為很簡單，就不多說了，直接放程式碼。

import basics
import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import json
import urllib
from
 
 urllib.parse import quote
import string
from shapely.geometry import Point 
from shapely.geometry import Polygon
from shapely.geometry import LineString
import pyproj

# 功能：採集行政區域邊界
# 返回儲存邊界資訊的BoundryWithAttr物件
def getDistrictBoundry(ak,citycode):
    districtBoundryUrl = "http://restapi.amap.com/v3/config/district?key="
 
+ak+"&keywords="+citycode+"&subdistrict=0&extensions=all"
    print(districtBoundryUrl)
    districtBoundryUrl = quote(districtBoundryUrl, safe=string.printable)
    json_obj = urllib.request.urlopen(districtBoundryUrl)
#     json_obj = response.read().decode('utf-8','ignore').replace(u'\xa9', u'')
 

    json_data=json.load(json_obj)
    districts=json_data['districts']
    try:
        polyline = districts[0]['polyline']
        centerLon=districts[0]['center'].split(',')[0]
        centerLat=districts[0]['center'].split(',')[1]
    except Exception as e:
        print("error!")
    pointscoords=polyline.split(';')
    point=basics.PointWithAttr(0,centerLon,centerLat,"行政區域",citycode)
    districtBoundry=basics.BoundryWithAttr(point,pointscoords)
    return districtBoundry

ak = "放上你的ak就行"
citycodes={'上城區':'330102','下城區':'330103','江乾區':'330104',"拱墅區":'330105',"西湖區":'330106','濱江區':'330108',
           '蕭山區':'330109','餘杭區':'330110','富陽區':'330111','臨安區':'330112','桐廬縣':'330122','淳安縣':'330127',
          '建德市':'330182'}
polygon_list = []
centerX_list = []
centerY_list = []
for citycode in citycodes.keys():
    districtBoundry = getDistrictBoundry(ak,citycode)
    tmp = districtBoundry.boundrycoords
    for i in range(len(tmp)):
        tmp[i] = [float(tmp[i].split(',')[0]) , float(tmp[i].split(',')[1])]
    tmp = Polygon(tmp)
    polygon_list.append(tmp)
    centerX_list.append(districtBoundry.point.lon)
    centerY_list.append(districtBoundry.point.lat)

# 將爬取得到的杭州市行政區儲存為shp檔案
gdf = {'geometry':polygon_list , 'centerLon':centerX_list , 'centerLat':centerY_list , 
       'name':citycodes.keys(),'adacode':citycodes.values()}
gdf = gpd.GeoDataFrame(gdf , crs = None)
gdf

gdf.plot()

gdf.to_file('.\\杭州市行政規劃\\杭州市行政區劃.shp', encoding="utf-8")

有一個問題無法解決，網上有人遇到這個問題，但是沒有給出解決方案，一旦我把投影更改成epsg:4326之後，就無法成功將gdf寫入到shp or geojson等檔案中去了。若有大佬解決該問題還請教一下解決方案，謝謝。

print(gdf.crs)

None

gdf.crs = 'epsg:4326'
gdf.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

gdf.to_file('.\\杭州市行政規劃\\杭州市行政區劃.shp', encoding="utf-8")

---------------------------------------------------------------------------

CRSError                                  Traceback (most recent call last)

<ipython-input-92-231f1b2525d8> in <module>
----> 1 gdf.to_file('.\\杭州市行政規劃\\杭州市行政區劃.shp', encoding="utf-8")


E:\Anaconda\lib\site-packages\geopandas\geodataframe.py in to_file(self, filename, driver, schema, index, **kwargs)
    744         from geopandas.io.file import _to_file
    745 
--> 746         _to_file(self, filename, driver, schema, index, **kwargs)
    747 
    748     def set_crs(self, crs=None, epsg=None, inplace=False, allow_override=False):


E:\Anaconda\lib\site-packages\geopandas\io\file.py in _to_file(df, filename, driver, schema, index, mode, crs, **kwargs)
    253             crs_wkt = crs.to_wkt("WKT1_GDAL")
    254         with fiona.open(
--> 255             filename, mode=mode, driver=driver, crs_wkt=crs_wkt, schema=schema, **kwargs
    256         ) as colxn:
    257             colxn.writerecords(df.iterfeatures())


E:\Anaconda\lib\site-packages\fiona\env.py in wrapper(*args, **kwargs)
    398     def wrapper(*args, **kwargs):
    399         if local._env:
--> 400             return f(*args, **kwargs)
    401         else:
    402             if isinstance(args[0], str):


E:\Anaconda\lib\site-packages\fiona\__init__.py in open(fp, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt, **kwargs)
    275             c = Collection(path, mode, crs=crs, driver=driver, schema=this_schema,
    276                            encoding=encoding, layer=layer, enabled_drivers=enabled_drivers, crs_wkt=crs_wkt,
--> 277                            **kwargs)
    278         else:
    279             raise ValueError(


E:\Anaconda\lib\site-packages\fiona\collection.py in __init__(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, ignore_fields, ignore_geometry, **kwargs)
    153             self._check_schema_driver_support()
    154             if crs_wkt or crs:
--> 155                 self._crs_wkt = crs_to_wkt(crs_wkt or crs)
    156 
    157         self._driver = driver


fiona\_crs.pyx in fiona._crs.crs_to_wkt()


CRSError: Invalid input to create CRS: GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["unknown"],AREA["World"],BBOX[-90,-180,90,180]],ID["EPSG",4326]]