Python获取省（或市、县、企业）地理数据

引言

如今，将地理数据变量作为经济学论文中的变量已是家常便饭，如研究气候变化影响的论文需要在结构化的气候地理数据（如降水、湿度、风速）中得到，计算作为污染的工具变量的逆温数据需要的原始数据也保存在结构化的温度数据中。大部分经济学论文以省份、城市、企业为研究对象，而地理数据则存放在结构化数据，如栅格数据（.tif数据格式）和分层数据格式（Hierarchical Data Formats，HDF，.netcdf数据格式）（还有一种是站点数据，我们在此不讨论这种情况）。如何将省份等经济数据与地理数据合并，对于初学的经济学的人来说，是一个很大的问题。本推文分享如何处理这一问题，一来帮助初学者厘清思路，二来供自己复习之用。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

对于省份（或城市）数据，将气候数据与其合并的一种方法是使用掩膜的方法将省份（或城市）数据所在区域的矢量数据与气候数据进行合并，并将合并得到的数据在省份内取均值来得到省份（或城市）对应的地理数据。这种方法的思想较为简单。然而，需要注意的是，气候数据的空间分辨率往往较大，对于面积较小的省份（或城市），可能最后得到的结果为空值。因此，为了避免这种问题，往往会先采取空间插值的方法（如逆距离加权）将高空间分辨率插值为低空间分辨率的数据，然后再与表示省份（或城市）区域的矢量数据进行合并并取均值。然而，这对数据进行了两次计算，导致数据的variation过小。为了避免这种问题，可以使用下面的另一种方法。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

另一种方法是使用省份（或城市）的行政中心的气象数据作为省份（或城市）的气象数据（Li, 2021）。本推文讲解如何使用Python获得省份（或城市）行政中心的气象数据。本推文的结构包含以下方面，首先是获取省份（或城市）的行政中心的WGS84坐标，其次是根据行政中心的WGS84坐标来获得对应位置的气象数据。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

获得行政中心的WGS84坐标

当前有两种获得某个地点地理坐标的方式（即两种地理编码方式），一种是使用百度的API，另一种是使用高德的API。但是我们需要注意的是，这两种方式得到的坐标都不是WGS84坐标，高德API获得的是GCJ02坐标，百度API获得的是BD09坐标。这两种坐标都是通过加密的方式来使WGS84坐标进行一定程度偏移，因此我们不能直接使用这两种坐标。因此，以高德为例，在使用API时，分为两步，第一步是使用高德获得GCJ02坐标，第二步是将GCJ02坐标转换为WGS84坐标。下面演示这两个步骤的做法。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

使用高德API获得行政中心的GCJ02坐标

获得省份行政中心的坐标很简单，在高德中输入省份名称即可（将在后面进行验证）。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

代码如下：文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

import requests
import pandas as pd
import json

def parse(prov_namedf,prov_namecol):
 '''
 prov_namedf: 要解析的包含省份或城市名称的数据框
 prov_namecol: 省份名称所在列
 '''
    prov_nameList = prov_namedf[prov_namecol].to_list()
    lists = [prov_nameList[i:i + 10] for i in range(0, len(prov_nameList), 10)]
    strlist = ["|".join(lists[i])  for i in range(0, len(lists))]

    dfz = pd.DataFrame()
    ak = 'c8f79e53c469a99dfb4b8dbf'  # your key
    for cityName in strlist:
        base = "http://restapi.amap.com/v3/geocode/geo?key={}&address={}&batch=true" .format(ak, cityName)
        response = requests.get(base)
        answer = response. json()
        # print(answer)
        dff = pd.DataFrame(answer['geocodes'])[['formatted_address','location','adcode']]
        dff['lon'] = dff['location'].apply(lambda x :x.split(",")[0])
        dff['lat'] = dff['location'].apply(lambda x :x.split(",")[1])

        dfz = pd.concat([dfz,dff])
        dfz = dfz.reset_index(drop=True)
    return dfz

首先打开一个要解析的包含省份或城市名称的数据框，如果是包含城市的，可以先将省份和城市的列加起来，形成如“江苏省连云港市”这样的表示省份的字符串。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

prov_name = pd.read_excel("d:\Desktop\公众号\地理编码\省份名称.xlsx")
prov_name

其次，将数据框名称和表示省或市或县或公司名称的字符串的列名传递给parse函数，得到使用高德API得到的经纬度文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

dfz = parse(prov_name,"prov_name")
dfz

最后，使用索引位置将原数据与得到的数据合并文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

prov_name_loc = pd.concat([prov_name,dfz],axis=1)  ## 不要变换prov_name 的index位置
prov_name_loc

导出为csv文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

prov_name_loc.to_csv(r"prov_name_loc_gcj.csv",index=False)

转换GCJ02坐标为WGS84坐标

使用现成的python代码，可以将GCJ02坐标转换为WGS84坐标。使用方法如下：文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

首先，在https://github.com/wandergis/coordTransform_py中下载python代码，并解压到自己指定的位置，我这里解压到C:\Users\10197\coordTransform_py。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

其次，调用函数。其中，-i后面的参数传入上面写入的要处理的csv的路径，这里传入相对路径 prov_name_loc_gcj.csv。-o后面的参数传入改变为WGS84坐标的csv的路径，即prov_name_loc_wgs.csv。-t后面的参数 g2w表示GCJ02转为WGS84。-n传入经度名称。-a 传入纬度名称。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

!python C:\Users\10197\coordTransform_py\coord_converter.py -i prov_name_loc_gcj.csv -o prov_name_loc_wgs.csv -t g2w -n lat -a lon

最后，打开转换完成的数据。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

验证一下，是否是行政中心。可以发现，121.473701,31.230416对应的地点为上海市黄浦区人民大道，而上海市黄浦区人民大道就是上海市人民政府所在地。因此，可以直接搜索省份或城市的名称得到行政中心位置。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

提取县中心的平均风速

首先，读取netcdf格式的数据，以温度数据为例，注意路径中不能有中文名。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

import xarray
from pyproj import Transformer

tem = xarray.open_dataset(r"d:\Desktop\example.nc",mask_and_scale=True) \
    .sel(time = slice("1996-01-01","2022-12-31"))['t2m'].rio.write_crs("4326")  ## 切记路径中不能有中文名
tem

接着，获得WGS84坐标。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

lonlt,latlt=prov_name_loc2['lon'], prov_name_loc2['lat']
lonlt,latlt

最后，使用下面的代码获得坐标位置的温度值。文章源自菜鸟学院-https://www.cainiaoxueyuan.com/bc/54124.html

'''
lonlt: 想要转换的点的经度的列表
latlt: 想要转换的点的维度的列表
raster: 想要提取的netcdf数据或raster数据
raster_lonlat_name: 想要提取的netcdf数据或raster数据的经度和维度坐标的名称的列表
'''
rasterr = tem           ## 改 想要提取的netcdf数据或raster数据
transformer = Transformer.from_crs("EPSG:4326", rasterr.rio.crs, always_xy=True)
xx, yy = transformer.transform(lonlt, latlt)


point_valuelt = []
for lon,lat in zip(xx,yy):
    value = rasterr.sel(longitude=lon, latitude=lat, method="nearest").values     ## 改 rasterr的经纬度坐标名称
    point_valuelt.append(value)

si10_df= pd.DataFrame(point_valuelt,columns = tem['time'].values )   ## 改要命名的df的名称和df的列名  
si10_df