同じ名前の画像から単一のVRTを作成する

Aug 23 2020

同じ日付のラスターのリストがあり、同じ日付のVRTに変換する必要があります。問題は、Gdalリストにループを作成すると、同じ名前のVRTが上書きされることです。同じVRTに追加したいのですが。例えば:

root = 'F:...\\MSAVI'
out_vrt = 'F:...\\out_vrt\\'

raster = [os.path.join(root, file) for root, dirs, files in os.walk(root) 
               for file in files if os.path.splitext(file)[1] == '.tif' ]
raster 
['.../33TVE6/S2_20190712_079_33TVE6_B.tif',
 '.../32TQQ3/S2_20190712_079_33TQQ3_B.tif',
  ...] 

dates = []
for file in raster:
    dates.append(file.split('\\')[6].split('_')[1])

dates
['20190817',
 '20190901',
 ...]


for day in dates:
    for file in raster:
        data_path = file.split('\\')[6].split('_')[1]
        name = 'S2x2A_' + day + '_MSAVI.vrt'
        tile = file.split('\\')[6].split('_')[3]
        
        if day in data_path:

            my_vrt = gdal.BuildVRT(out_vrt+ name, file)
            my_vrt = None

出力は、ラスターが1つしかないVRTです。

<VRTDataset rasterXSize="4001" rasterYSize="4001">
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["WGS 84 / UTM zone 33N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32633"]]</SRS>
  <GeoTransform>  4.3500000000000000e+05,  1.0000000000000000e+01,  0.0000000000000000e+00,  4.4983300000000000e+06,  0.0000000000000000e+00, -1.0000000000000000e+01</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1">
    <NoDataValue>-32768</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">F:...\33TVE6\S2_20190712_079_33TVE6_B_MSAVI.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="4001" RasterYSize="4001" DataType="Int16" BlockXSize="4001" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="4001" ySize="4001" />
      <DstRect xOff="0" yOff="0" xSize="4001" ySize="4001" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
  </VRTRasterBand>
</VRTDataset>

以前のVRTを上書きせずに、同じ日付のラスターを追加して、より多くのラスターを持つVRTを取得するにはどうすればよいですか?

回答

3 IanTurton Aug 23 2020 at 17:49

BuildVRT 2番目(ソース)パラメーターとしてファイルのリストを取得できるため、日付に一致するすべてのファイルを(フルパスを使用して)検索し、そのリストを関数に渡す必要があります。

だから次のようなもの:

for day in dates:
  f = [file for file in rasters if rasters.contains(day)]
  my_vrt = gdal.BuildVRT(out_vrt+name, f)