code: def get_res(ds, t_srs=None, square=False):
      
      gt = ds.GetGeoTransform()
      ds_srs = get_ds_srs(ds)
      #This is Xres, Yres
      res = [gt[1], np.abs(gt[5])]
      if square:
          res = [np.mean(res), np.mean(res)]
      if t_srs is not None and not ds_srs.IsSame(t_srs):
          if True:
              #This diagonal approach is similar to the approach in gdaltransformer.cpp
              #Bad news for large extents near the poles
              #ullr = get_ullr(ds, t_srs)
              #diag = np.sqrt((ullr[0]-ullr[2])**2 + (ullr[1]-ullr[3])**2)
              extent = ds_extent(ds, t_srs)
              diag = np.sqrt((extent[2]-extent[0])**2 + (extent[3]-extent[1])**2)
              res = diag / np.sqrt(ds.RasterXSize**2 + ds.RasterYSize**2)
              res = [res, res]
          else:        
              #Compute from center pixel
              ct = osr.CoordinateTransformation(ds_srs, t_srs)
              pt = get_center(ds)
              #Transform center coordinates
              pt_ct = ct.TransformPoint(*pt)
              #Transform center + single pixel offset coordinates
              pt_ct_plus = ct.TransformPoint(pt[0] + gt[1], pt[1] + gt[5])
              #Compute resolution in new units 
              res = [pt_ct_plus[0] - pt_ct[0], np.abs(pt_ct_plus[1] - pt_ct[1])]
      return res
docs: Get GDAL Dataset raster resolution