前言

需要獲取的是目標(biāo)柵格區(qū)域內(nèi)的人口分布（密度）數(shù)據(jù)。本文從數(shù)據(jù)獲取，到處理方式上一步步詳細(xì)進(jìn)行闡述，借助工具有：arcgis10.7，python3.7，matlabR2018b。

一、人口分布數(shù)據(jù)下載和格式轉(zhuǎn)換

資源環(huán)境科學(xué)與數(shù)據(jù)中心，網(wǎng)址: https://www.resdc.cn/Default.aspx
依次選擇左側(cè)的數(shù)據(jù)集（庫）目錄----->>社會(huì)經(jīng)濟(jì)數(shù)據(jù)----->>中國人口空間分布公里網(wǎng)格數(shù)據(jù)集，如下圖：

下載數(shù)據(jù)需要注冊(cè)和登錄，拉到網(wǎng)頁最下端，有如下下載列表，選擇2019年數(shù)據(jù)進(jìn)行下載。

下載下來的數(shù)據(jù)是.adf格式數(shù)據(jù)，需要轉(zhuǎn)tif，shp，或netcdf格式數(shù)據(jù)的同學(xué)，請(qǐng)參照https://blog.csdn.net/rqjabc/article/details/124771567
該數(shù)據(jù)為柵格數(shù)據(jù)類型，每個(gè)柵格代表該網(wǎng)格范圍（1平方公里）內(nèi)的人口數(shù)，單位為人/平方公里，數(shù)據(jù)格式為gird，數(shù)據(jù)以Krassovsky橢球?yàn)榛鶞?zhǔn)，投影方式為Albers投影。

二、坐標(biāo)系統(tǒng)轉(zhuǎn)換

所需要的是基于經(jīng)緯度的地理坐標(biāo)系，然而，現(xiàn)在拿到的數(shù)據(jù)是投影坐標(biāo)系，具體是krasovsky_1940_Albers到GCS_WGS_1984 關(guān)于地理坐標(biāo)系統(tǒng)和投影坐標(biāo)系統(tǒng)想了解的可以查看https://blog.csdn.net/weixin_43465015/article/details/110875759，說白了，投影坐標(biāo)系=地理坐標(biāo)系+投影方法，兩種坐標(biāo)系統(tǒng)的轉(zhuǎn)換使用arcgis工具來實(shí)現(xiàn)。

• 首先，需要定義投影（如果投影坐標(biāo)系和要轉(zhuǎn)換的地理坐標(biāo)系是同一個(gè)地理坐標(biāo)系，則直接調(diào)至下一步），打開arcgis，系統(tǒng)工具箱------->>數(shù)據(jù)管理工具------->>投影和變換------->>創(chuàng)建自定義地理（坐標(biāo)）變換。如下圖：
  
  地理（坐標(biāo)）變換名稱自己起，方便識(shí)別就行，輸入地理坐標(biāo)系和輸出地理坐標(biāo)系在右側(cè)的小手圖標(biāo)打開，可通過搜索選擇。

• 其次，柵格投影，系統(tǒng)工具箱------->>數(shù)據(jù)管理工具------->>投影和變換------->>柵格------->>投影柵格，如下圖：
  
  輸入柵格為需要轉(zhuǎn)換的數(shù)據(jù)，輸出數(shù)據(jù)集我這里選擇的默認(rèn)地址，輸出坐標(biāo)系和上面一樣，此時(shí)地理坐標(biāo)變換會(huì)自動(dòng)識(shí)別出上面我們定義的投影，確定后即可完成轉(zhuǎn)換。

三、數(shù)據(jù)解析和插值

我這里轉(zhuǎn)換為tif格式數(shù)據(jù)，以下為該格式數(shù)據(jù)解析的python代碼，需要gdal包的加持。

from osgeo import gdal import numpy as npdef read_data(data_path):dataset = gdal.Open(data_path) # 打開tifprojection = dataset.GetProjection() # 地理/投影坐標(biāo)系im_geotrans = dataset.GetGeoTransform() # 獲取仿射矩陣，含有 6 個(gè)元素的元組# [左上角的x坐標(biāo), 像素寬度, 行旋轉(zhuǎn)（通常為零）, 左上角的y坐標(biāo), 列旋轉(zhuǎn)（通常為零）, 像素高度(北半球上圖像為負(fù)值)]。im_width = dataset.RasterXSize # 獲取寬度，數(shù)組第二維，左右方向元素長度，代表經(jīng)度范圍im_height = dataset.RasterYSize # 獲取高度，數(shù)組第一維，上下方向元素長度，代表緯度范圍band = dataset.RasterCount # 波段數(shù)im_data = dataset.GetRasterBand(1).ReadAsArray(xoff=0, yoff=0, win_xsize=im_width, win_ysize=im_height) # 數(shù)據(jù)im_lon = [im_geotrans[0] + i * im_geotrans[1] for i in range(im_width)] # 經(jīng)度列表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

由于想要更精細(xì)的數(shù)據(jù)，因此需要對(duì)數(shù)據(jù)進(jìn)行插值處理，這里應(yīng)用的是matlab中的interp2函數(shù)進(jìn)行二維插值，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');

END

參考資料

• https://blog.csdn.net/rqjabc/article/details/124771567
• https://blog.csdn.net/weixin_43465015/article/details/110875759
• https://blog.csdn.net/qq_33657870/article/details/103647090
• https://ww2.mathworks.cn/help/matlab/ref/interp2.html?s_tid=srchtitle_interp2_1

總結(jié)

以上是生活随笔為你收集整理的栅格区域人口分布数据获取及坐标系转换的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。