栅格数据处理

最新推荐文章于 2024-06-21 12:07:56 发布

裘千仞不会水上漂

最新推荐文章于 2024-06-21 12:07:56 发布

阅读量3.8k

点赞数 3

分类专栏：技术交流之GIS 技术交流之python 文章标签：栅格数据 modis rasterio

本文链接：https://blog.csdn.net/xtfge0915/article/details/86201619

版权

技术交流之python 同时被 2 个专栏收录

19 篇文章 13 订阅

订阅专栏

技术交流之GIS

5 篇文章 1 订阅

订阅专栏

栅格数据处理工具及方法推荐

1.前言

最近受好友之邀写一些地理数据处理的内容，虽然觉得水平有限不想误人子弟，奈何抵不住好友的柔情攻势，只能尽可以把自己觉得有用的东西写下来，希望能帮到需要的人。因为矢量数据
和栅格数据存储
方式差异较大，用到的处理工具也有所不同，因此笔者打算分两篇写，这篇文章的主要内容是栅格数据的处理。

2.软件安装和环境配置

鉴于windows系统较多，因此本教程以wiondows系统为例展开说明，在其它系统中操作有任何疑问都可以联系作者。

2.1 andconda安装

ps : 虽然这是一个初级教程，但也需要一定的基础知识储备，包括python基础、GIS基础、程序设计基础，
当然用到的都是一些很简单的知识，没有也不影响大家学习和使用

Anaconda是一个开源的Python发行版本,使用它的主要原因是它的包管理器非常好用。像下面用到的这些库如果用pip安装能把你折腾到怀疑人生。

下载地址：

