## 用ogr和PIL把矢量数据转化成栅格图像

### 求索

PIL和GDAL可以互相共享数据，互相协调互相补充，这很好。现在再加上PIL对ogr的矢量绘图支持，真的是－－

－－“Perfect！”

(似乎这一篇的名称可以改成“江湖”，哈哈，说笑，小说看多了)

### 身手

```   1 Example: Draw a Grey Cross Over an Image
2 import Image, ImageDraw
3 im = Image.open("lena.pgm")
4 draw = ImageDraw.Draw(im)
5 draw.line((0, 0) + im.size, fill=128)
6 draw.line((0, im.size[1], im.size[0], 0), fill=128)
7 del draw
8 # write to stdout
9 im.save(sys.stdout, "PNG")
```

```   1 im.save("lena.png","PNG")
```

```   1 import Image, ImageDraw
2 import sys
3 import ogr
4 from Numeric import *
```

```   1 ds = ogr.Open("e:/gisdata/data/streets.shp")
2 layer = ds.GetLayer()
3 extent = layer.GetExtent()
4 geowidth = extent[1]-extent[0]
5 geoheight = extent[3]-extent[2]
6 width = 1000
7 height = 1000/geowidth*geoheight
8 scale = geowidth/width
```

```   1 # create a new png file
2 imn = Image.new("RGB",(width,height),(255,255,255))
3 imn.save("test.pgm","PNG")
4 del imn
```

Image.new并不能写入数据到硬盘上，所以要用imn.save写入，这里保存的是PNG的格式。

```   1 import Image, ImageDraw
2 import sys,math
3 import ogr
4 from Numeric import *
5 class PILImgDC:
6     def __init__(self,layer):
7         self.layer = layer
8         self.features = [layer.GetFeature(i) for i in range(layer.GetFeatureCount())]
9         FeatureDrawMap = {
10             #ogr.wkbPoint:self.__DrawPoint,
11             #ogr.wkbMultiPoint:self.__DrawPoints,
12             ogr.wkbLineString:self.__DrawLine,
13             ogr.wkbMultiLineString:self.__DrawLines,
14             #ogr.wkbPolygon:self.__DrawPolygon,
15             #ogr.wkbMultiPolygon:self.__DrawPolygons
16             }
17         self.featureDrawMap = FeatureDrawMap
18         self.fill = None
19         self.iimgname = None
20     def SetFillColour(self,color):
21         self.fill = color
22     def SetInputImg(self,imgname):
23         self.iimgname = imgname
24     def DrawAll(self,width,ofname):
25         if self.fill is None:
26             self.fill = (128,128,128)
27         extent = layer.GetExtent()
28         geowidth = math.fabs(extent[1]-extent[0])
29         geoheight = math.fabs(extent[3]-extent[2])
30         self.minx = min(extent[0],extent[1])
31         self.maxy = max(extent[3],extent[2])
32         height = width/geowidth*geoheight
33         self.scale = width/geowidth
34         im = Image.new("RGB",(width,height),(255,255,255))
35         self.draw = ImageDraw.Draw(im)
36         for feature in self.features:
37             geomref = feature.GetGeometryRef()
38             geomreftype = geomref.GetGeometryType()
39             drawFun = self.featureDrawMap[geomreftype]
40             drawFun(geomref)
41         del self.draw
42         im.save(ofname, "PNG")
43     def __GetPoints(self,georef):
44         pointcount = georef.GetPointCount()
45         xs,ys = [],[]
46         for i in range(pointcount):
47             xs.append(georef.GetX(i))
48             ys.append(georef.GetY(i))
49         xs = array(xs,Float); ys = array(ys,Float)
50         xs = (xs-self.minx)*self.scale
51         ys = (self.maxy-ys)*self.scale
52         xs.shape=(-1,1);ys.shape=(-1,1)
53         points = concatenate((xs,ys),1)
54         return points
55     def __DrawLine(self,georef):
56         points = self.__GetPoints(georef)
57         xys = reshape(points, (1,-1))[0]
58         self.draw.line(xys,fill=self.fill)
59     def __DrawLines(self,georef):
60         sublinecount = georef.GetGeometryCount()
61         for i in range(sublinecount):
62             self.__DrawLine(georef.GetGeometryRef(i))
63 if __name__=="__main__":
64     ds = ogr.Open("e:/gisdata/shp/lines.shp")
65     layer = ds.GetLayer()
66     dc = PILImgDC(layer)
67     dc.SetFillColour((100,100,0))
68     dc.DrawAll(1000,"output.png")
```