sdtaylor / GrassForecasts

0 stars 0 forks source link

Using daymet in MAP calculation instead of bioclim #3

Open sdtaylor opened 4 years ago

sdtaylor commented 4 years ago

lots of issues here since daymet netCDF files are in a funky CRS, here is progress

precip = xr.open_mfdataset('./data/annual_precip/*_na.nc4', 
                      combine = 'by_coords',
                      chunks={'x':1000, 'y':1000})
precip_long_term_average = precip.mean('time')

daymet_crs = pyproj.Proj('+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +a=6378137 +b=6356752.314706705 +units=m +no_defs')
#daymet_crs = pyproj.CRS.from_proj('+proj=lcc +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +a=6378137 +rf=298.257223563 +lat_1=25 +lat_2=60')
out_proj = pyproj.CRS.from_string('epsg:4326')

# a = width of a pixel, 1000
# b = row rotation (typically zero), 0
# c = x-coordinate of the upper-left corner of the upper-left pixel -4560750
# d = column rotation (typically zero), 0
# e = height of a pixel (typically negative), -1000
# f = y-coordinate of the of the upper-left corner of the upper-left pixel, 4984500
# Affine(a,b,c,d,e,f,g)
daymet_affine = Affine(1000, 0, -4560750, 0, -1000, 4984500)

our_affine = Affine(1000preci)
longitude, latitude = pyproj.transform(daymet_crs, out_proj, precip.x.values[0], precip.y.values[0])

# Coursen the 1km daymet precip just a tad to better match the 1/8 deg cmip5
# note that here it's going from a DataSet to a DataArray
precip_long_term_average = precip_long_term_average.prcp.isel(x=range(7800), y=range(8070)).coarsen({'x':10,'y':10}).mean().compute()

precip_long_term_average = precip_long_term_average.rename({'x':'longitude','y':'latitude'})

precip_long_term_average = precip_long_term_average.compute()
p = spatial_downscale(precip_long_term_average, target_array = target_array)