Python中用GDAL實現向量對柵格的切割
概述:
本文講述如何在Python中用GDAL實現根據輸入向量邊界對柵格資料的裁剪。
效果:
裁剪前
向量邊界
裁剪後
實現程式碼:
# -*- coding: utf-8 -*- """ @author lzugis @date 2017-06-02 @brief 利用shp裁剪影像 """ from osgeo import gdal, gdalnumeric, ogr from PIL import Image, ImageDraw import os import operator gdal.UseExceptions() # This function will convert the rasterized clipper shapefile # to a mask for use within GDAL. def imageToArray(i): """ Converts a Python Imaging Library array to a gdalnumeric image. """ a=gdalnumeric.fromstring(i.tobytes(),'b') a.shape=i.im.size[1], i.im.size[0] return a def arrayToImage(a): """ Converts a gdalnumeric array to a Python Imaging Library Image. """ i=Image.frombytes('L',(a.shape[1],a.shape[0]), (a.astype('b')).tobytes()) return i def world2Pixel(geoMatrix, x, y): """ Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate the pixel location of a geospatial coordinate """ ulX = geoMatrix[0] ulY = geoMatrix[3] xDist = geoMatrix[1] pixel = int((x - ulX) / xDist) line = int((ulY - y) / xDist) return (pixel, line) # # EDIT: this is basically an overloaded # version of the gdal_array.OpenArray passing in xoff, yoff explicitly # so we can pass these params off to CopyDatasetInfo # def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ): ds = gdal.Open( gdalnumeric.GetArrayFilename(array) ) if ds is not None and prototype_ds is not None: if type(prototype_ds).__name__ == 'str': prototype_ds = gdal.Open( prototype_ds ) if prototype_ds is not None: gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff ) return ds def histogram(a, bins=range(0,256)): """ Histogram function for multi-dimensional array. a = array bins = range of numbers to match """ fa = a.flat n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins) n = gdalnumeric.concatenate([n, [len(fa)]]) hist = n[1:]-n[:-1] return hist def stretch(a): """ Performs a histogram stretch on a gdalnumeric array image. """ hist = histogram(a) im = arrayToImage(a) lut = [] for b in range(0, len(hist), 256): # step size step = reduce(operator.add, hist[b:b+256]) / 255 # create equalization lookup table n = 0 for i in range(256): lut.append(n / step) n = n + hist[i+b] im = im.point(lut) return imageToArray(im) def main( shapefile_path, raster_path ): # Load the source data as a gdalnumeric array srcArray = gdalnumeric.LoadFile(raster_path) # Also load as a gdal image to get geotransform # (world file) info srcImage = gdal.Open(raster_path) geoTrans = srcImage.GetGeoTransform() # Create an OGR layer from a boundary shapefile shapef = ogr.Open(shapefile_path) lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] ) poly = lyr.GetNextFeature() # Convert the layer extent to image pixel coordinates minX, maxX, minY, maxY = lyr.GetExtent() ulX, ulY = world2Pixel(geoTrans, minX, maxY) lrX, lrY = world2Pixel(geoTrans, maxX, minY) # Calculate the pixel size of the new image pxWidth = int(lrX - ulX) pxHeight = int(lrY - ulY) clip = srcArray[:, ulY:lrY, ulX:lrX] # # EDIT: create pixel offset to pass to new image Projection info # xoffset = ulX yoffset = ulY print "Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ) # Create a new geomatrix for the image geoTrans = list(geoTrans) geoTrans[0] = minX geoTrans[3] = maxY # Map points to pixels for drawing the # boundary on a blank 8-bit, # black and white, mask image. points = [] pixels = [] geom = poly.GetGeometryRef() pts = geom.GetGeometryRef(0) for p in range(pts.GetPointCount()): points.append((pts.GetX(p), pts.GetY(p))) for p in points: pixels.append(world2Pixel(geoTrans, p[0], p[1])) rasterPoly = Image.new("L", (pxWidth, pxHeight), 1) rasterize = ImageDraw.Draw(rasterPoly) rasterize.polygon(pixels, 0) mask = imageToArray(rasterPoly) # Clip the image using the mask clip = gdalnumeric.choose(mask, \ (clip, 0)).astype(gdalnumeric.uint8) # This image has 3 bands so we stretch each one to make them # visually brighter for i in range(3): clip[i,:,:] = stretch(clip[i,:,:]) # Save new tiff # # EDIT: instead of SaveArray, let's break all the # SaveArray steps out more explicity so # we can overwrite the offset of the destination # raster # ### the old way using SaveArray # # gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path) # ### # gtiffDriver = gdal.GetDriverByName( 'GTiff' ) if gtiffDriver is None: raise ValueError("Can't find GeoTiff Driver") gtiffDriver.CreateCopy( "beijing.tif", OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset ) ) # Save as an 8-bit jpeg for an easy, quick preview clip = clip.astype(gdalnumeric.uint8) gdalnumeric.SaveArray(clip, "beijing.jpg", format="JPEG") gdal.ErrorReset() if __name__ == '__main__': #shapefile_path, raster_path shapefile_path = 'beijing.shp' raster_path = 'world.tif' main( shapefile_path, raster_path )
-------------------------------------------------------------------------------------------------------------
技術部落格
CSDN:http://blog.csdn.NET/gisshixisheng
部落格園:http://www.cnblogs.com/lzugis/
線上教程
http://edu.csdn.Net/course/detail/799
Github
https://github.com/lzugis/
聯絡方式
q q:1004740957
e-mail:[email protected]
公眾號:lzugis15
Q Q 群:452117357(webgis)
337469080(Android)
相關推薦
Python中用GDAL實現向量對柵格的切割
概述:本文講述如何在Python中用GDAL實現根據輸入向量邊界對柵格資料的裁剪。效果:裁剪前向量邊界裁剪後實現程式碼:# -*- coding: utf-8 -*- """ @author lzugis @date 2017-06-02 @brief 利用shp裁剪影像 "
python中用list實現queue
python中用list實現queue class myqueue(): def push(self,a): if self.count>=self.maxsize: print('The queue has been full!') re
python中用logging實現日誌滾動和過期日誌刪除
用python中的logging庫實現日誌滾動和過期日誌刪除。 logging庫提供了兩個可以用於日誌滾動的class(可以參考https://docs.python.org/2/library/logging.handlers.html),一個是Rota
Arcgis+Python實現對柵格歸一化處理
影象歸一化就不多說了,就是(數值-min)/(max-min),把結果都劃歸到0-1範圍,便於不同變數之間的比較,取消了不同數量差別。 第一個方法,需要對柵格資料預先知道取值範圍。 第二種方法,比較好點,直接讀取屬性
python學習--如何實現可叠代對象(itearable)和叠代器(iterator)
dict 作用 pri 返回 -- 生成器 ble ear item 關於可叠代對象Iterable 和叠代器對象iterator 可叠代對象:可以直接作用於for循環的對象統稱為可叠代對象:Iterable。 可叠代對象包含一個__iter__方法,或__getitem_
python裝飾器實現對異常代碼出現進行監控
args lin sha lines 監控腳本 一秒 readline utf 發送 異常,不應該存在,但是我們有時候會遇到這樣的情況,比如我們監控服務器的時候,每一秒去采集一次信息,那麽有一秒沒有采集到我們想要的信息,但是下一秒采集到了, 而
[Python Study Notes]實現對鼠標控制
移動 bin log ont 實例化 html 實現 坐標 tool ‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘ &g
[Python Study Notes]實現對鍵盤控制與監控
實現 博客 pac rip art 字符串 line 一個 sys ‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘ >
IDL實現向量(shp)裁剪柵格TASK(一)
隨著ENVI/IDL版本的更新,IDL對向量和柵格資料的處理也變得越來越簡單化。其提供了很多方便的介面,使得使用者呼叫和學習練習便捷成為了可能。 最近接觸IDL,發現好多網上的程式碼都是延後的,新的
運維中的日誌切割操作梳理(Logrotate/python/shell指令碼實現)
對於Linux系統安全來說,日誌檔案是極其重要的工具。不知為何,我發現很多運維同學的伺服器上都執行著一些諸如每天切分Nginx日誌之類的CRON指令碼,大家似乎遺忘了Logrotate,爭相發明自己的輪子,這真是讓人沮喪啊!就好比明明身邊躺著現成的性感美女,大家卻忙著自娛自樂,罪過!logrotate程式是一
Python中使用面狀向量裁剪柵格影像,並依據Value值更改向量屬性
本文整體思路:在Python中使用Geopandas庫,依次讀取shp檔案的每一個面狀要素,獲取其空間邊界資訊並裁剪對應的柵格影像,計算所裁剪影像Value值的眾數,將其設定為對應面狀要素的NewTYPE值,所有要素屬性值都改好之後儲存為新的shp檔案。 使用Python處理空間資料確實用的不多,所以一個星
【初學Python】關於python實現Cameo對視訊的擷取和錄製
學習《OpenCV3計算機視覺,Python語言實現》時讀書筆記。 最近對計算機視覺產生了濃厚的興趣,又剛好在學校的圖書館看到了這本書,那可剛好滿足了我的興趣,今天先敲一個Cameo。 儲存圖片視訊時,直接通過當前日期進行儲存。 先建立一個manager.py imp
python中用裝飾器實現單例模式
def singleton(cls,*args,**kwargs): instances = {} def get_instance(*args,**kwargs): if cls not in instances: i
在Python中用while迴圈實現姓名的錄入
程式碼為 def get_formatted_name(first_name, last_name): full_name = first_name + ' ' + last_name return full_name.title() while True
SVM引數引數介紹以及python實現GA對SVM引數的優化
最近開始玩起了機器學習,以前都是用matlab做一些機器學習的東西,畢竟要真正放到工程上應用還是python用起來比較好,所以今天就開始學習下使用SVM進行迴歸(分類)預測。 SVM 使用的一般步驟是: 1)準備資料集,轉化為 SVM支援的資料格式 : [label] [ind
使用結巴分詞(jieba)對自然語言進行特徵預處理(Python、Java 實現)
一、前言 之前使用基於 Python 語言的 Spark 進行機器學習,程式設計起來是十分簡單。 ① 但是演算法部署到雲伺服器上,是一個障礙。 ② 得藉助 Flask/Django 等 Python W
利用Python網路爬蟲實現對網易雲音樂歌詞爬取
今天小編給大家分享網易雲音樂歌詞爬取方法。 本文的總體思路如下: 找到正確的URL,獲取原始碼; 利用bs4解析原始碼,獲取歌曲名和歌曲ID; 呼叫網易雲歌曲API,獲取歌詞; 將歌詞寫入檔案,並存入本地。 本文的目的是獲取網易雲音樂的歌詞,並將歌詞存入到本地檔案。整
Python:用Numpy來實現向量的各種運算
首先要寫上這一句: from numpy import * (寫上這句的前提也得你已經安了numpy) (1) 定義一個零向量(4維): >>>a=zeros(4) >>>a array([0.,0.,0.,0.]) 定義一個List:
python練習:實現一個整數數組裡面兩個數之和為183的所有整數對
1 l1 = [183,0,1,2,-184,367] 2 3 num = [] 4 5 for i in range (0,len(l1)): 6 7 for l in range (i+1,len(l1)): 8 9 if l1[i]+l1[l]==
[Python程式設計]PyMongo實現對JSON的匯入和匯出
JSON匯入 1. 開啟Collection import json import pymongo client = pymongo.MongoClient('localhost') db = cli