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

Python 获取省（或市、县、公司）地理数据。例如，研究气候变化影响的研究应该在结构化的气候中进行。从地理数据（如降水量、湿度、风速）中我们得到了温度反演数据作为计算污染工具变量所需的原始数据：还存储了结构化温度数据。大多数经济学文章关注的是省、市和企业，而地理数据则以结构化数据存储，如栅格数据（.tif 数据格式）和分层数据格式（HDF、.netcdf 数据格式）（还有位置数据，此处不涉及，我们正在讨论）。对于刚接触经济学的人来说，将省份等经济数据与地理数据合并是一个大问题。这条推文分享了如何处理这个问题，既是为了阐明初学者的想法，也是为了阐述他自己的观点。

对于省（市）数据，气候数据合并的一种方法是使用掩模法将省（市）数据所在区域的矢量数据与气候数据合并，然后进行合并。将所得数据在范围内取平均值，即可得到该范围（或城市）对应的地理数据。这种方法背后的想法相对简单。但需要注意的是，气候数据的空间分辨率往往较高，对于面积较小的省（市），最终结果可能为零。因此，为了避免这一问题，通常采用空间插值方法（如反距离加权）将高空间分辨率数据插值到低空间分辨率数据中，然后将其与代表省（或市）的矢量数据相结合。合并并平均。但是，这样计算了两次数据，所以数据的偏差会太小。为了避免这个问题，您可以使用下面的另一种方法。

另一种方法是使用省（市）行政中心的气象数据作为省（市）的气象数据（李，2021）。这篇推文解释了如何使用Python从一个省（或市）的行政中心获取气象数据。这条推文的结构包括以下几个方面。首先是获取省（市）行政中心的WGS84坐标。二是根据行政中心WGS84坐标获取相应位置的气象数据。

获取行政中心WGS84坐标

目前获取特定位置的地理坐标有两种方式（即使用两种地理编码方法），一种是使用百度API，另一种是使用高德地图API。但要注意的是，这两种方法得到的坐标都不是WGS84坐标。高德API获取GCJ02和百度API BD09坐标。这两个坐标都经过加密，在一定程度上抵消了WGS84坐标，所以我们不能直接使用这两个坐标。因此，以高德为例，使用API​​时有两个步骤。第一步是使用高德获取GCJ02坐标，第二步是将GCJ02坐标转换为WGS84坐标。下面介绍这两个步骤。

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

获取省级行政中心坐标非常简单，只需在高德中输入省份名称即可（稍后查看）。

代码如下：

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

首先打开一个包含要分析的省份或城市名称的数据框。如果包含城市，可以先添加省份和城市列，例如“江苏省连云港市”就是代表省份的字符串。

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

其次，将数据框名称和代表省、市、县、公司名称的字符串列名称传递给parse函数，即可通过Amap API获取经纬度。

dfz = parse(prov_name,"prov_name")
dfz

最后使用索引位置合并原始数据 得到的数据将被合并

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

导出为csv

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

将GCJ02坐标转换为WGS84坐标

可以使用完成的python代码将GCJ02坐标转换为WGS84 坐标。使用方法如下：

首先，从https://github.com/wandergis/coordTransform_py下载Python代码并将其解压到您指定的位置。我将其解压到文件夹C:\Users\10197\coordTransform_py。

再次调用该函数。其中，待处理的csv的路径中

后面的参数写在上面，这里传递的是 prov_name_loc_gcj.csv的相对路径。 -o更改为csv路径中的WGS84坐标，即prov_name_loc_wgs.csv。 -t g2w后面的参数表示GCJ02转换为WGS84。 -n 输入长度名称。 -a 输入宽度名称。

!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

最后，打开转换后的数据。

检查是否是行政中心。可以发现，数字121.473701和31.230416对应的地点是上海市黄浦区人民大道，上海黄浦区人民大道是上海市人民政府所在地。因此，可以直接代表省或市搜索，得到行政中心的位置。

提取县城中心的平均风速

首先读取netcdf格式的数据，以温度数据为例，注意路线不能包含中文名称。

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坐标。

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

最后使用以下代码获取坐标位置的温度值。

'''
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