栅格区域人口分布数据获取及坐标系转换

最新推荐文章于 2024-06-05 09:15:23 发布

windawdaysss

最新推荐文章于 2024-06-05 09:15:23 发布

阅读量4.8k

点赞数 3

本文链接：https://blog.csdn.net/weixin_40356612/article/details/127570668

版权

Python 同时被 3 个专栏收录

22 篇文章 0 订阅

订阅专栏

天气雷达

3 篇文章 0 订阅

订阅专栏

MATLAB

3 篇文章 1 订阅

订阅专栏

前言

需要获取的是目标栅格区域内的人口分布（密度）数据。本文从数据获取，到处理方式上一步步详细进行阐述，借助工具有：arcgis10.7，python3.7，matlabR2018b。

一、人口分布数据下载和格式转换

资源环境科学与数据中心，网址: https://www.resdc.cn/Default.aspx
依次选择左侧的数据集（库）目录----->>社会经济数据----->>中国人口空间分布公里网格数据集，如下图：
在这里插入图片描述
下载数据需要注册和登录，拉到网页最下端，有如下下载列表，选择2019年数据进行下载。

下载下来的数据是.adf格式数据，需要转tif，shp，或netcdf格式数据的同学，请参照https://blog.csdn.net/rqjabc/article/details/124771567
该数据为栅格数据类型，每个栅格代表该网格范围（1平方公里）内的人口数，单位为人/平方公里，数据格式为gird，数据以Krassovsky椭球为基准，投影方式为Albers投影。

二、坐标系统转换

所需要的是基于经纬度的地理坐标系，然而，现在拿到的数据是投影坐标系，具体是krasovsky_1940_Albers到GCS_WGS_1984 关于地理坐标系统和投影坐标系统想了解的可以查看https://blog.csdn.net/weixin_43465015/article/details/110875759，说白了，投影坐标系=地理坐标系+投影方法，两种坐标系统的转换使用arcgis工具来实现。

首先，需要定义投影（如果投影坐标系和要转换的地理坐标系是同一个地理坐标系，则直接调至下一步），打开arcgis，系统工具箱------->>数据管理工具------->>投影和变换------->>创建自定义地理（坐标）变换。如下图：

地理（坐标）变换名称自己起，方便识别就行，输入地理坐标系和输出地理坐标系在右侧的小手图标打开，可通过搜索选择。
其次，栅格投影，系统工具箱------->>数据管理工具------->>投影和变换------->>栅格------->>投影栅格，如下图：

输入栅格为需要转换的数据，输出数据集我这里选择的默认地址，输出坐标系和上面一样，此时地理坐标变换会自动识别出上面我们定义的投影，确定后即可完成转换。

三、数据解析和插值

我这里转换为tif格式数据，以下为该格式数据解析的python代码，需要gdal包的加持。

from osgeo import gdal
import numpy as np

def read_data(data_path):
    dataset = gdal.Open(data_path)  # 打开tif
    projection = dataset.GetProjection()  # 地理/投影坐标系

    im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵，含有 6 个元素的元组
    # [左上角的x坐标, 像素宽度, 行旋转（通常为零）, 左上角的y坐标, 列旋转（通常为零）, 像素高度(北半球上图像为负值)]。

    im_width = dataset.RasterXSize  # 获取宽度，数组第二维，左右方向元素长度，代表经度范围
    im_height = dataset.RasterYSize  # 获取高度，数组第一维，上下方向元素长度，代表纬度范围
    band = dataset.RasterCount  # 波段数
    im_data = dataset.GetRasterBand(1).ReadAsArray(xoff=0, yoff=0, win_xsize=im_width, win_ysize=im_height)  # 数据

    im_lon = [im_geotrans[0] + i * im_geotrans[1] for i in range(im_width)]  # 经度列表
    im_lat = [im_geotrans[3] + i * im_geotrans[5] for i in range(im_height)]  # 纬度列表

    return np.array(im_lon), np.array(im_lat), im_data

由于想要更精细的数据，因此需要对数据进行插值处理，这里应用的是matlab中的interp2函数进行二维插值，matlab代码如下：

[lon_data1, lat_data1] = meshgrid(lon_data, lat_data);
[lon_u1, lat_u1] = meshgrid(lon_u, lat_u);
inter_data= interp2(lon_data1, lat_data1, double(origin_data), lon_u1, lat_u1, 'spline'); %三次样条插值
% 'linear' ：双线性插值算法(缺省算法);
% 'nearest' ：最临近插值;
% 'spline' ：三次样条插值;
% 'cubic' ：双三次插值。
save('C:\Users\ADMIN\PycharmProjects\pythonProject\inter_data.mat', 'inter_data');