GDAL: Clip Rasters

Introduction

GDAL利用shapefile文件裁剪栅格影像，并配置裁剪边界之外的背景数值。注意shapefile文件与栅格影像具有完全一致的投影信息。

Example

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

% Created by LI Xu % Version 1.0 % April 19, 2016 % Description: % Clip rasters from a shapefile boundary % If you have any question about this code, % please do not hesitate to contact me via E-mail: % jeremy456@163.com % Blog: % http://blog.sciencenet.cn/u/lixujeremy % http://lixuworld.blogspot.com/ clear; clc; % Source Directory SouDir='input'; % Destination Directory DesDir='output'; % All Input Files files=dir(fullfile(SouDir, '*.tif')); % Vector File vecpath='shapefile.shp'; [~, layername, ~]=fileparts(vecpath); % Output No Data NoData=-2; % Loop for ii=1:length(files) filename=files(ii).name; filepath=fullfile(SouDir, filename); otpath=fullfile(DesDir, filename); % gdalwarp -of GTiff -cutline INPUT.shp -cl shapefile -crop_to_cutline INPUT.tif OUTPUT.tif strcmd=['gdalwarp -of GTiff -cutline ', vecpath, ' -cl ', layername, ' -crop_to_cutline ']; strcmd=[strcmd, '-dstnodata ' num2str(NoData),' ', filepath, ' ', otpath]; py.os.system(strcmd); disp([num2str(ii), ':', filename]); end disp('************************************************');

References

[1] gdalwarp home page.

[2] gdalwarp changing value in nodata pixels.

Matlab: Call User-defined Python Functions

Introduction

Matlab可以调用Python自定义函数，但似乎这个函数还不太好用。这里介绍两种调用方式供参考。

Example

自定义的Python函数：

1 2 3 4 5 6 7 8 9

import sys def squared(x): y = x * x return y if __name__ == '__main__': x = float(sys.argv[1]) sys.stdout.write(str(squared(x)))

第一种方式，应用python（这是M函数）函数，该函数来源请见参考文献[2]、[3]。

1 2 3 4 5 6

% Method One result=python('sqd.py', '5'); disp('Method One'); fprintf('Variable result returns %s ', result); fprintf('with the type of %s.\n', class(result)); disp('-------------------------------------');

第二种方式，应用system函数。

1 2 3 4 5

% Method Two disp('Method Two'); [~, cmdout]=system('python sqd.py 5'); fprintf('Variable cmdout returns %s ', cmdout); fprintf('with the type of %s.\n', class(cmdout));

以上两种方式返回结果完全一致。Matlab调用方法选择任意一种方式，请仔细核对参数传递过程中的对应关系，正确对应才会返回理想结果。File.rar.

References

[1] Call Python Libraries.

[2] Call Python function from MATLAB? Answered by Adrian, kyphos.

[3] python函数参考Adrian的回复.

Code: Archiver Code

[1] 人口压力密度分布.

[2] GDP密度分布.

[3] MODIS NDVI切换为AVHRR NDVI代码.

[4] 牲畜密度分布.

[5] 典型草原1989~2014年NDVI数据.

[6] 典型草原日照时数的处理代码。

[7] 二级土地利用图转一级土地利用图处理代码。

[8] 土地利用变化检测.

[9] 土地利用类型转换分析.

[10] 草地覆盖度空间分布.

[11] 典型草原1989~2014年生长季NDVI数据.

博士论文数据

[1] 指标数据.

[2] 检测Geotiff格式数据数值区间，单波段或多波段，加强版【2023-08】.

[3] 土地利用类型结合行政区划提取数值.

[4] 年度评价代码及归纳评价代码.

[5] 年度评价之全局评价靶心度贡献率代码（纵向评价靶心度贡献率代码）.

[6] 拆分多波段文件.

[7] 矩阵分级代码.

[8] 等级面积比重代码.（面向归纳评价）

[9] 提取生态安全指数（年度、归纳）.

[10] 归纳评价生态安全指数贡献率（全局、纵向）.

[11] 归纳权重识别（全局、纵向）.

[12] 重心识别（年度、归纳）.

[13] 指标分量的线性趋势及其显著性.

[14] 分组，返回索引位置，不同分组依次位于元胞之中.

[15] 归纳评价得分之子系统分量占比（全局、纵向）.

[16] 归纳评价之指标分量一元二次曲线变化率的变化率.

[17] 年度评价之全局评价公式修正.

[18] 指标评估得分.

25指标

[1] 删除部分靶心系数.

[2] 年度指标分量贡献提取（全局、纵向）.

[3] 栅格年度指标分量（全局、纵向）.

[4] 栅格多年指标分量（全局、纵向）.

[5] 多年指标分量贡献提取（全局、纵向）.

[6] 逐年准则层指标分量贡献（全局、纵向）.

[7] 多年准则层指标分量贡献（全局、纵向）.

[8] 重命名文件，达到重新排序的目的.

[9] 多年全局视角的权重.

[10] 添加指标代码.

[11] 趋势百分比.

[12] 逐年栅格线性回归参数a和b，显著性是1.

[13] 多年栅格一元二次回归参数a/b/c.

[14] 逐年趋势预测情景.

[15] 修剪情景数值分布在0~1之间.

[16] 矩阵分级（6等级）.

[17] 等级比重6等级（逐年、多年）.

[18] 一元二次情景计算ABC.

[19] 识别起止年份.

