docs: def extract(self, px, py, radius=1, nodata=None):
          """
          extract weighted average of pixels intersecting with a defined radius to a point.
  
          Parameters
          ----------
          px: int or float
              the x coordinate in units of the Raster SRS
          py: int or float
              the y coordinate in units of the Raster SRS
          radius: int or float
              the radius around the point to extract pixel values from; defined as multiples of the pixel resolution
          nodata: int
              a value to ignore from the computations; If `None`, the nodata value of the Raster object is used
  
          Returns
          -------
          int or float
              the the weighted average of all pixels within the defined radius
  
          """
          if not self.geo['xmin'] <= px <= self.geo['xmax']:
              raise RuntimeError('px is out of bounds')
          
          if not self.geo['ymin'] <= py <= self.geo['ymax']:
              raise RuntimeError('py is out of bounds')
          
          if nodata is None:
              nodata = self.nodata
          
          xres, yres = self.res
          
          hx = xres / 2.0
          hy = yres / 2.0
          
          xlim = float(xres * radius)
          ylim = float(yres * radius)
          
          # compute minimum x and y pixel coordinates
          xmin = int(floor((px - self.geo['xmin'] - xlim) / xres))
          ymin = int(floor((self.geo['ymax'] - py - ylim) / yres))
          
          xmin = xmin if xmin >= 0 else 0
          ymin = ymin if ymin >= 0 else 0
          
          # compute maximum x and y pixel coordinates
          xmax = int(ceil((px - self.geo['xmin'] + xlim) / xres))
          ymax = int(ceil((self.geo['ymax'] - py + ylim) / yres))
          
          xmax = xmax if xmax <= self.cols else self.cols
          ymax = ymax if ymax <= self.rows else self.rows
          
          # load array subset
          if self.__data[0] is not None:
              array = self.__data[0][ymin:ymax, xmin:xmax]
              # print('using loaded array of size {}, '
              #       'indices [{}:{}, {}:{}] (row/y, col/x)'.format(array.shape, ymin, ymax, xmin, xmax))
          else:
              array = self.raster.GetRasterBand(1).ReadAsArray(xmin, ymin, xmax - xmin, ymax - ymin)
              # print('loading array of size {}, '
              #       'indices [{}:{}, {}:{}] (row/y, col/x)'.format(array.shape, ymin, ymax, xmin, xmax))
          
          sum = 0
          counter = 0
          weightsum = 0
          for x in range(xmin, xmax):
              for y in range(ymin, ymax):
                  # check whether point is a valid image index
                  val = array[y - ymin, x - xmin]
                  if val != nodata:
                      # compute distances of pixel center coordinate to requested point
                      
                      xc = x * xres + hx + self.geo['xmin']
                      yc = self.geo['ymax'] - y * yres + hy
                      
                      dx = abs(xc - px)
                      dy = abs(yc - py)
                      
                      # check whether point lies within ellipse: if ((dx ** 2) / xlim ** 2) + ((dy ** 2) / ylim ** 2) <= 1
                      weight = sqrt(dx ** 2 + dy ** 2)
                      sum += val * weight
                      weightsum += weight
                      counter += 1
          
          array = None
          
          if counter > 0:
              return sum / weightsum
          else:
              return nodata