Python|10m土地覆盖怎么找到自己需要的

2021-12-19 Python

背景

由清华大学地球系统科学系宫鹏教授团队研发的重大成果世界首套“2017年10米分辨率全球土地覆盖产品（FROM-GLC10）（2017年的为http://data.ess.tsinghua.edu.cn/fromglc2017v1.html，2015年的为http://data.ess.tsinghua.edu.cn/fromglc2015_v1.html）是很不错的土地覆盖产品。打开2017年数据的网站内容如下，可以看出是按照经纬度分块提供的。如果手动找自己需要的那几块比较费时，本文使用python的ogr库根据输入的矢量边界寻找对应的土地覆盖。

文章图片

步骤1：获取每个tif文件对应的矢量

首先想到的是将每个影像下载下来，然后批量提取其影像边界，但这样又要耽误不少功夫。下载几个文件发现，每个文件基本是按照10度进行分块，文件名记录了左上角的经纬度，像元分辨率为0.00025度，一般情况下一个文件的行列数为40000行40000列，但是在北纬85度时，行数为20000行。根据这些信息，将网页上的文件名列表复制为文本文件，通过Python的ogr库将生成这些文件的边界的矢量。将生成的矢量（fromglc2017.shp和fromglc2015.shp）叠加到地图上显示。
【Python|10m土地覆盖怎么找到自己需要的】2017年的覆盖图

文章图片
2015年的影像覆盖图

文章图片
2015年的没有85N的

步骤2：根据矢量筛选相应的文件名

有了上一步的每个文件对应的矢量，再筛选就好办了，一个简单思路是直接用ArcGIS将研究区的矢量和土地覆盖的范围进行按位置选择，再进行相应的下载即可。

文章图片

文章图片
将选择的内容导出到txt或excel中，再通过excel将网址提取出来下载。将url粘贴到批量下载的软件中就可以下载了。

文章图片
还有一个思路就是使用Python的ogr库做相交，再提取出响应的文件名。使用方法为python test2.py your_shp_path.shp。需要注意的是这里的shp必须是不带投影的wgs84坐标系。运行结果入下图，将控制台的地址粘贴到批量下载的软件中就可以下载了。
test2.py只针对2017年的有效，对于2015年的，读者可自行修改

文章图片

数据获取

获取上文生成的矢量和脚本的方法，关注微信公众号小猿猴GISer，回复2017下载

文章图片

推荐阅读

上一篇：Java笔记-13 http服务器

下一篇：Google earth engine对Landsat8辐射定标