是否可以使用gdalbuildvrt和python将元数据存储在geotiff的.vrt文件中?

2024-06-26 13:44:39 发布

您现在位置:Python中文网/ 问答频道 /正文

我编写了一个python脚本,从.aux.xml文件读取元数据标记值,并使用gdalinfo将它们添加到相应的tiff文件中。然后使用gdalinfo,我可以在元数据中找到它们

现在我需要将这个tif文件转换为.vrt格式。使用gdalbuidlvrt,我意识到我缺少这些自定义字段

在转换为.vrt文件时,是否有办法保留它们

在下面,您可以找到我的代码的简化版本


import os
import gdal
from lxml import etree


dataset = gdal.Open('source_file.tif',gdal.GA_Update)



filename='C:\\Users\\' \
         '\\source_file.tif.aux.xml'
doc = etree.parse(filename)

CLASS_0_NO_DATA_xml = doc.xpath('.//MDI[@key="CLASS_0_NO_DATA[%]"]//text()')
CLASS_1_RAIN_xml = doc.xpath('.//MDI[@key="CLASS_1_RAIN[%]"]//text()')


band2_categories_xml = doc.xpath('.//PAMRasterBand[@band="2"]/CategoryNames/Category//text()')
band2_nodata_xml = band2_categories_xml[0]
band2_quality_index_xml = band2_categories_xml[1]

##### Note GetRasterBand() takes band number starting from 1 not 0
band = dataset.GetRasterBand(1)
CLASS_0_NO_DATA = 'CLASS_0_NO_DATA'
CLASS_0_NO_DATA_value = CLASS_0_NO_DATA_xml[0] # str(dataset.GetMetadataItem('CLASS_0_NO_DATA[%]'))

CLASS_1_RAIN = 'CLASS_1_RAIN'
CLASS_1_RAIN_value = str(CLASS_1_SNOW_xml[0]) #str(dataset.GetMetadataItem('CLASS_1_RAIN[%]'))


band.SetMetadata( {CLASS_0_NO_DATA: CLASS_0_NO_DATA_value, CLASS_1_RAIN: CLASS_1_RAIN_value
                   } )



###### All values must be strings
dataset.SetMetadata( {'Band_1': 'MAP', 'Band_2': 'QUALITY_FLAG'} )


metadata = os.popen('gdalinfo Output_file.tif').read()
print(metadata) 

这段代码的输出是

Driver: GTiff/GeoTIFF
Files: source_file.tif
Size is 4901, 2867
Coordinate System is:
LOCAL_CS["ETRS89 / ETRS-LAEA",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4258"]],
    AUTHORITY["EPSG","3035"],
    UNIT["metre",1]]
Origin = (3847097.795985520351678,2872447.744972674641758)
Pixel Size = (231.656358262999987,-231.656358262999987)
Metadata:
  AREA_OR_POINT=Area
  Band_1=MAP  
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( 3847097.796, 2872447.745)
Lower Left  ( 3847097.796, 2208288.966)
Upper Right ( 4982445.608, 2872447.745)
Lower Right ( 4982445.608, 2208288.966)
Center      ( 4414771.702, 2540368.355)
Band 1 Block=4901x1 Type=Byte, ColorInterp=Palette
  Metadata:
    CLASS_0_NO_DATA=12.388736109961544
    CLASS_1_RAIN=3.151026530394237

如您所见,我添加了两个类:CLASS_0_NO_DATA和CLASS_1_RAIN及其百分比

当我使用

gdabuildvrt -resolution highest -srcnodata 0 -vrtnodata 0 Output_file.vrt Output_file.tif

我明白了

 Driver: VRT/Virtual Raster
Files: Output_file.vrt
       Output_file.tif
Size is 4901, 2867
Coordinate System is:
LOCAL_CS["ETRS89 / ETRS-LAEA",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4258"]],
    AUTHORITY["EPSG","3035"],
    UNIT["metre",1]]
Origin = (3847097.795985520351678,2872447.744972674641758)
Pixel Size = (231.656358262999987,-231.656358262999987)
Corner Coordinates:
Upper Left  ( 3847097.796, 2872447.745)
Lower Left  ( 3847097.796, 2208288.966)
Upper Right ( 4982445.608, 2872447.745)
Lower Right ( 4982445.608, 2208288.966)
Center      ( 4414771.702, 2540368.355)
Band 1 Block=128x128 Type=Byte, ColorInterp=Palette
  NoData Value=0

其中有关这两个类的信息丢失

有什么想法吗


Tags: 文件nooutputdatabandxmlepsgdataset