Resample 1km grid to 25km grid¶

In [1]:

import fsspec
import rioxarray as rxr
import hvplot.xarray

In [2]:

fs_osn = fsspec.filesystem('s3', anon=True, client_kwargs={'endpoint_url': 'https://mghp.osn.xsede.org'})

In [3]:

fs_osn.ls('s3://rsignellbucket1/testing/ndvi/')

Out[3]:

['rsignellbucket1/testing/ndvi/2018062.1_km_VIIRS_NDVI_cog.tif']

In [4]:

public_url ='https://mghp.osn.xsede.org/rsignellbucket1/testing/ndvi/2018062.1_km_VIIRS_NDVI_cog.tif'

In [5]:

da_1km = rxr.open_rasterio(public_url, masked=True).squeeze('band', drop=True)

In [6]:

da_1km

Out[6]:

<xarray.DataArray (y: 14505, x: 37298)>
[541007490 values with dtype=float32]
Coordinates:
  * x            (x) float64 -180.0 -180.0 -180.0 -180.0 ... 180.0 180.0 180.0
  * y            (y) float64 80.0 79.99 79.98 79.97 ... -59.98 -59.99 -60.0
    spatial_ref  int64 0
Attributes:
    scale_factor:  1.0
    add_offset:    0.0

In [7]:

%%time
da_1km = da_1km.load()

CPU times: user 9.18 s, sys: 3.9 s, total: 13.1 s
Wall time: 22.9 s

In [8]:

da_1km.nbytes/1e9   # GB 

Out[8]:

2.16402996

In [9]:

res = 25

In [10]:

%%time
da_25km = da_1km.coarsen(x=res, y=res, boundary='pad').mean()

CPU times: user 3.21 s, sys: 1.3 s, total: 4.52 s
Wall time: 4.51 s

In [11]:

da_25km.spatial_ref.attrs

Out[11]:

{'crs_wkt': 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]',
 'semi_major_axis': 6378137.0,
 'semi_minor_axis': 6356752.314245179,
 'inverse_flattening': 298.257223563,
 'reference_ellipsoid_name': 'WGS 84',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'WGS 84',
 'grid_mapping_name': 'latitude_longitude',
 'spatial_ref': 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]',
 'GeoTransform': '-179.99999550299998 0.0096519999 0.0 79.99999999300007 0.0 -0.0096519999'}

In [12]:

da_25km.hvplot(x='x', y='y', rasterize=True, geo=True, tiles='OSM', cmap='turbo', clim=(-1000,8000))

WARNING:param.Image02918: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.
WARNING:param.Image02918: Image dimension(s) x and y are not evenly sampled to relative tolerance of 0.001. Please use the QuadMesh element for irregularly sampled data or set a higher tolerance on hv.config.image_rtol or the rtol parameter in the Image constructor.

Out[12]: