title: Calculating spatial averages for each country after spatial join
body: <p>Hello I am using the following code at the bottom to extract countries from coordinates. Please see the following url which provides a more detailed explanation of the code: <a href="https://stackoverflow.com/questions/69979666/extracting-countries-from-netcdf-data-using-geopandas">Extracting countries from NetCDF data using geopandas</a>.</p>
  <p>My main variable/value is the monthly mean pdsi value from: <a href="https://psl.noaa.gov/data/gridded/data.pdsi.html" rel="nofollow noreferrer">https://psl.noaa.gov/data/gridded/data.pdsi.html</a>. The image below represents a portion of the visualization created by the code below. The shaded squares represent the spatial regions of pdsi values, which is overlapping a shapefile of the world.</p>
  <p><a href="https://i.stack.imgur.com/yLtu3.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/yLtu3.png" alt="Country example" /></a></p>
  <p>From the image of Belgium, you can see that the 4 squares that touch the land area of Belgium are also touching other countries. If I attribute the base values to the Belgium, I believe this overestimates the mean pdsi values. Especially when considering the bottom two squares barely touch Belgium, the weight of these values when calculating the mean should be significantly lower. Thus, is there a way to incorporate some sort of weighted average where the area of each square within a country can be used as the weight to adjust each pdsi value? Additionally, I would like to standardize this process not only for Belgium, but for all countries as well.</p>
  <p>Any help would be greatly appreciated!</p>
  <pre><code>import geopandas as gpd
  import numpy as np
  import plotly.express as px
  import requests
  from pathlib import Path
  from zipfile import ZipFile
  import urllib
  import shapely.geometry
  import xarray as xr
  
  # download NetCDF data...
  # fmt: off
  url = &quot;https://psl.noaa.gov/repository/entry/get/pdsi.mon.mean.selfcalibrated.nc?entryid=synth%3Ae570c8f9-ec09-4e89-93b4-babd5651e7a9%3AL2RhaV9wZHNpL3Bkc2kubW9uLm1lYW4uc2VsZmNhbGlicmF0ZWQubmM%3D&quot;
  f = Path.cwd().joinpath(Path(urllib.parse.urlparse(url).path).name)
  # fmt: on
  
  if not f.exists():
      r = requests.get(url, stream=True, headers={&quot;User-Agent&quot;: &quot;XY&quot;})
      with open(f, &quot;wb&quot;) as fd:
          for chunk in r.iter_content(chunk_size=128):
              fd.write(chunk)
  ds = xr.open_dataset(f)
  pdsi = ds.to_dataframe()
  pdsi = pdsi.reset_index().dropna()  # don't care about places in oceans...
  
  # use subset for testing... last 5 times...
  pdsim = pdsi.loc[pdsi[&quot;time&quot;].isin(pdsi.groupby(&quot;time&quot;).size().index[-5:])]
  
  # create geopandas dataframe
  gdf = gpd.GeoDataFrame(
      pdsim, geometry=pdsim.loc[:, [&quot;lon&quot;, &quot;lat&quot;]].apply(shapely.geometry.Point, axis=1)
  )
  
  # make sure that data supports using a buffer...
  assert (
      gdf[&quot;lat&quot;].diff().loc[lambda s: s.ne(0)].mode()
      == gdf[&quot;lon&quot;].diff().loc[lambda s: s.ne(0)].mode()
  ).all()
  # how big should the square buffer be around the point??
  buffer = gdf[&quot;lat&quot;].diff().loc[lambda s: s.ne(0)].mode().values[0] / 2
  gdf[&quot;geometry&quot;] = gdf[&quot;geometry&quot;].buffer(buffer, cap_style=3)
  
  # Import shapefile from geopandas
  path_to_data = gpd.datasets.get_path(&quot;naturalearth_lowres&quot;)
  world_shp = gpd.read_file(path_to_data)
  
  # the solution... spatial join buffered polygons to countries
  # comma separate associated countries
  gdf = gdf.join(
      world_shp.sjoin(gdf.set_crs(&quot;EPSG:4326&quot;))
      .groupby(&quot;index_right&quot;)[&quot;name&quot;]
      .agg(&quot;,&quot;.join)
  )
  gdf[&quot;time_a&quot;] = gdf[&quot;time&quot;].dt.strftime(&quot;%Y-%b-%d&quot;)
  
  # simplest way to test is visualise...
  px.choropleth_mapbox(
      gdf,
      geojson=gdf.geometry,
      locations=gdf.index,
      color=&quot;pdsi&quot;,
      hover_data=[&quot;name&quot;],
      animation_frame=&quot;time_a&quot;,
      opacity=.3
  ).update_layout(
      mapbox={&quot;style&quot;: &quot;carto-positron&quot;, &quot;zoom&quot;: 1},
      margin={&quot;l&quot;: 0, &quot;r&quot;: 0, &quot;t&quot;: 0, &quot;b&quot;: 0},
  )
  </code></pre>