title: Calculate mean value for each pixel of a sum of Xarray DataArrays
body: <p>I am trying to calculate a fog frequency map based on a number of geoTIFFs that I have read as Xarray DataArrays using the <code>rioxarray.open_rasterio</code> function in Python 3.10. Each &quot;pixel&quot; can have one of the following values: 1 = fog, 0 = no fog, -9999 = no data. The end goal is to calculate a new DataArray that contains the ratio of the number &quot;fog&quot; pixels to the number of pixel with either &quot;fog&quot; or &quot;no fog&quot;.</p>
  <p>For this I want to write a for-loop that creates the sum of &quot;fog&quot; and &quot;no_fog&quot; entries per pixel while excluding the &quot;no data&quot; pixels. Then it should divide the pixel values of the sum DataArray by the number of pixels that were used in the calculation of each individual sum. So, if for a single pixel there are the following values: 0, 1, 1, -9999, 0, and -9999 the loop should create a sum of 2 and divide it by 4, creating a fog frequency of 0.5 or 50%.</p>
  <p>So far, I have only been able to calculate the sum of all input DataArrays, without excluding the &quot;no data&quot; pixels using this code:</p>
  <pre><code># open all fog maps and create a list:
  folder = &quot;E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency&quot;
  list_of_maps = glob.glob(folder + '/fog_map*.tif', recursive=True)  # all files that start with &quot;fog_map&quot;
  
  # make list with all different filenames (dates) in this folder:
  maps = []  # initialize empty list for all file names
  for i in range(0, np.size(list_of_maps)):
      # files naming convention &quot;fog_map_YYYYMMDD_HHMMSS.tif&quot;:
      maps.append(list_of_maps[i].split('fog_map_')[1][0:8])
  
  # find out how many dates are in the folder:
  maps = np.unique(maps)  # remove duplicates from array
  print(maps)
  print('\ndata from {} different dates in this folder\n'.format(np.size(maps)))
  
  # create fog_sum xarray dataArray to have something to start out with and later subtract it again:
  fog_sum = rioxarray.open_rasterio(&quot;E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency/fog_map_20210601.tif&quot;)
  fog_sum_subtract = rioxarray.open_rasterio(&quot;E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency/fog_map_20210601.tif&quot;)
  
  # add all fog maps:
  for i in range(0, np.size(list_of_maps)):
      # open data sets:
      fog_map = rioxarray.open_rasterio(list_of_maps[i], engine='rasterio')
      # fog_map = fog_map.where(fog_map &gt;= 0)
      fog_sum = fog_sum + fog_map
  
  # subtract original fog map and export as geoTIFF:
  fog_sum = fog_sum - fog_sum_subtract
  fog_sum.rio.to_raster(&quot;E:/Jasper/Studium/BA_Thesis/MODIS_data/MODIS_2021_data/2021_06/fog_frequency/fog_sum.tif&quot;,
                        driver=&quot;GTiff&quot;)
  </code></pre>
  <p>I tried to exclude the &quot;no data&quot; values using <code>fog_map = fog_map.where(fog_map &gt;= 0)</code>, but this left me with a fog_sum geoTIFF, where each pixel had the value 1.79769e+308</p>
  <p>This is an example of what the output of a <code>fog_map_YYYYMMDD.tif</code> DataArray looks like, before applying the <code>fog_map.where(fog_map &gt;= 0)</code> function:</p>
  <pre><code>&lt;xarray.DataArray (band: 1, y: 412, x: 388)&gt;
  [159856 values with dtype=float32]
  Coordinates:
    * band         (band) int32 1
    * x            (x) float64 -92.49 -92.49 -92.48 ... -89.02 -89.02 -89.01
    * y            (y) float64 2.0 1.991 1.982 1.973 ... -1.674 -1.683 -1.692
      spatial_ref  int32 0
  Attributes:
      STATISTICS_MAXIMUM:        1
      STATISTICS_MEAN:           -858.62379891903
      STATISTICS_MINIMUM:        -9999
      STATISTICS_STDDEV:         2801.4551987932
      STATISTICS_VALID_PERCENT:  100
      scale_factor:              1.0
      add_offset:                0.0
  </code></pre>
  <p>And after applying the function:</p>
  <pre><code>&lt;xarray.DataArray (band: 1, y: 412, x: 388)&gt;
  array([[[ 0.,  0.,  0., ...,  0.,  0.,  0.],
          [ 0.,  0.,  0., ...,  0., nan,  0.],
          [ 0.,  0.,  0., ...,  0.,  0.,  0.],
          ...,
          [ 0.,  0.,  0., ..., nan, nan, nan],
          [ 0.,  0.,  0., ..., nan, nan, nan],
          [ 0.,  0.,  0., ..., nan, nan, nan]]], dtype=float32)
  Coordinates:
    * band         (band) int32 1
    * x            (x) float64 -92.49 -92.49 -92.48 ... -89.02 -89.02 -89.01
    * y            (y) float64 2.0 1.991 1.982 1.973 ... -1.674 -1.683 -1.692
      spatial_ref  int32 0
  Attributes:
      STATISTICS_MAXIMUM:        1
      STATISTICS_MEAN:           -858.62379891903
      STATISTICS_MINIMUM:        -9999
      STATISTICS_STDDEV:         2801.4551987932
      STATISTICS_VALID_PERCENT:  100
      scale_factor:              1.0
      add_offset:                0.0
  </code></pre>
  <p>Any help is greatly appreciated!</p>