Odd cubing behavoir with extract

[mmann1123]

I am again seeing some oddities with extract. Wondering if you have any ideas for problem solving it, just not sure where to start. I have a large vector dataset of 70k administrative units in Ethiopia, trying to get monthly precip means for each unit. Precip is from CHIRPs and is 0.05 dd res.

Example image, resampled to 0.01dd.

Open and extract code:

with gw.config.update(ref_bounds=bounds, ref_res=(0.01, 0.01)):
    with gw.open(
        f_list,
        band_names=["ppt"],
        time_names=dates,
        nodata=-9999,
        resampling="bilinear",
    ) as ds:
        print(ds)
        df = ds.gw.extract(
            aoi=eas,
            all_touched=True,
            band_names=ds.band.values.tolist(),
            time_names=ds.time.values.tolist(),
            n_threads=4,  # n_jobs creates memory error that is uncaught
            verbose=2,
        )
        print(df.head)

        ag = df.groupby(by=["id"])   # ive also tried aggregating by the administrative code, no change
            .agg("mean")
        )

1st issue, without resampling to a higher resolution say 0.01 dd, 30k+ administrative units return no data, even with 'all touched' as true.

2nd issue, if I do resample, my outputted admin units take on a really weird block pattern. I have also tried the groupby based on the unique code for the admin unit, although id seems to correspond to the index on the features. This is the same for resample = bilinear and nearest.

@jgrss

[jgrss]

@mmann1123 can you send me small versions of this data to test?

[mmann1123]

@jgrss I just shared the dropbox folders with you. Keep in mind also that I am implementing my fix to issue #75, which I am 95% sure is unrelated to this issue. Let me know if you want a more compact example to work with.

[jgrss]

Hey, I don't have an active Dropbox account. Is there another way you can send the data? Are they small enough to email?

[jgrss]

Hi @mmann1123 if your vector data are polygons geowombat.extract will do:

df = gpd.overlay(df,
                 gpd.GeoDataFrame(data=[0],
                                  geometry=[data.gw.geometry],
                                  crs=df_crs),
                 how='intersection').drop(columns=[0])

I tried to clip the shapefile that you sent by approximating the image bounds in your graphic above. When I did that I received a geopandas error. Are you able to clip the shapefile to a smaller bounding box by gpd.overlay or by:

gpd.clip(esa.to_crs('epsg:4326'), bounds)