前言
GWR软件是实现地理加权回归建模的专业软件。
GWR软件下载地址Click here。
软件下载后自带了User guidance,也可以参照这篇博客。
软件要求输入的数据格式为.csv。数据中要包括投影坐标,因变量,解释变量。
数据制作
我这里以地表温度(LST)做因变量,NDVI做解释变量为例,使用Arcgis制作.csv输入文件。
LST和NDVI的原始数据是两幅投影坐标的栅格影像。
首先在Arcgis中打开影像,然后使用栅格转点工具将栅格转换为点状shp文件。
得到每个像素点的shp文件。
原先栅格每个像素的值会自动在shp文件属性表中生成一个字段GRID_CODE。
接下来为属性表添加投影坐标数据。
打开shp图层的属性表,添加字段,并计算几何。同理添加Y坐标也这样操作。
成功添加XY坐标后,属性表左上角导出,选择保存类型为文本文件,然后手动将文件后缀改为.csv。
至此得到GWR要求的.csv输入文件。
问题
地理加权回归一般会在省市间进行,元素数量级一般在几十,几百。但是我的需求是对栅格点进行地理加权回归,数量级比较大。输入数据比较大时,GWR软件会报错。所以我考虑使用pysal
这个python库进行GWR。
pysal.model.mgwr.gwr的Github地址
代码如下:
import pysal.lib as ps
from pysal.model.mgwr.gwr import GWR
import numpy as np
import pandas as pd
excel_file = r'F:\Py_project\GWR\LST_excelNo-1.csv'
dataset = pd.read_csv(excel_file).loc[:, ['LST', 'NDVI', 'X', 'Y']]
coords = list(zip(dataset['X'], dataset['Y']))
y = np.array(dataset['LST'], dtype='float32').reshape(-1, 1)
# y = (y - y.mean(axis=0)) / y.std(axis=0)
X = np.array(dataset['NDVI'], dtype='float32').reshape(-1, 1)
# X = (X - X.mean(axis=0)) / X.std(axis=0)
del dataset
model = GWR(coords, y, X, bw=44, fixed=False, kernel='gaussian')
results = model.fit()
params = results.params
报错如:MemoryError: Unable to allocate 35.9 GiB for an array with shape (69420, 69420) and data type float64
。
同样是因为数据量过大,运行内存不足导致的。
试过修改源码中开辟nunpy数组的数据类型为float16,没有成功。
我看文献中的研究区数据量应该比我这大得多,可能是工作站跑的。也有明确说是使用GWR软件的,这就很疑惑了- _ - ,GWR软件不支持大数据量。