安装方式并没有什么特别的，一直点下一步就行了(记住安装目录）。
安装完成后修改condal镜像，因为我发现默认的镜像有很多包都没有。国内比较优秀的镜像有清华镜像、阿里镜像、华为镜像等，我们修改成清华镜像，这个它包含我们要用的所有包，其它镜像我没测试过。具体操作是在命令行依次执行：

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/main/
conda config --set show_channel_urls yes

2.2 安装栅格数据处理库

gdal
rasterio

安装方式：

打开命令行（开始->windows系统->命令提示符）输入conda install gdal rasterio
在这里插入图片描述
看到这个界面输入y回车等待安装完成。
ps:如果提示找不到需要的包请检查前面镜像修改是否成功。

3.数据准备

说到数据准备再给大家安利一个指批量下载landsat和modis影像的工具。

3.1 landsat数据批量下载工具

（1）安装landsat-util

在命令行执行pip install landsat-util
如果出现下载错误请参考这篇文章
在这里插入图片描述
如果用pip安装失败只能从源码安装了，虽然conda很好用，但它的包目前还是没法和pip比（不得不说pip真的很坑，反正我用pip没有安装成功）。

首先，从github下载源码，需要用到git工具，git安装命令conda install git，安装完成用在命令行执行git clone https://github.com/developmentseed/landsat-util.git下载landsat-util源码。
其次，编译安装，依次执行以下命令：
```
cd landsat-util
python setup.py install
```
安装完成后在命令行输入landsat -h，如果没有报错那么万事大吉，已经安装成功了，如果提示缺少什么包，那么很遗憾，你还需要继续安装依赖包。
最后，安装依赖，landsat-util依赖的包比较多，需要依次安装，以下是我在操作过程中需要安装的包
1. six:安装命令conda install six
2. usgs:usgs在conda和pip都找不到，需要从源码安装，具体如下：
```
git clone https://github.com/kapadia/usgs.git
cd usgs
python setup.py install
```
3. future:pip install future
4. humanize：conda install humanize
  如何知道缺少什么依赖包呢，那就是不停的执行landsat -h它会提示你缺少哪个包，比如我执行landsat -h后出现
  
  错误信息提示缺少usgs那么你就要安装usgs，安装后继续执行landsat -h直到没有错误。安装一个包的方式较多，优先考虑使用conda，如果conda找不到使用pip，如果pip也找不到或安装出错，使用源码安装。
  
  ps:包比较多，安装起来还是比较麻烦的，如果大家能用pip安装最好，它会自动安装所有依赖。

（2）使用landsat-util检索数据
常用的有两种方式，行列号和经纬度

按行列号搜索：landsat search --cloud 4 --start "january 1 2014" --end "january 10 2014" -p 009,045
- –cloud 4指动量小于4%的影像，也可以写成-c 4，如果不指定默认20%。
- –start "january 1 2014"指搜索的起始时间,–end "january 10 2014"指结束时间
- -p 009,045指pathwrow，前面的是path后面的是row，关于遥感数据path,row有疑问的请自行查阅相关资料。
- 更多参数请用landsat -h查看。
按经纬度搜索：landsat search --lat 38.9004204 --lon -77.0237117
搜索结果是一个json字符串，我们需要的是其中的scenID，下载的时候要用到它。

（3）landsat-util下载数据
```
landsat download LC80090452014008LGN00// LC80090452014008LGN00指的是scenID
landsat download LC80090452014008LGN00 --bands 432//只下载432波段
landsat download LC80090452014008LGN00 LC80090452015008LGN00    LC80090452013008LGN00//下载多个影像
```
ps:scenID并非只能用landsat搜索到，也可以USGS网站检索想到的数据，然后复制scenID用landsat下载，landsat是在google镜像下载数据，速度比USGS快。

3.2 modis数据批量下载工具

（1）安装pymodis：

下载git clone https://github.com/lucadelu/pyModis.git
进行pymodis目录:cd pyModis
安装:python setup.py install
或着pip install pyModis

(2)使用pymodis下载数据
pymodis包含 modis_download.py, modis_download_from_list.py,modis_parse.py,modis_multiparse.py, modis_mosaic.py, modis_convert.py,modis_quality.py共7个模块。其作用看名称就可以知道，我就不赘述了。

#下载数据代码示例
from pymodis import downmodis
# destination foldert
dest = "/tmp"
# tiles to download
tiles = "h18v04,h19v04"#条带号
# starting day
day = "2014-08-14"
# number of day to download
delta = 1
user='xtfge'
passwd='ABCabc123'#这是我刚注册的账号，大家可以随便用。
modisDown = downmodis.downModis(destinationFolder=dest, tiles=tiles, today=day, delta=delta,user=user,password=passwd)
modisDown.connect()
modisDown.downloadsAllDay()

或着也可以在命令行下载

modis_download.py -I -r -t h18v03,h18v04 -f 2008-01-01 -e 2008-01-31 lst_terra

更多内容请在这里下载教程或访问官网。
ps: 不管是下载landsat还是modis数据，可用的工具很多这里抛砖引玉只举了两个例子，想了解更多请自行查阅相关资料。

3.3 下载数据

本章的内容是数据准备，叽叽歪歪说了半天我们还一幅影像都没下载呢，比如我们要下载2018年北京市的landsat8影像，按照这个要求我们肯定能搜索到一堆数据，所以对数据的要求再做细化：同一地方只下载云量最少的一景，那么这个代码应该怎么写呢，最笨的方法当然是全部下载下来慢慢挑（手动滑稽），对于以上要求，设计思路应该是这样的：

确定要下载影像的条带号，最笨但最简单的方法是直接查。

从图中可以看到，北京的条带号为path=[130,131,132,133,134]，row=[30,31,32,33,34,34]
搜索2018年北京市所有影像
官网说是可以用landsat search --cloud 4 --start "january 1 2018" --end "December 31 2018" -p 130,30,130,31这样的形式同时搜索多个条带，但我试了下好像不行。因此就只能一个一搜了。别看path,row只有这么几个，但组合起来可就多了因此写个简单和程序生成搜索代码：
```
rst='landsat search  --start "january 1 2018" --end "December 31 2018" -p'
path=[130,131,132,133,134]
row=[30,31,32,33,34,34]
fp=open('E://search_code.bat','a')
for p in path:
   for r in row:
       fp.write("%s %d,%d > E://%d%d.txt\n" % (rst,p,r,p,r))
fp.write("@echo off\necho 搜索完成，按任意键退出\npause>nul（>nul)\nexit )")
fp.close()
```
如果想看刚才生成的代码，在命令行输入notepad E://search_code.bat，或着直接用记事本打开该文件。
ps:因为生成的是一个.bat文件，这段代码的执行结果可能被360等杀毒软件拦截，请选择不处理或信任。

大家可能对就段代码有疑问，对照结果解释一下它的作用

结果的每一行对应一条搜索语句，可能引起疑惑的最后面的类似 >E://13030.tx这样t的代码，它的作用是将搜索结果保存到E://13030.txt这个文件中，至于 fp.write("@echo off\necho 搜索完成，按任意键退出\npause>nul（>nul)\nexit )")这句与搜索无关，大家就不要深究了，不求甚解也是必要的。
找到刚才生成的search_code.bat文件，右击选择以管理员身份运行。

找到云量最少的影像
现在在我的E盘有一大堆txt文件，它们存储了北京市2018年所有的landsat影像，我们同一坐标只需要云量最少的一景。其本质就是在每一个txt文件中找到cloud最小的那个影像的sceneID

import glob,re,json
filepath='E:'
ls=glob.glob('%s/*.txt' % filepath)#遍历E盘所有txt文件，建议把文件保存到一个新建的目录中，避免遍历到其它文件
r=re.compile('.*?({.*}).*?',re.DOTALL)#python正则表达式，感兴趣可以查阅相关资料
IDs=[]
for l in ls:
   min_cloud=100
   min_cloud_obj=None
   with open(l) as fp:
       rst=json.loads(r.match(fp.read()).group(1))
       for result in rst['results']:
           if result['cloud']<min_cloud:
               min_cloud=result['cloud']
               min_cloud_obj=result
   IDs.append(min_cloud_obj['sceneID'])
   print(l,min_cloud_obj['cloud'],min_cloud_obj['sceneID'])

下载影像
现在符合条件的影像sceneID已经存储到IDs这个数组中了，直接把IDs的值复制到landsat download后面就行了。
```
import os
download_code='landsat download %s' % " ".join(IDs)
# os.system(download_code)
print(download_code)
```
把结果复制到命令行就可以下载了。

3.gdal

关于gdal的作用和强大之处网络上实在太多，这里不再赘述，事物都有双面性，所以我们要用辩证的思想看问题，gdal虽然强大，但也很复杂，它里面的大多数功能我们平时都用不到，所以下面我只介绍3个常用的命令（其实我只会这三个）。

3.1 gdalinfo

作用：查看栅格数据的基本信息
参数：gdalinfo
用法 gdalinfo filename.tif
在这里插入图片描述
结果包含该数据的所有信息，如坐标系、原点、像元大小、元数据、范围，NoData、数据类型等。
如果要将某个数据的信息保存到文件gdalinfo filename > info.txt

3.2 gdal_merge

作用：影像拼接。
参数：gdal_merge
这里我偷个懒，关于gdal_merge这篇文章写的很详细，大家可以参考。

3.3 gdalwarp

作用：数据投影和变形
参数：gdalwarp
用法1：影像裁剪

//规则切割
gdalwarp -te xmin,ymin,xmax,ymax input output 
//不规则切割
gdalwarp -cutline bound.shp input output

比如我现在有这样的数据：
在这里插入图片描述我要用这个多边形切割这个dem数据，写法如下：

gdalwarp -cutline bound.shp fdem.tif newdem.tif

在这里插入图片描述 如果要批量切割多个影像，可以新建一个txt文档然后在里面写入gdalwarp语句，比如：

然后把这个文件扩展名修改为.bat，右击这个文件选择以管理员身份运行就可以完成切割了。

ps:windows系统默认不显示已知文件的扩展名，需要进行设置显示已知文件的扩展名才可以修改txt的扩展名。

用法2：投影转换

gdalwarp -s_src srs_def -t_srs srs_def input output

- s_srs 原数据的坐标
-t_srs要转换的坐标系
坐标的定义方式是epsg、proj4或wkt格式（如果对此有疑问请自行查阅相关资料），通常来说我记住某个坐标的EPSG码或proj4内容是不容易，其实也并不需要记忆，因为它们很容易可以查到。如果你熟悉QGIS，那么通过QGIS不论是EPSG还是proj4都可以查到，另外，通过EPSG文件（下载地址）也可以查到。

比如我们要将newdem.tif的坐标转换为Pseudo-Mercator坐标系（web数据常用的坐标），首先我们要找到原坐标的定义（EPSG或Proj4），上面用gdalinfo看过fdem.tif的信息，知道它的坐标是UTM_Zone_50N（如果没记下翻回去看），所以我们先找到UTM_Zone_50N这个坐标的定义，用任意文本编辑器（记事本）打开前面下载的epsg文件，查找UTM Zone 50N和Pseudo-Mercator
在这里插入图片描述
其中椭圆内的是EPSG码，矩形内的是proj4定义。所以投影转换的代码应该写成：

gdalwarp -s_srs EPSG:23850 -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs" newdem.tif projdem.tif

ps:一般情况下-s_srs可以不用定义。
gdal的功能非常强大，我上面写的还不到九牛一毛，就算是这三个工具，我写的也只是皮毛，如果想进一步了解gdal，网上的学习资料很多。

4.rasterio

rasterio是一个栅格数据的读写工具，其本质还是基于gdal写的，但是它的功能比较简单易懂，平时操作栅格数据用就它足够了。

4.1 打开数据

import rasterio
rst=rasterio.open('fdem.tif')

4.2 数据浏览和编辑

在这里插入图片描述

4.3更多

关于rasterio的更多内容大家还是访问官网吧，实在写不动了。

5.后记

关于以上内容大家都什么疑问可以通过下载的邮箱联系我，最后给大家安利一个写python的工具Jupyter Notebook，anaconda默认安装，大家可以在开始菜单找到。推荐推荐它的原因在于它可以实时显示代码结果（上图），这对于数据处理是非常有帮助的。

作者：zhang
邮箱：xtfge_0915@163.com
博客：https://blog.csdn.net/xtfge0915

裘千仞不会水上漂

关注

3
点赞
踩
22

收藏

觉得还不错? 一键收藏
打赏
0
评论
栅格数据处理

栅格数据处理工具及方法推荐1.前言2.软件安装和环境配置2.1 andconda安装2.2 安装栅格数据处理库3.数据准备3.1 landsat数据批量下载工具3.2 modis数据批量下载工具3.3 下载数据1.前言最近受好友之邀写一些地理数据处理的内容，虽然害怕水平有限写的东西误人子弟，浪费大学时间，奈何抵不住好友的柔情攻势，只能尽可以把自己觉得有用的东西写下来，希望能帮到需要的人。因为...
复制链接

扫一扫