5.3. 创建数据集

gdal不单是读取，还可以创建数据集。而且创建的过程也是很迅速快捷的。可以参考一下[[http://www.gdal.org/gdal_tutorial.html|这里]]的快速指南（在Techniques for Creating Files,Using !CreateCopy()和Using Create()这三段）

gdal创建数据集有两种方式。一种用Create方法，一种用!CreateCopy方法。所有支持创建新文件的驱动都支持!CreateCopy方法，但只有少数支持Create方法。判断是否支持Create或者!CreateCopy方法可以在驱动对象原数据当中检查DCAP_CREATE和DCAP_CREATECOPY的值。当然，要验证，还需要创建一个Driver。当然，要创建数据集更需要创建一个Driver！

{{{ #! python format = “GTiff” driver = gdal.GetDriverByName( format ) metadata = driver.GetMetadata() if metadata.has_key(gdal.DCAP_CREATE)

and metadata[gdal.DCAP_CREATE] == ‘YES’:

print ‘Driver %s supports Create() method.’ % format

if metadata.has_key(gdal.DCAP_CREATECOPY)

and metadata[gdal.DCAP_CREATECOPY] == ‘YES’:

print ‘Driver %s supports CreateCopy() method.’ % format

}}} 上面是例子，先创建一个!GeoTiff的数据集格式驱动driver，然后提取出驱动的原数据，看原数据中是否有DCAP_CREATE或者DCAP_CREATECOPY的属性，通过判断这两个属性的值，就可以指导这个驱动是否是支持用Create创建数据集还是支持用!CreateCopy创建数据集。

原文中有一条注意：许多驱动是只读的并且不支持Create和!CreateCopy 创建驱动时要创建的驱动名是什么，可以参考[[http://www.gdal.org/formats_list.html|这里]]。第二列就是驱动名。第三列是是否支持创建数据集。如果不支持创建，那么也就无法调用Create或者!CreateCopy了。 使用!CreateCopy比较简单。就是把一个格式的图像直接转化为另一个格式的图像。

看看!CreateCopy的原型：

{{{ #! python >>> help(driver.CreateCopy) Help on method CreateCopy in module gdal: CreateCopy(self, filename, source_ds, strict=1, options=[], callback=None, callb ack_data=None) method of gdal.Driver instance >>> }}} 第一个参数是源文件的名称，第二个是目标文件名称。第三以后的参数都可选。如果不输入也可以。程序按照默认的方式运行。第三个参数取值是0或者1(True或者False)。取值为非的时候说明即使不能精确匹配地由原数据转化为目标数据，程序也照样执行!CreateCopy方法，不会产生致命错误。这种错误有可能是输出格式不支持输入数据格式象元的数据类型，或者是目标数据不支持写入空间参考等等。 文中举个例子:

{{{ #! python src_ds = gdal.Open( src_filename ) dst_ds = driver.CreateCopy( dst_filename, src_ds, 0 ) }}} 文中接着又举了个更复杂的例子。来说明第四个参数，第四个参数是指在格式转换中所要用到的一些特殊的参数。

{{{ #! python src_ds = gdal.Open( src_filename ) dst_ds = driver.CreateCopy( dst_filename, src_ds, 0,

[ ‘TILED=YES’, ‘COMPRESS=PACKBITS’ ] )

}}} 这个例子就是说明在转换Tiff的过程中用瓦片式存储代替条带式存储。用的压缩方法是PACKBITS。第四个参数根据各种格式都可能有不同的选项。所以无法统一列出，在[[http://gdal.maptools.org/frmt_gtiff.html|这里]]可以看到GTiff的全部创建参数。只有在根据各种不同格式时参考各自的参数和写法。 第五个和第六个参数是登记一个挂钩函数。使得可以把转换过程反映出来（比如画个进度条什么的）。这两个函数就不多说了。

如果你不存在一个已有的文件来创建新的文件，就需要用到Create方法了。它可以把建立在内存中的虚拟数据集创建到实际文件。

{{{ #! python >>> help(driver.Create) Help on method Create in module gdal: Create(self, filename, xsize, ysize, bands=1, datatype=1, options=[]) method of gdal.Driver instance >>> }}} 可以看出来，这个函数和!CreateCopy很像，不过它多了几个参数，xsize，ysize是图像大小，bands是图像的波段(通道)数，datatype就是图像的象元数据类型了。 dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte ) 这句就创建了一个图像。宽度是512*512,单波段，数据类型是Byte这里要注意，它少了源数据，因为这里用不到源数据。它创建了一个空的数据集。 要使它有东西，需要用其他步骤往里头塞东西。

{{{ #! python import Numeric, osr dst_ds.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] ) srs = osr.SpatialReference() srs.SetUTM( 11, 1 ) srs.SetWellKnownGeogCS( ‘NAD27’ ) dst_ds.SetProjection( srs.ExportToWkt() ) raster = Numeric.zeros( (512, 512) ) dst_ds.GetRasterBand(1).WriteArray( raster ) }}} 这里的例子往里面塞了个空间范围和坐标系(关于这些下面的文章再重点介绍)，还有一个512*512的全部都是0的数组数据。当然这样做，你在打开数据的时候只能看到一片黑乎乎。如果你要使你的数据有点实质性参考价值的话，你要把一个丰富多彩的数组塞到图片当中。当然这最好的办法就是从已有的数据中读取了。当然，出了从原数据中老老实实读出数据后，还可以进行矩阵魔术。使其满足我们的需要。然后再写到磁盘上。

不是每个数据都支持Create，比如JPG，PNG。所以，要写入内存中的数据，[[http://lists.maptools.org/pipermail/gdal-dev/2002-January/003254.html|官方建议]]用Create先建立一个过渡用的tiff，再用!CreateCopy来创建JPG或者PNG。(有点麻烦哦)

== 注意 == 在转换GTiff的过程中，需要注意一点。利用默认的参数在生成gtiff的过程中有可能出问题。比如用

{{{ #! python dataset = gdal.Open(“e:/gisdata/gtif/sd.tif”) width = dataset.RasterXSize height = dataset.RasterYSize data = dataset.ReadAsArray(0,0,width,height) driver = gdal.GetDriverByName(“GTiff”) driver.CreateCopy(“e:/sd.tif”,dataset,0) }}} 代码转出的GTiff文件虽然在ESRI的系列软件和其他一些看图程序中可以正常显示，但是在windows图像浏览器中不能正常显示，更重要的是在java的jai中不能正常显示。究其原因，是GDAL在导出的时候把284号域（!PlanarConfiguration域）设成了2，也就是RRRR……，GGGG……，BBBB……模式显示。但是在一些软件中只认值1，也就是RGBRGBRGBRGB……，所以上面的代码需要修改成

{{{ #! python dataset = gdal.Open(“e:/gisdata/gtif/sd.tif”) width = dataset.RasterXSize height = dataset.RasterYSize data = dataset.ReadAsArray(0,0,width,height) driver = gdal.GetDriverByName(“GTiff”) driver.CreateCopy(“e:/sd.tif”,dataset,0,[“INTERLEAVE=PIXEL”]) }}} 默认的INTERLEAVE是BAND（tags[284]=2），我们把它改成PIXEL(tags[284]=1)。这样就可以正常显示了。 更多的创建参数看[[http://gdal.maptools.org/frmt_gtiff.html|这里]]。

== 小例子 == 下面介绍一个建立3波段GTiff的小例子。当然凭空想个数据出来是很难的一件事情，除非我可以飞到几万米高空瞻仰地球母亲:)，于是我只好从另一个GTiff中读数据出来然后保存成另一个GTiff文件，做个意思吧。

{{{ import gdal import Numeric dataset = gdal.Open(“e:/gisdata/gtif/sd.tif”) width = dataset.RasterXSize height = dataset.RasterYSize datas = dataset.ReadAsArray(0,0,width,height) driver = gdal.GetDriverByName(“GTiff”) tods = driver.Create(“e:/gisdata/sd2.tif”,width,height,3,options=[“INTERLEAVE=PIXEL”]) tods.WriteRaster(0,0,width,height,datas.tostring(),width,height,

band_list=[1,2,3])

5.3.1. }}}

这是一个很简单的另存遥感图像的方法（不包括空间信息）。这里尤其注意Create函数中的options= [ “ INTERLEAVE=PIXEL “ ] 参数。没有这个参数，波段像素组织会错。保存出的图像只有横向的1/3。而且彩色完全不对 。

当然datas可以另外组织，比如这样也可以

{{{ datas = dataset.ReadAsArray(0,0,width,height) datas = Numeric.concatenate(datas) }}} —-

当然从波段里读取数据再拼接成完整的图像更是可以的。

{{{ datas = [] for i in range(3):

band = dataset.GetRasterBand(i+1) data = band.ReadAsArray(0,0,width,height) datas.append(Numeric.reshape(data,(1,-1)))

datas = Numeric.concatenate(datas) }}} —-

注意：从波段中读取的数组要拼接组织成RRR…GGG…BBB…这种形式才可以正确导出。不然整个图像看起来就像冥王星上的地形图。(这里有一点很奇怪，既然是PIXEL组织的，居然是RRR…GGG…BBB…形式而不是通常认为的RGB,RGB,RGB…形式)。如果你需要各个波段输入的话，可以循环到各个波段中，然后用Band对象的!WriteRaster来操作，而非在Dataset中调用!WriteRaster。

下面再看看空间参考如何导入。比如我们导入一个NAD27的空间参考，我们可以这样写

{{{ import osr tods.SetGeoTransform( [ 444720, 30, 0, 3751320, 0, -30 ] ) srs = osr.SpatialReference() srs.SetUTM( 11, 1 ) srs.SetWellKnownGeogCS( ‘NAD27’ ) tods.SetProjection( srs.ExportToWkt() ) }}} —-

呵呵，于是这个Tiff就变成了GTiff。空间参考系统是NAD27，起点（444720，3751320），像素大小30的TIFF了。当然直接这样写肯定是不对的。空间转换参数要进行配准运算，然后确定。我们这里就只是写个意思，说明可以这样写入文件罢了。我这里用的空间参考是美国常用的。至于中国的空间参考，你要西安还是北京，就看你高兴了。当然前提是先配准。

= 反馈 = 如果您发现我写的东西中有问题，或者您对我写的东西有意见，请登陆[[http://groups.google.com/group/geosings?hl=zh-CN|这个论坛]]。