In [627]:
%matplotlib inline
from __future__ import print_function, division
import seaborn as sns
sns.set(rc={'image.cmap': 'cubehelix'})
sns.set_context('poster')
import pandas as pd
import os
import glob
import pip
import netCDF4, numpy
import numpy as np
np.random.seed(0)
from numpy import linspace, meshgrid
# GIS modules
import geopandas as gp
import folium
import postgis
import imageio
from osgeo import gdal
import rasterstats
import fiona
import vincent
vincent.core.initialize_notebook()
"""
from osgeo import gdal
from osgeo import gdal_array
from osgeo import osr"""
# Plotting functions
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection
from matplotlib.patches import PathPatch
from matplotlib.colors import rgb2hex, Normalize
from matplotlib.cm import ScalarMappable
from matplotlib.colorbar import ColorbarBase
import matplotlib.colors
from matplotlib.mlab import griddata
import geopandas
from rasterio import features
from affine import Affine
import xarray as xray
import sys
sys.path.insert(0,'..')
In [3]:
# Open the Counties shape file and inspect the contents
counties_file = "/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties.shp"
counties_gp = gp.GeoDataFrame.from_file(counties_file)
counties_gp.head(2)
print(counties_gp.columns)
lts = []
lns = []
for i in range(len(counties_gp["INTPTLAT"])):
lts.append(float(counties_gp["INTPTLAT"][i]))
lns.append(float(counties_gp["INTPTLON"][i]))
#len(counties_gp["INTPTLON"]))
print(lts[0:5], lns[0:5])
In [404]:
# The NetCDF I/O files from VIC model archive
filenames = glob.glob("/home2/svimal/Data/VIC_Fluxes/*/*")
fname = filenames[0]
print(fname)
data = netCDF4.Dataset(fname, "r")
#print(data.variables)
# Also Xarray library can be used. Use Xarray preferably.
In [392]:
# Inspect the variables inside data
print("Time, Lat, Lon, Soil_liquid")
print(len(data["time"][:]), len(data["lat"][:]), len(data["lon"][:]), len(data["Soil_liquid"][:]))
data["lat"][0:5], data["lon"][0:5] # subset of lat lon data
Out[392]:
In [383]:
# data.variables["Soil_liquid"][29][2][349][309] # #data.variables["Soil_liquid"][time][lat][lon]
data_lists = []
for soil_layer in [0]:#,1,2]:
for time in range(len(data["time"][:])):
data_list = []
if time==0:
for lat in reversed(range(len(data["lat"][:]))): # reversed to get order of tiff format lat-lon correct
for lon in range(len(data["lon"][:])):
d = float(data.variables["Soil_liquid"][time][soil_layer][lat][lon])
if numpy.isnan(float(d)) == True:
data_list.append(0)
else:
data_list.append(d)
data_lists.append(data_list)
else:
pass
In [354]:
# Length of the data lists ()
len(data_lists), len(data_lists[0])
Out[354]:
In [351]:
# Make vectors for lat and lon
lats, lons = [], []
for lat in range(350):
for lon in range(310):
lats.append(lat)
lons.append(lon)
lats = list(reversed(lats))
len(lats), len(lons)
Out[351]:
In [352]:
# Interpolation surface that uses lat, lon and Z values
def grid(x, y, z, resX=100, resY=100):
"Convert 3 column data to matplotlib grid"
xi = linspace(min(x), max(x), resX)
yi = linspace(min(y), max(y), resY)
Z = griddata(x, y, z, xi, yi, interp='linear')
X, Y = meshgrid(xi, yi)
return X, Y, Z
In [353]:
# Create the interpolated surface of the VIC variables for each lat and lon and Z (soil layer example)
month, year = "Dec", "2015"
filenames = []
for i in range(len(data_lists)):
day = str(i+1) + ""
z = data_lists[i]
x = lons
y = lats
X, Y, Z = grid(x, y, z)
plt.figure(i)
plt.contourf(X, Y, Z)
plt.colorbar()
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title("Soil Liquid: " + str(day) + "-" + str(month) + "-" + str(year))
#fname = "swe_" + str(i) + ".png"
fname = "Soil_Liquid_" + str(i) + ".png"
plt.savefig(fname)
filenames.append(fname)
plt.show()
else:
pass
print(fname)
os.getcwd()
Out[353]:
In [387]:
def write_asc(fname, data_list, filename):
"""
input: netCDF file, and list containing the x,y,z data, and filename
output: filename
creates an asc raster file and clips it to california size
"""
data = netCDF4.Dataset(fname, "r")
ncols = len(data.variables["lon"][:])
nrows = len(data.variables["lat"][:])
xllcorner = float(min(data.variables["lon"][:])) + (-0.0625/2.0)
yllcorner = float(min(data.variables["lat"][:])) + (-0.0625/2.0)
cell_size = abs(data.variables["lat"][1] - data.variables["lat"][0])
#nodata_value = -9999
myArray = np.reshape(data_list, (nrows, ncols))
with open(filename, "w") as f:
f.write("NCOLS " + str(ncols) + "\n")
f.write("NROWS " + str(nrows) + "\n")
f.write("XLLCORNER " + str(xllcorner) + "\n")
f.write("YLLCORNER " + str(yllcorner) + "\n")
f.write("CELLSIZE " + str(cell_size) + "\n")
f.write("NODATA_VALUE = -9999\n")
data_2d = np.reshape(data_list, (nrows, ncols))
for line in data_2d:
f.write(str(list(line)).strip("[").strip("]")[1:].strip("\n") + "\n")
return filename
In [388]:
fname
Out[388]:
In [ ]:
f = "/home2/svimal/Github/UCLA-Hydro/MRPI/JupyterNotebooks/sm.asc"
os.chdir("/home2/svimal/Github/UCLA-Hydro/MRPI/JupyterNotebooks/")
#os.system('gdal_translate -of "GTiff" sm.asc sm.tif')
#os.system("gdalinfo sm.tif") # https://gis.stackexchange.com/questions/65998/how-can-i-use-gdal-to-batch-define-a-projection
# Assign projection
#os.system("gdalwarp -t_srs EPSG:4326 -te -124.8437500 31.1562500 -105.4687500 53.0312500 sm.tif sm_wgs84.tif")
os.system("gdalwarp -t_srs EPSG:4326 -te -124.8437500 31.1562500 -105.4687500 53.0312500 sm.asc sm_wgs84.tif")
# Clip to boundary
#os.system("gdalwarp -cutline /home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp -crop_to_cutline -dstalpha sm_wgs84.tif sm_clipped.tif")
os.system("gdalwarp -cutline /home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp -crop_to_cutline -dstalpha sm_wgs84.tif sm_clipped_2.tif")
In [356]:
# Create a GIF from the different time steps of the VIC I/O variables
images = []
for filename in filenames:
images.append(imageio.imread(filename))
# Create GIF files from the snow water equivalent
#imageio.mimsave('SWE_GIF.gif', images, duration=0.5)
imageio.mimsave('Soil_Liquid_GIF.gif', images, duration=0.5)
os.getcwd()
Out[356]:
In [357]:
# View the GIF
from IPython.display import Image
Image(url='Soil_Liquid_GIF.gif')
Out[357]:
In [360]:
# Other input and uutput fluxes from VIC: Precipitation, Runoff, Baseflow, Air Temperature:
fname = "/home2/svimal/Data/VIC_Fluxes/1920s/fluxes.1920-12.nc"
data = netCDF4.Dataset(fname, "r")
variables = ["Prec", "Evap", "Runoff", "Baseflow", "Tair", "Soil_liquid", "SWE"]
day=1
data_lists = []
filenames=[]
for i, var in enumerate(variables):
for time in range(len(data["time"][0:1])):
data_list = []
if time==0:
try:
for lat in range(len(data["lat"][:])):
for lon in range(len(data["lon"][:])):
d = float(data.variables[var][time][lat][lon])
if numpy.isnan(d) == True:
data_list.append(0)
else:
data_list.append(d)
data_lists.append(data_list)
z = data_lists[i]
x = lons
y = lats
X, Y, Z = grid(x, y, z)
plt.figure(i)
plt.contourf(X, Y, Z)
plt.colorbar()
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.title(str(var))
#fname = "swe_" + str(i) + ".png"
#fname = str(var)+ "_" + str(i) + ".png"
#plt.savefig(fname)
#filenames.append(fname)
plt.show()
except:
pass
else:
pass
In [409]:
# define your directories and file names
dir_input = '/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/'
name_in = 'CA_counties_WGS.shp'
dir_output = '/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/'
name_out = 'CA.shp'
# create a dictionary
states = {}
# open your file with geopandas
counties = gp.GeoDataFrame.from_file(dir_input + name_in)
for i in range(len(counties)):
state_id = counties.at[i, 'STATEFP']
county_geometry = counties.at[i, 'geometry']
# if the feature's state doesn't yet exist, create it and assign a list
if state_id not in states:
states[state_id] = []
# append the feature to the list of features
states[state_id].append(county_geometry)
# create a geopandas geodataframe, with columns for state and geometry
states_dissolved = gp.GeoDataFrame(columns=['state', 'geometry'], crs=counties.crs)
# iterate your dictionary
for state, county_list in states.items():
# create a geoseries from the list of features
geometry = gp.GeoSeries(county_list)
# use unary_union to join them, thus returning polygon or multi-polygon
geometry = geometry.unary_union
# set your state and geometry values
states_dissolved.set_value(state, 'state', state)
states_dissolved.set_value(state, 'geometry', geometry)
# save to file
states_dissolved.to_file(dir_output + name_out, driver="ESRI Shapefile")
In [18]:
# 3D Globe View
plt.figure(dpi=80)
map = Basemap(projection='ortho', lat_0=30, lon_0=-90)
#llcrnrx=min(data.variables["lon"][:]),llcrnry=min(data.variables["lat"][:]),urcrnrx=max(data.variables["lon"][:]),urcrnry=max(data.variables["lat"][:])
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='green',lake_color='blue')
map.readshapefile("/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS", "counties")
#map.drawcoastlines()
x, y = map(-118, 38)
#x2, y2 = (-90, 10)
#plt.annotate('Barcelona', xy=(x, y), xycoords='data', xytext=(x2, y2), textcoords='offset points',color='r', arrowprops=dict(arrowstyle="fancy", color='g'))
x2, y2 = map(-70, 30)
plt.annotate('California \n Counties \n Shapefile', xy=(x, y), xycoords='data', xytext=(x2, y2), textcoords='data',
arrowprops=dict(arrowstyle="->"))
map.contourf(data=Z,x=X, y=Y)
plt.show()
In [412]:
# Water Bodies on the map
# California Water Bodies
water = 'lightskyblue'
earth = 'cornsilk'
juneau_lon, juneau_lat = -124.4167, 38.3
fig, ax1 = plt.subplots(figsize=(10, 10))
mm = Basemap(llcrnrlon=min(data.variables["lon"][:]), llcrnrlat=min(data.variables["lat"][:])+1,
urcrnrlon=max(data.variables["lon"][:])-8.5, urcrnrlat=max(data.variables["lat"][:])-10.5)
shp_info = mm.readshapefile(shp1[:-4],'Counties',drawbounds=True)
coast = mm.drawcoastlines()
rivers = mm.drawrivers(color=water, linewidth=1.5)
continents = mm.fillcontinents(color=earth,lake_color=water)
bound= mm.drawmapboundary(fill_color=water)
countries = mm.drawcountries()
merid = mm.drawmeridians(np.arange(-180, 180, 2), labels=[False, False, False, True])
parall = mm.drawparallels(np.arange(0, 80), labels=[True, True, False, False])
#x, y = mm(juneau_lon, juneau_lat)
#juneau = mm.scatter(x, y, 80, label="Juneau", color='red', zorder=10)
plt.title("Map of California Water Bodies")
Out[412]:
In [51]:
# Open Shapefile
shp = counties
print(shp.columns)
# Look at one of the polygon and shapefile attributes by index
shp.geometry[1]
Out[51]:
In [68]:
# View the shapefile attribute table
shp[0:10]
Out[68]:
In [378]:
# Raster Stats from the county shapefiles and the Soil Moisture Layer
shp1 = '/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp'
ras = '/home2/svimal/Github/UCLA-Hydro/MRPI/JupyterNotebooks/sm3.asc'
#ras = '/home2/svimal/Github/UCLA-Hydro/MRPI/Images/sm.asc'
stats = rasterstats.zonal_stats(shp1, ras, stats="min max mean median majority sum")
#print(stats)
# Medians
medians = []
for i in range(len(stats)):
medians.append(stats[i]["median"])
plt.ylabel("Median value of Soil Moisture"); plt.xlabel("County ID")
plt.plot(medians, label="Median")
# Means
means = []
for i in range(len(stats)):means.append(stats[i]["mean"])
plt.ylabel("Mean value of Soil Moisture"); plt.xlabel("County ID")
plt.plot(means, label="Mean")
# Maxs
maxs = [];
for i in range(len(stats)):maxs.append(stats[i]["max"])
# Mins
mins = []
for i in range(len(stats)):mins.append(stats[i]["min"])
plt.plot(mins, label="Min")
# Plot properties
plt.title("Min, Mean, Median and Max values of Soil Moisture for all counties"); plt.ylabel("County ID")
plt.plot(maxs, label="Max")
plt.legend()
Out[378]:
In [381]:
# List of Properties of the shapefile that can be queried
len(medians)
shp["NAME"][i]
# Searching for attributes
for i in range(len(shp)):
if shp.NAME[i] == "Sacramento":
print("found it!")
# Dictionary of the medians of soil moisture value in each county
medians_dict = {}
for i in range(len(shp)):
medians_dict[str(shp["NAME"][i])] = medians[i]
In [62]:
# Dictionary of Counties and their corresponding median values
print(medians_dict)
In [65]:
# Map of median values of soil moisture by counties
m = Basemap(llcrnrlon=min(data.variables["lon"][:]), llcrnrlat=min(data.variables["lat"][:])+1,
urcrnrlon=max(data.variables["lon"][:])-8.5, urcrnrlat=max(data.variables["lat"][:])-10.5)
ax = plt.gca()
fig = plt.gcf()
shp_info = m.readshapefile(shp1[:-4],'Counties',drawbounds=True)
nodata_color = "darkorange"
colors={}
County_names=[]
patches = []
cmap = plt.cm.OrRd #binary, summer https://matplotlib.org/users/colormaps.html
vmin = min(medians); vmax = max(medians)
norm = Normalize(vmin=vmin, vmax=vmax)
# color mapper to covert values to colors
mapper = ScalarMappable(norm=norm, cmap=cmap)
for shapedict in m.Counties_info:
County_name = shapedict['NAMELSAD']
County_name = County_name[:-7]
if County_name in medians_dict:
value = medians_dict[County_name]
colors[County_name] = mapper.to_rgba(value)
County_names.append(County_name)
else:
County_names.append(County_name)
colors[County_name] = nodata_color
for nshape,seg in enumerate(m.Counties):
color = rgb2hex(colors[County_names[nshape]])
poly = Polygon(seg,facecolor=color,edgecolor=color)
if (colors[County_names[nshape]] == nodata_color):
p_no = poly
ax.add_patch(poly)
plt.title('Median Values of Soil Moisture per County')
# put legend for no data states
#if p_no is not None:
# plt.legend((p_no,), ('No data',))
# construct custom colorbar
cax = fig.add_axes([0.27, 0.1, 0.5, 0.05]) # posititon
cb = ColorbarBase(cax,cmap=cmap,norm=norm, orientation='horizontal')
cb.ax.set_xlabel('Soil Moisture')
plt.show()
In [ ]:
In [750]:
def transform_from_latlon(lat, lon):
lat = np.asarray(lat)
lon = np.asarray(lon)
trans = Affine.translation(lon[0], lat[0])
scale = Affine.scale(lon[1] - lon[0], lat[1] - lat[0])
return trans * scale
def rasterize(shapes, coords, latitude='latitude', longitude='longitude',
fill=np.nan, **kwargs):
"""Rasterize a list of (geometry, fill_value) tuples onto the given
xray coordinates. This only works for 1d latitude and longitude
arrays.
"""
transform = transform_from_latlon(coords[latitude], coords[longitude])
out_shape = (len(coords[latitude]), len(coords[longitude]))
raster = features.rasterize(shapes, out_shape=out_shape,
fill=fill, transform=transform,
dtype=float, **kwargs)
spatial_coords = {latitude: coords[latitude], longitude: coords[longitude]}
return xray.DataArray(raster, coords=spatial_coords, dims=(latitude, longitude))
# this shapefile is from natural earth data
# http://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-1-states-provinces/
counties = geopandas.read_file("/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp")
county_ids = {k: i for i, k in enumerate(counties.NAME)}
shapes = zip(counties.geometry, range(len(counties)))
ds = data = xray.open_dataset("/home2/svimal/Data/VIC_Fluxes/1920s/fluxes.1924-11.nc")
# xray.open_mfdataset("/home2/svimal/Data/VIC_Fluxes/1920s/fluxes.????-??.nc") #
ds["counties"] = rasterize(shapes, ds.coords, longitude="lon", latitude="lat")
# Plot of all counties
(ds.counties.sel(lat=slice(32.5, 42.5), lon=slice(-125.2, -114.0)).plot())
plt.title("California")
plt.show()
print("All time steps in a month for Orange County")
# Make plots
(ds.Tair
#.isel(time=slice(30))
.where(ds.counties == county_ids['Orange'])
#.sel(lat=slice(33.4, 34.0), lon=slice(-118.2, -117.4))
.plot.imshow(col='time', col_wrap=4))
plt.title("All time steps in a month for Orange County")
plt.show()
# Plot of time series
(ds.Tair
#.isel(time=slice(30))
.where(ds.counties == county_ids['Orange'])
.sel(lat=slice(32, 43), lon=slice(-122, -112))
.mean(['lat', 'lon'])
.plot())
plt.title("Time series of mean for Orange County")
plt.show()
series = (ds.Tair
#.isel(time=slice(30))
.where(ds.counties == county_ids['Orange'])
.mean(['lat', 'lon'])).sel(lat=slice(32, 43), lon=slice(-122, -112))
In [ ]:
In [748]:
(ds.Tair.where(ds.counties == county_ids['Sonoma']).mean(['lat', 'lon']))
Out[748]:
In [ ]:
In [334]:
# Use GDAL functions to convert .asc to .tif, set projection, clip raster with polygon
f = "/home2/svimal/Github/UCLA-Hydro/MRPI/JupyterNotebooks/sm.asc"
os.chdir("/home2/svimal/Github/UCLA-Hydro/MRPI/JupyterNotebooks/")
#os.system('gdal_translate -of "GTiff" sm.asc sm.tif')
#os.system("gdalinfo sm.tif") # https://gis.stackexchange.com/questions/65998/how-can-i-use-gdal-to-batch-define-a-projection
# Assign projection
#os.system("gdalwarp -t_srs EPSG:4326 -te -124.8437500 31.1562500 -105.4687500 53.0312500 sm.tif sm_wgs84.tif")
os.system("gdalwarp -t_srs EPSG:4326 -te -124.8437500 31.1562500 -105.4687500 53.0312500 sm.asc sm_wgs84.tif")
# Clip to boundary
#os.system("gdalwarp -cutline /home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp -crop_to_cutline -dstalpha sm_wgs84.tif sm_clipped.tif")
os.system("gdalwarp -cutline /home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp -crop_to_cutline -dstalpha sm_wgs84.tif sm_clipped_2.tif")
Out[334]:
In [746]:
# Data aggregation
# Path to folders
forcings = glob.glob("/home2/svimal/Data/VIC_Forcing/*/*.nc") # 2 datasets because its divided into cali1 and cali2
fluxes = glob.glob("/home2/svimal/Data/VIC_Fluxes/*/*.nc") # Monthly for 86 years
# No. of files each
len(forcings), len(fluxes) # Forcing files are stored as two parts for california
# Start with fluxes
# Select all the files of one year and work on that.
years_sorted=[]
for year in range(1920,2015+1):
#print(year)
for flux_file in fluxes:
year_c = str(flux_file.split(".")[1].split("-")[0])
#print(year_c)
if str(year)==year_c:
years_sorted.append(flux_file)
year1 = years_sorted[0:12]
year1_sorted = sorted(year1, key=lambda x: x.split(".")[1].split("-")[1])
counties = geopandas.read_file("/home2/svimal/Github/CA_drought/data/Spatial/CA_counties/CA_counties_WGS.shp")
county_ids = {k: i for i, k in enumerate(counties.NAME)}
shapes = zip(counties.geometry, range(len(counties)))
for fname in year1_sorted:
#data = netCDF4.Dataset(fname, "r")
ds = data = xray.open_dataset(fname)
# xray.open_mfdataset("/home2/svimal/Data/VIC_Fluxes/1920s/fluxes.????-??.nc") #
ds["counties"] = rasterize(shapes, ds.coords, longitude="lon", latitude="lat")
#print(ds)
ds.time
In [ ]:
In [744]:
# Leap year values are in total 18 in the time period - can this be correct?
df = pd.DataFrame.from_dict(data, orient='index')
for i in range(df.shape[0]):
if numpy.isnan(float(df.sort_index().T.iloc[[53], [i]].as_matrix())):
pass
else:
print(float(df.sort_index().T.iloc[[53], [i]].as_matrix()))
In [ ]:
index = range(len(ds.time))
stacked = vincent.StackedArea(data, iter_idx='index')
stacked.axis_titles(x='Days', y='Mean Value of Variable')
stacked.legend(title='Flux Variables')
stacked.colors(brew='Spectral')
stacked
We want something like one file for each week (1-53) for each year (1920-2015) which has the format:
File name will be, e.g. 2015_week1.csv (following Karina's format in /Github/CA_drought/data/processed_data/caseCounts_westnile/weekly/2015/)
So the hydrology fluxes dataset would be (for means):
County, Precipitation, Evaporation, Surface Runoff, Baseflow, Air Temperature, Soil Liquid (layers), Snow Water Equivalent
Alameda,0,0,0,0,0,0,0
Butte,0,0,0,0,0,0,0
Calaveras,0,0,0,0,0,0,0
A similar file will be made for max, min, median as well
In [670]:
ds = ds.resample("W", "time", how="mean")
variables = ["Prec", "Evap", "Runoff", "Baseflow", "Tair", "SWE"]
index = range(len(ds.time))
data = {'index': index}
for v in variables:
series = eval("(ds." + v + ".where(ds.counties == county_ids['Orange']).mean(['lat', 'lon']))")
#.isel(time=slice(30))#.sel(lat=slice(32, 43), lon=slice(-122, -112)).
series = list(np.array(series))
data[v] = series
stacked = vincent.StackedArea(data, iter_idx='index')
stacked.axis_titles(x='Days', y='Mean Value of Variable')
stacked.legend(title='Flux Variables')
stacked.colors(brew='Spectral')
stacked
Out[670]:
In [476]:
def getXY(pt):
return (pt.x, pt.y)
centroidseries = counties['geometry'].centroid
In [755]:
series = (ds.Tair.isel(time=slice(30)).where(ds.counties == county_ids['Orange'])
#.sel(lat=slice(32, 43), lon=slice(-122, -112))
.mean(['lat', 'lon']))
buoy_map = folium.Map(location=[33.6835811,-117.9551712,10.5], zoom_start=5,tiles='Stamen Terrain')
for i, name in enumerate(counties["NAME"]):
centroid = centroidseries[i]
variables = ["Prec", "Evap", "Runoff", "Baseflow", "Tair", "SWE"]
index = range(1,31)
data = {'index': index}
for v in variables:
series = eval('(ds.' + v + '.isel(time=slice(30)).where(ds.counties == county_ids["' + name + '"]).mean(["lat", "lon"]))')
# .sel(lat=slice(32, 43), lon=slice(-122, -112))
series = list(np.array(series))
data[v] = series
stacked = vincent.StackedArea(data, iter_idx='index', width=400, height=250)
stacked.axis_titles(x='Days', y='Mean Value of Variable')
stacked.legend(title='Flux Variables for '+name + ' County')
stacked.colors(brew='Spectral')
stacked.to_json("vis" + str(i) + ".json")
"""
series = (ds.Tair.isel(time=slice(30)).where(ds.counties == county_ids[name])
.sel(lat=slice(32, 43), lon=slice(-122, -112)).mean(['lat', 'lon']))
plot_data = list(np.array(series))
#print(i, name, centroid)
#print(list(centroid.xy)[0][0] , list(centroid.xy)[1][0])
vis1 = vincent.Line(plot_data, width=400, height=250)
vis1.axis_titles(x='Day', y='Mean Values')
vis1.legend(title=name)
vis1.to_json("vis" + str(i) + ".json")
"""
popup1 = folium.Popup(max_width=600).add_child(folium.Vega(json.load(open("vis" + str(i) + ".json")), width=600, height=300))
folium.Marker([list(centroid.xy)[1][0], list(centroid.xy)[0][0]], popup=popup1).add_to(buoy_map)
buoy_map
Out[755]:
In [597]:
interactive_map = buoy_map
interactive_map.save('interactive_map2.html')
In [ ]: