while循环中函数参数未更新引发无限循环的问题排查

问题描述 #

执行内嵌带参数函数的while循环时，发现循环内的参数未更新，进而导致无限循环。该while循环应在pointInRect函数返回false时终止，但却重复执行相同操作，疑问是gt值未更新还是被重新初始化了？

相关代码 #

import shapefile
from osgeo import gdal
#import rasterio

print(gdal.VersionInfo())
def pointInRect(x, y, x1, y1, w, h): # check if a raster point is in another raster
    x2, y2 = x1+w, y1+h
    if (x1 < x and x < x2):
        if (y1 < y and y < y2):
            return True
    return False

# Open the shapes centroids
shp_cntrds = 'Path to centroids'
sf_cntrds = shapefile.Reader(shp_cntrds)
shapes_cntrds = sf_cntrds.shapes()
records_cntrds = sf_cntrds.records()

# adjust labels position according to its shapes centroids position
for i in range(len(records_cntrds)):
    print(i)

    tods = gdal.Open(str(records_cntrds[i][1]))

    width = tods.RasterXSize
    height = tods.RasterYSize

    tods.SetGeoTransform([shapes_cntrds[i].points[0][0] - (width * 0.005), 0.01, 0,
                          shapes_cntrds[i].points[0][1] + (height * 0.005), 0, -0.01])

    gt = tods.GetGeoTransform()
    left = gt[0]
    bottom = gt[3] + width * gt[4] + height * gt[5]
    right = gt[0] + width * gt[1] + height * gt[2]
    top = gt[3]
    srs = osr.SpatialReference()
    srs.SetUTM(32, 1)  # set crs
    srs.SetWellKnownGeogCS('WGS84')  # set crs
    tods.SetProjection(srs.ExportToWkt())  # set Projection and save file
    print(width, height)
    tods = None

# iterate through Labels and move labels away from each others  if they overlapp
for i in range(len(records_cntrds)):
    tods1 = gdal.Open(str(records_cntrds[i][1])) # records of the centroid shapefile contains the raster file path

    width = tods1.RasterXSize
    height = tods1.RasterYSize
    gt = tods1.GetGeoTransform()
    left = gt[0]
    bottom = gt[3] + width * gt[4] + height * gt[5]
    right = gt[0] + width * gt[1] + height * gt[2]
    top = gt[3]
    face = [x for x in list(range(len(records_cntrds))) if x != i]
    tods1 = None
    for j in face:
        if str(records_cntrds[i][1]) == str(records_cntrds[j][1]):
            pass
        else:
            ds_raster_face = gdal.Open(str(records_cntrds[j][1]))
            gt_face = ds_raster_face.GetGeoTransform()
            width_face = ds_raster_face.RasterXSize
            height_face = ds_raster_face.RasterYSize
            left_face = gt_face[0]
            bottom_face = gt_face[3] + width_face * gt_face[4] + height_face * gt_face[5]
            right_face = gt_face[0] + width_face * gt_face[1] + height_face * gt_face[2]
            top_face = gt_face[3]

            width1 = width
            left1 = left
            height1 = height
            bottom1 = bottom
            while pointInRect(left_face, bottom_face, left1, bottom1, width1*0.01, height1*0.01) :
                tods2 = gdal.Open(str(records_cntrds[i][1]))
                gt = tods2.GetGeoTransform()
                width1 = tods2.RasterXSize
                height1 = tods2.RasterYSize
                left1 = gt[0]
                bottom1 = gt[3] + width1 * gt[4] + height1 * gt[5]
                print("while executed")
                tods2.SetGeoTransform([(shapes_cntrds[i].points[0][0] - (width1 * 0.005)) - 2.7, 0.01, 0,
                                  (shapes_cntrds[i].points[0][1] + (height1 * 0.005)) - 2.8, 0, -0.01])
                print("coordinates changed to",(i, left1, bottom1, width1, height1))
                tods2 = None

问题原因 #

1. GDAL文件打开模式错误：默认gdal.Open()以只读模式打开文件，调用SetGeoTransform()后修改不会写入磁盘，下次打开仍读取旧的GeoTransform值，导致参数无法真正更新。
2. 循环内参数更新逻辑错误：在while循环中，修改GeoTransform后没有重新获取最新的参数更新left1、bottom1，而是使用修改前读取的gt值，导致循环条件始终满足。
3. 未确保修改写入磁盘：修改GeoTransform后没有调用FlushCache()强制写入，即使开启更新模式，也可能因缓存机制导致磁盘文件未及时更新。
4. 缺失必要导入：代码中使用了osr.SpatialReference()但未导入osr模块，会导致运行报错。

修复方案 #

1. 添加缺失的导入 #

在代码开头补充：

from osgeo import osr

2. 修改GDAL文件打开模式 #

所有需要修改文件的地方，使用gdal.GA_Update参数以更新模式打开：

# 第一个循环中修改投影和GeoTransform时
tods = gdal.Open(str(records_cntrds[i][1]), gdal.GA_Update)
# while循环中修改GeoTransform时
tods2 = gdal.Open(str(records_cntrds[i][1]), gdal.GA_Update)

3. 修正while循环内参数更新逻辑 #

修改GeoTransform后，重新获取最新的GeoTransform值并更新循环条件依赖的变量：

while pointInRect(left_face, bottom_face, left1, bottom1, width1*0.01, height1*0.01) :
    tods2 = gdal.Open(str(records_cntrds[i][1]), gdal.GA_Update)
    width1 = tods2.RasterXSize
    height1 = tods2.RasterYSize
    # 计算新的GeoTransform参数
    new_gt_x = (shapes_cntrds[i].points[0][0] - (width1 * 0.005)) - 2.7
    new_gt_y = (shapes_cntrds[i].points[0][1] + (height1 * 0.005)) - 2.8
    tods2.SetGeoTransform([new_gt_x, 0.01, 0, new_gt_y, 0, -0.01])
    # 强制写入磁盘
    tods2.FlushCache()
    # 更新循环条件需要的参数
    gt_updated = tods2.GetGeoTransform()
    left1 = gt_updated[0]
    bottom1 = gt_updated[3] + width1 * gt_updated[4] + height1 * gt_updated[5]
    print("while executed")
    print("coordinates changed to",(i, left1, bottom1, width1, height1))
    tods2 = None

内容的提问来源于stack exchange，提问作者Houssem