基于矢量数据的大数据量遥感影像裁剪方法


哈哈，还是师兄比较猛！（说实话，师兄还不知道我转他的呢。算了，算是帮他宣传下咯，也学习了！）转http://user.qzone.qq.com/604888793ptlang=2052#!app=2&via=QZ.HashRefresh&pos=1322410365

Ben

rs.hxd@163.com

题记：如果一个遥感影像数据处理的程序不能够满足大影像处理的话，那么，基本可以认定这个程序是废的。

最近，在帮一位师姐做图像处理时，遇到了不少的麻烦，她的要求是处理大数据量的影像，而且速度不能太慢。对于我们这些业余的“码奴”来讲的确是一件很不容易的事情。她提出了三个需求，其中的一个需求就是使用矢量数据对遥感影像进行裁剪。矢量数据使用的是中国行政区图，而要裁剪的影像的大小是要1G以上的。刚看到这个数据的大小的时候，首先让我感受的不是信心倍增而是心理充满着畏惧。我也试想，对于我这种业余的“程序猿”来说，怎么能够完成？？。可是，既然答应了，那么就要一定做好。好吧！！
那就从最简单的开始吧。

我们知道中国政区图使用的是矢量图，里面存储的是多边形信息，也就是每个省的矢量化后的范围，那么最简单的裁剪办法是什么呢？就是对影像的每个点来判断是否在这所有的多边形内，这个最容易不过了，呵呵，以为找到了答案，欣喜若狂！！！，可是等把代码写完，满怀信心的双击执行文件后，我才发现，得到的不是程序运行结束的提示符，而是不知道要运行到什么时候的程序等待符。OMG！！！，我几乎抓狂，别说把代码给师姐了，就连我自己都觉的这代码可以说是“垃圾”。不管了，让它运行吧，看它能支持多久，午饭过后，回来再看一下，天哪，还在运行。睡一觉吧，午觉醒来后，在看一眼，终于运行完了，整整用了40多分钟啊！！！没办法了，那只能寻求高效的方法了。

提到这个高效的方法，我还真是有些“胆大”，那么，我是怎么“胆大”的呢？就让我们开始吧，呵呵！！

既然一副数据量那么大，那么何必不按块进行处理呢？干脆影像直接按照512*512的块大小处理吧，呵呵，跨出了程序高效地第一步。那么下一步怎么样呢？我思考了一下，按点进行处理的低效的原因在于大量的点与多边形拓扑关系的运算，那么提高效率的方法就在于尽量避免大量的点与多边形的拓扑关系的运算。基于该出发点，首先是让影像块与中国政区矢量图的多边形进行做包含运算，如果该影像块的在政区图的多边形内的话，那么直接跳过该数据的处理，进行下一个影像块的处理。其次，对于不包含在政区图内的影像块进行一次求交运算，如果没有一个行政图的多边形与其相交，那么，直接跳过，进行下一个影像块的处理。第三，如果行政区内的影像有一个或多个与该影像块相交的话，那么，求取该影像块与该多边形的交。第四，虽然相交，但是还是需要将影像块的所有点与多边形进行运算，这样的话，内部还是包含了大量的点。呵呵，没办法，继续优化吧。谁让咱们的终止是减少与点与多边形的运算呢？该怎么处理呢？，这样吧，我们取求交后的多边形的外包矩形吧。呵呵，这样数据量少了吧？只对外包矩形内的点进行与多边形的求交运算。哈哈......，所有的思路理清楚了。可以完事了吗？

根本不可能，为什么？这里面又是多边形求交，又是点与多边形的处理，我一个学遥感的，怎么能在那么短时间里写这样的拓扑关系的代码？没办法了，上网搜吧。还好，我对GDAL很熟悉，哈哈....,优势来了，GDAL里面的OGR是矢量模块的运算库，但是，让我失望的是，OGR只是提供了数据读取的一半接口，而矢量处理却依赖于GEOS库。那就用采用GEOS库吧，事情依然不是那么顺利，一个小小的库居然花了我1天的时间来编译，黄天不负有心人啊，终于大功告成！！！事情还没完呢？haha。。。。下一步工作是要写代码了，逻辑清楚了，代码还在话下吗？？？？？

测试一下效率吧，惊人啊！！！！！！！！！！！！！！！！！！！！！！！！！，7800*2500的影像的裁剪（包括数据读写）不到10秒钟就裁剪完了。不仅你不相信，就是我也有些不相信自己的眼睛......。效率总算提上去了，虽然程序还是有些小BUG，慢慢地在维护吧...



PS：

1. 该思路没有参考任何软件，纯自己的想法（可能会其他的软件于此雷同吧，这个我也不清楚）

2. 有同学建议我使用某国产GIS软件的二次开发库，为什么不用它的开发包？ 这个“你懂得”！！



3. GDAL(栅格处理库) + OGR（矢量处理库） + GEOS（矢量关系库） + PROJ4（地图投影库）， 四库合璧，天下无敌....