docs: def union(rasters):
      """
      Union of rasters
  
      Usage:
          union(rasters)
      where:
          rasters is a list of GeoRaster objects
      """
      if sum([rasters[0].x_cell_size == i.x_cell_size for i in rasters]) == len(rasters) \
         and sum([rasters[0].y_cell_size == i.y_cell_size for i in rasters]) == len(rasters)\
         and sum([rasters[0].projection.ExportToProj4() == i.projection.ExportToProj4() for i in rasters]) == len(rasters):
          if sum([rasters[0].nodata_value == i.nodata_value for i in rasters]) == len(rasters):
              ndv = rasters[0].nodata_value
          else:
              ndv = np.nan
          if ndv == None:
              ndv = np.nan
          if sum([rasters[0].datatype == i.datatype for i in rasters]) == len(rasters):
              datatype = rasters[0].datatype
          else:
              datatype = None
          projection = rasters[0].projection
          lonmin = min([i.xmin for i in rasters])
          lonmax = max([i.xmax for i in rasters])
          latmin = min([i.ymin for i in rasters])
          latmax = max([i.ymax for i in rasters])
          shape = (np.abs(np.floor((latmax-latmin)/rasters[0].y_cell_size)).astype(int), np.floor((lonmax-lonmin)/rasters[0].x_cell_size).astype(int))
          out = ndv*np.ones(shape)
          outmask = np.ones(shape).astype(bool)
          for i in rasters:
              (row, col) = map_pixel(i.xmin, i.ymax, rasters[0].x_cell_size, rasters[0].y_cell_size, lonmin, latmax)
              out[row:row+i.shape[0], col:col+i.shape[1]] = np.where(i.raster.data != i.nodata_value, i.raster.data,\
                                                           out[row:row+i.shape[0], col:col+i.shape[1]])
              outmask[row:row+i.shape[0], col:col+i.shape[1]] = np.where(i.raster.mask == False, False,\
                                                           outmask[row:row+i.shape[0], col:col+i.shape[1]])
          out = np.ma.masked_array(out, mask=outmask, fill_value=ndv)
          return GeoRaster(out, (lonmin, rasters[0].x_cell_size, 0.0, latmax, 0.0, rasters[0].y_cell_size), nodata_value=ndv, projection=projection, datatype=datatype)
      else:
          raise RasterGeoError('Rasters need to have same pixel sizes. Use the aggregate or dissolve functions to generate correct GeoRasters')