code: def ds_getma_sub(src_ds, bnum=1, scale=None, maxdim=1024., return_ds=False):    
      
      #print src_ds.GetFileList()[0]
      b = src_ds.GetRasterBand(bnum)
      b_ndv = get_ndv_b(b)
      ns, nl, scale = get_sub_dim(src_ds, scale, maxdim)
      #The buf_size parameters determine the final array dimensions
      b_array = b.ReadAsArray(buf_xsize=ns, buf_ysize=nl)
      bma = np.ma.masked_values(b_array, b_ndv)
      out = bma
      if return_ds:
          dtype = src_ds.GetRasterBand(1).DataType
          src_ds_sub = gdal.GetDriverByName('MEM').Create('', ns, nl, 1, dtype)
          gt = np.array(src_ds.GetGeoTransform())
          gt[[1,5]] = gt[[1,5]]*scale
          src_ds_sub.SetGeoTransform(list(gt))
          src_ds_sub.SetProjection(src_ds.GetProjection())
          b = src_ds_sub.GetRasterBand(1)
          b.WriteArray(bma)
          b.SetNoDataValue(b_ndv)
          out = (bma, src_ds_sub)
      return out
docs: Load a subsampled array, rather than full resolution
  
      This is useful when working with large rasters
  
      Uses buf_xsize and buf_ysize options from GDAL ReadAsArray method.
  
      Parameters
      ----------
      ds : gdal.Dataset 
          Input GDAL Datset
      bnum : int, optional
          Band number
      scale : int, optional
          Scaling factor
      maxdim : int, optional 
          Maximum dimension along either axis, in pixels
      
      Returns
      -------
      np.ma.array    
          Masked array containing raster values