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