Code前端首页关于Code前端联系我们

Python获取省(或市、县、企业)地理数据

terry 2年前 (2023-09-24) 阅读数 133 #后端开发
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
Python获取省(或市、县、企业)地理数据

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

dfz = parse(prov_name,"prov_name")
dfz
Python获取省(或市、县、企业)地理数据

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

prov_name_loc = pd.concat([prov_name,dfz],axis=1)  ## 不要变换prov_name 的index位置
prov_name_loc
Python获取省(或市、县、企业)地理数据

导出为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_pyPython获取省(或市、县、企业)地理数据

再次调用该函数。其中,待处理的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

最后,打开转换后的数据。 Python获取省(或市、县、企业)地理数据

检查是否是行政中心。可以发现,数字121.473701和31.230416对应的地点是上海市黄浦区人民大道,上海黄浦区人民大道是上海市人民政府所在地。因此,可以直接代表省或市搜索,得到行政中心的位置。Python获取省(或市、县、企业)地理数据Python获取省(或市、县、企业)地理数据

提取县城中心的平均风速

首先读取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
Python获取省(或市、县、企业)地理数据

然后获取WGS84坐标。

lonlt,latlt=prov_name_loc2['lon'], prov_name_loc2['lat']
lonlt,latlt
Python获取省(或市、县、企业)地理数据

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

'''
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
Python获取省(或市、县、企业)地理数据

版权声明

本文仅代表作者观点,不代表Code前端网立场。
本文系作者Code前端网发表,如需转载,请注明页面地址。

发表评论:

◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。

热门