code: def get_xy_grids(ds, stride=1, getval=False):
      
      gt = ds.GetGeoTransform()
      #stride = stride_m/gt[1]
      pX = np.arange(0, ds.RasterXSize, stride)
      pY = np.arange(0, ds.RasterYSize, stride)
      psamp = np.meshgrid(pX, pY)
      mX, mY = pixelToMap(psamp[0], psamp[1], gt)
      return mX, mY
docs: Return 2D arrays of x and y map coordinates for input GDAL Dataset