[20] 移动平均求算多年合成归纳数值.

[21] 移动平均三子系统指数分量（全局、纵向）.

[22] 移动平均指标层指数分量（全局、纵向）.

[23] 识别（移动平均）三子系统分地区指数分量.

[24] 识别（移动平均）指标层分地区指数分量.

[23] 移动平均时序演化.

[24] 逐年经济演化情景.

[25] 移动平均经济演化情景.

[26] 逐年情景渲染Package.

[27] 移动平均情景渲染Package.

[28] 逐年多情景生态安全数值分区提取.

[29] 移动平均多情景生态安全数值分区提取.

[30] 敏感度系数（逐年、移动平均）.

[31] 移动平均的指标评估得分.

[32] 时序规律情景计算.

[33] 情景抽取态势指数.

[34] 情景抽取等级比重.

[35] 情景渲染和构图成片.

[36] 斜率截距渲染和构图成片.

Matlab+Python: Theil–Sen estimator

Summary

In non-parametric statistics, there is a method for robust simple linear regression that chooses the median slope among all lines through pairs of two-dimensional sample points. It has been called the Theil–Sen estimator, Sen's slope estimator, slope selection, the single median method, the Kendall robust line-fit method, and the Kendall–Theil robust line. It is named after Henri Theil and Pranab K. Sen, who published papers on this method in 1950 and 1968 respectively.

Example

本例以Anscombe's Quartet为例，检验Theil–Sen estimator对趋势的估计效果。关于Anscombe's Quartet四组数据的特征请参考这里的描述。这里的代码整合Matlab+Python，Theil–Sen estimator计算调用Python Scipy。

Fig. 1

观察Fig. 1，蓝线是简单线性回归，红线是Theil–Sen estimator，两条灰虚线之间表示α=0.9水平下的置信区间。A与D，蓝线与红线基本重合，A的蓝线基本可以反映数据的趋势，但D的蓝线显然不能反映数据分布的趋势，此处的红线也还是不能很好的剔除异常值的影响。B组数据并非呈现线性趋势，所以红线并不会比蓝线更好。C组红线非常好地剔除一个异常值的影响，很好地诠释其他10个数据分布趋势。综上，Theil–Sen estimator的确比简单线性回归优秀，但是也同样存在着无效的情况。

以下Fig. 2~5是不同α水平下的置信区间范围，很显然，随着α增大区间范围增大，更多的数据进入置信区间。

Fig. 2

Fig. 3

Fig. 4

Fig. 5

Fig. 6

MatThie.rar.

栅格处理代码.

多波段栅格处理代码.

References

[1] scipy.stats.mstats.theilslopes.

[2] Theil–Sen estimator, from Wikipedia, the free encyclopedia.

Python+GDAL: 文件格式的转换

空间数据的格式千变万化，常见的有ENVI Standard、GEOTIFF、Erads Image等等（栅格文件列表，矢量文件列表），以下应用GDAL工具转换文件的格式，省时高效。

举个栗子，ENVI Standard转换为GEOTIFF，参考How to call gdal_translate from Python code? Answered by Max。

Note

os.system (command): Execute the command (a string) in a subshell. This is implemented by calling the Standard C function system(). On Windows, the return value is that returned by the system shell after running command. 实际上，示例是Python调用C+GDAL的代码。

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

##Created by LI Xu ##Version 1.0 ##7 October, 2014 ##Convert ENVI Standard to GEOTIFF ##http://gis.stackexchange.com/questions/42584/how-to-call-gdal-translate-from-python-codearray-returns-just-nan-values-when-trying-to-read-envi-file import gdal import ogr import os import datetime import time def IsSubString(SubStrList,Str): flag=True for substr in SubStrList: if not(substr in Str): flag=False return flag #~ #---------------------------------------------------------------------- def GetFileList(FindPath,FlagStr=[]): import os FileList=[] FileNames=os.listdir(FindPath) if (len(FileNames)>0): for fn in FileNames: if (len(FlagStr)>0): if (IsSubString(FlagStr,fn)): fullfilename=os.path.join(FindPath,fn) FileList.append(fullfilename) else: fullfilename=os.path.join(FindPath,fn) FileList.append(fullfilename) if (len(FileList)>0): FileList.sort() return FileList def convert_format(infile, otfile, oFormat): os.system("gdal_translate -of" +" " +oFormat+" "+ infile + " " + otfile) ##Main begintime=time.strftime("%Y-%m-%d %H:%M:%S") print 'Time Begin:'+begintime starttime = datetime.datetime.now() #Source Directory SouDir=r'E:\Tools\rad\NEW' #Destination Directory DesDir=r'E:\Tools\rad\Rad' #Retrieve files FlagStr='.hdr' files=GetFileList(SouDir,FlagStr) for file in files: #print file str_file=file.split('.') infile=str_file[0]+'.tif' #print infile oformat="GTiff" ##Output Image filename=infile.split('\\') filename=filename[len(filename)-1] filename=filename.split('.') otimg=DesDir+'\\'+filename[0]+'.tif' print otimg convert_format(infile, otimg, oformat) print "END" endtime=time.strftime("%Y-%m-%d %H:%M:%S") print 'Time End:'+endtime finishtime=datetime.datetime.now() timespan=(finishtime-starttime).seconds timespan='%f' %timespan print 'Time Span:'+timespan+' s'

Attachment.

Matlab+GDAL版本.