In [8]:
%matplotlib inline
%qtconsole

import os
import numpy as np
import pandas as pd
import xarray as xr

import matplotlib
import matplotlib.pyplot as plt
from matplotlib.figure import Figure
from matplotlib.colors import from_levels_and_colors
from mpl_toolkits.basemap import Basemap
import plotly.plotly as py
py.sign_in('ctasich', 'fpoe1n01ek')

from datetime import datetime, timedelta
from scipy import stats
In [3]:
## Load Data

# GRACE data
nc = 'http://opendap.jpl.nasa.gov:80/opendap/GeodeticsGravity/tellus/L3/mascon/RL05/JPL/CRI/netcdf/GRCTellus.JPL.200204_201606.GLO.RL05M_1.MSCNv02CRIv02.nc'
grace = xr.open_dataset(nc)

# Well data
csv = 'https://www.hydroshare.org/django_irods/download/d3659dcf575d4090801a74d1ce096d7c/data/contents/WPDx_Well_Function_Upd_151224_xy161117.csv'
wells = pd.read_csv(csv)
In [4]:
## Preprocess data

## Wells
wells['color'] = np.where(wells['FUNC']=='Yes', '#2ECC71', '#E74C3C')

## GRACE
rmap = grace['lwe_thickness'][0,:,:]

# Extract Lat/Lon Metadata
lat_min = grace.geospatial_lat_min
lat_max = grace.geospatial_lat_max
lat_res = float(grace.geospatial_lat_resolution[0:3])

lon_min = grace.geospatial_lon_min
lon_max = grace.geospatial_lon_max
lon_res = float(grace.geospatial_lon_resolution[0:3])
In [6]:
## Plot GRACE data

# Build grid
lon_g = np.arange(lon_min,lon_max+lon_res,lon_res)
lat_g = np.arange(lat_min,lat_max+lat_res,lat_res)
x,y = np.meshgrid(lon_g[:], lat_g[:])

# Plot Fig
plt.figure(figsize=(20,10))
m = Basemap(projection='moll',llcrnrlat=-87,urcrnrlat=81,lon_0=0,\
            llcrnrlon=0,urcrnrlon=360,resolution='c')
# draw parallels and meridians.
parallels = np.arange(-89.75,89.75,30.)
# Label the meridians and parallels
m.drawparallels(parallels,labels=[False,True,True,False])
# Draw Meridians and Labels
meridians = np.arange(-180.,181.,30.)
m.drawmeridians(meridians)
m.drawmapboundary(fill_color='white')

ax = plt.gca()
masked_array = np.ma.array(rmap, mask=np.isnan(rmap))
cmap = matplotlib.cm.jet
cmap.set_bad('white',1.0)

im1 = m.pcolormesh(x,y,rmap,shading='flat',latlon=True);
im2 = m.pcolormesh(x,y,masked_array,shading='flat',latlon=True)
m.drawcoastlines();
cbar = plt.colorbar(fraction=0.046, pad=0.04)
cbar.set_label('Water thickness equivalent (cm)')
plt.title('GRACE initial measurement',size=20);
In [9]:
## Plot well data in plotly

data = [ dict(
    lat = wells.LAT_DD,
    lon = wells.LONG_DD,
    marker = dict(
        color = wells.color.tolist(),
        opacity = 0.7,
        size = 2,                
    ),
    type = 'scattergeo'
) ]

layout = dict(
    geo = dict(showland = True,
        landcolor = "rgb(212, 212, 212)",
        subunitcolor = "rgb(255, 255, 255)",
        countrycolor = "rgb(255, 255, 255)",
        showlakes = True,
        lakecolor = "rgb(255, 255, 255)",
        showsubunits = True,
        showcountries = True,
        resolution = 10,
        projection = dict(
            type = 'utm'),
        lonaxis = dict(
            showgrid = True,
            gridwidth = 0.5,
            range= [ -20, 80 ],
            dtick = 5
        ),
        lataxis = dict (
            showgrid = True,
            gridwidth = 0.5,
            range= [ -20, 40 ],
            dtick = 5
        )
    ),
    title = 'Wells from WPDx',
)
fig = { 'data':data, 'layout':layout }
py.iplot(fig, filename='wells')
Out[9]:
In [6]:
# Code that subselects regions of interest. This is for all of Africa, but will be used later to get individual time series

data = xr.open_dataset('http://opendap.jpl.nasa.gov:80/opendap/GeodeticsGravity/tellus/L3/mascon/RL05/JPL/CRI/netcdf/GRCTellus.JPL.200204_201606.GLO.RL05M_1.MSCNv02CRIv02.nc')

af = xr.concat( [data['lwe_thickness'].sel(lat=slice(-37.75,37.75)).sel(lon=slice(340.25,359.75)),
                  data['lwe_thickness'].sel(lat=slice(-37.75,37.75)).sel(lon=slice(0.25,50.75))],
                  dim='lon')

lonaf = xr.concat( [data['lon'].sel(lon=slice(340.25,359.75)),
                  data['lon'].sel(lon=slice(0.25,50.75))],
                  dim='lon')

lataf = data['lat'].sel(lat=slice(-37.75,37.75))
In [7]:
# Find nearest grid locations for all data
# lon_g and lat_g are the lons and lats of the gridded products, respectively
# nb this is only for Africa for now! Change things in the previous cell if you want to deal with the global GRACE dataset.

lon_g = lonaf
lat_g = lataf

xRes = np.median(np.diff(lon_g))
yRes = np.median(np.diff(lat_g))

# Define grid box centers
lon_c = lon_g[:-1]+xRes/2
lat_c = lat_g[:-1]+yRes/2

# Define a new metadata file that has grid coordinates for this resolution choice
wg = wells

wg.loc[:,'grid_lat'] = np.nan
wg.loc[:,'grid_lon'] = np.nan
wg.loc[:,'grace_mean'] = np.nan
wg.loc[:,'grace_std'] = np.nan
wg.loc[:,'grace_at_rpt_date'] = np.nan

## Determine grid_lat and grid_lon for every record

for index, row in wg.iterrows():
    lon_s = row[u'LONG_DD']
    lat_s = row[u'LAT_DD']
    # correct for wrapping
    if lon_s<0:
        lon_s = 360+lon_s
    glat = lat_g.values[np.argmin(np.abs(lat_c.values-lat_s))]
    glon = lon_g.values[np.argmin(np.abs(lon_c.values-lon_s))]
    wg.set_value(index,'grid_lat',glat)
    wg.set_value(index,'grid_lon',glon)

# Get all unique grid_lat and grid_lon pairs. Don't totally understand this bit of magic...
allpairs = wg[['grid_lat', 'grid_lon']].values
upairs = np.array(list(set(tuple(p) for p in allpairs)))

# GRACE at well locations. sel_points is necessary to get coordinate pairs.
wellG = data['lwe_thickness'].sel_points(lat=upairs[:,0],lon=upairs[:,1])
In [8]:
# Loop through the dataframe again and compute stats!

for index, row in wg.iterrows():
    glat = row[u'grid_lat']
    glon = row[u'grid_lon']

    # get the corresponding point
    pt = wellG[(wellG['lat']==glat).values & (wellG['lon']==glon).values].points.values
    allhere = wellG.sel(points=pt)
    wg.set_value(index,'grace_mean',np.mean(allhere.values))
    wg.set_value(index,'grace_std',np.std(allhere.values))

    dda = row[u'RPT_DATE']
    dgr = allhere['time'].values
    best_gr_ind = np.argmin(np.abs(pd.to_datetime(dgr)-pd.to_datetime(dda)))

    wg.set_value(index,'grace_at_rpt_date',np.squeeze(allhere.values)[best_gr_ind])
In [9]:
# Take a look
wg
Out[9]:
WELL_ID LAT_DD LONG_DD FUNC STATUS COD_FCN COD_QTY COD_RESRCE ADM1 ADM2 ... MGMT PAY SOURCE RPT_DATE color grid_lat grid_lon grace_mean grace_std grace_at_rpt_date
0 362092 5.982436 -8.180609 Yes Working but with problems. Well polluted|Under... 2 1 0 Grand Gedeh Tchien ... NaN No water committee WASH Liberia 21/01/2011 #2ECC71 5.75 351.75 1.420943 10.576535 -0.053641
1 362100 5.899207 -8.173315 Yes Working but with problems. Well polluted|Under... 2 1 0 Grand Gedeh Tchien ... NaN No water committee WASH Liberia 21/01/2011 #2ECC71 5.75 351.75 1.420943 10.576535 -0.053641
2 357349 5.802157 -9.645714 Yes Working but with problems. Not priming 2 1 0 Rivercess Norwein ... NaN No water committee WASH Liberia 02/02/2011 #2ECC71 5.75 350.25 1.083445 13.598019 -7.108154
3 489514 -0.541100 34.375820 No Drought|No operation in the dry season 999 0 1 Homa Bay Mbita ... 0 No payment system Engineering Sciences & Global Development 24/01/2011 #E74C3C -0.75 34.25 9.842547 12.380296 0.130227
4 357595 5.716055 -9.618187 No Broken Down System. low water table 0 1 1 Rivercess Norwein ... NaN No water committee WASH Liberia 03/02/2011 #E74C3C 5.25 350.25 0.111561 0.825407 -0.457256
5 489266 -0.730000 34.366000 No No fuel|No operation at least once a week 2 1 0 Homa Bay Ndhiwa ... Private Operator/Delegated Management Per Bucket Engineering Sciences & Global Development 18/02/2011 #E74C3C -0.75 34.25 9.842547 12.380296 0.130227
6 489625 -0.450333 34.009880 No Low yield|No operation in the dry season 2 1 1 Homa Bay Mbita ... Institutional Management No payment system Engineering Sciences & Global Development 02/10/2011 #E74C3C -0.75 33.75 9.842547 12.380296 0.130227
7 364570 5.231378 -9.141873 Yes Working but with problems. Well polluted|Under... 2 1 0 Sinoe Sanquin Dist#2 ... Community Management Yes but only in case of breakdown WASH Liberia 27/01/2011 #2ECC71 4.75 350.75 0.111561 0.825407 -0.457256
8 361779 5.225134 -8.121493 Yes Working but with problems. insufficient water 2 1 0 River Gee Karforh ... NaN No water committee WASH Liberia 18/02/2011 #2ECC71 4.75 351.75 1.420943 10.576535 -0.053641
9 361780 5.225422 -8.119787 Yes Working but with problems. insufficient water 2 1 0 River Gee Karforh ... NaN No water committee WASH Liberia 18/02/2011 #2ECC71 4.75 351.75 1.420943 10.576535 -0.053641
10 365177 5.190745 -9.043733 Yes Working but with problems. low pressure 2 1 0 Sinoe Butaw ... Community Management Yes but only in case of breakdown WASH Liberia 28/01/2011 #2ECC71 4.75 350.75 0.111561 0.825407 -0.457256
11 489624 -0.480917 34.009890 No Equipment not-function|No operation in the dry... 0 0 1 Homa Bay Mbita ... Institutional Management No payment system Engineering Sciences & Global Development 02/09/2011 #E74C3C -0.75 33.75 9.842547 12.380296 0.130227
12 359723 5.020838 -7.660045 No Broken Down System. low water 0 0 1 River Gee Tuobo ... NaN No water committee WASH Liberia 07/02/2011 #E74C3C 4.75 352.25 1.420943 10.576535 7.399117
13 365567 4.994706 -7.896630 Yes Working but with problems. Require Redigging|R... 2 0 1 Maryland Gwelekpoken ... Community Management Never WASH Liberia 29/01/2011 #2ECC71 4.75 351.75 1.420943 10.576535 -0.053641
14 359873 6.386340 -10.755959 Yes Working but with problems. low pressire 2 1 0 Montserrado Greater Monrovia ... Other No water committee WASH Liberia 07/06/2011 #2ECC71 6.25 348.75 0.111561 0.825407 0.757440
15 359934 6.379438 -10.775051 Yes Working but with problems. low pressure 2 1 0 Montserrado Greater Monrovia ... Community Management Yes but only in case of breakdown WASH Liberia 07/06/2011 #2ECC71 6.25 348.75 0.111561 0.825407 0.757440
16 363124 6.334777 -10.792682 Yes Working but with problems. water polluted 2 1 0 Montserrado Greater Monrovia ... Other No water committee WASH Liberia 24/05/2011 #2ECC71 6.25 348.75 0.111561 0.825407 -1.211578
17 362739 6.333671 -10.785141 Yes Working but with problems. water pressure away... 2 1 0 Montserrado Greater Monrovia ... Direct Government Operation No water committee WASH Liberia 23/05/2011 #2ECC71 6.25 348.75 0.111561 0.825407 -1.211578
18 366595 6.329771 -10.755218 Yes Working but with problems. low pressure 2 1 0 Montserrado Greater Monrovia ... Other No water committee WASH Liberia 31/05/2011 #2ECC71 6.25 348.75 0.111561 0.825407 -1.211578
19 366548 6.328716 -10.755411 Yes Working but with problems. Apron damaged|sticking 2 1 0 Montserrado Greater Monrovia ... Other No water committee WASH Liberia 31/05/2011 #2ECC71 6.25 348.75 0.111561 0.825407 -1.211578
20 363676 6.310975 -10.798569 Yes Working but with problems. no water 0 1 0 Montserrado Greater Monrovia ... Other No water committee WASH Liberia 25/05/2011 #2ECC71 6.25 348.75 0.111561 0.825407 -1.211578
21 364879 6.307755 -10.799409 Yes Working but with problems. dry|dry 2 1 0 Montserrado Greater Monrovia ... Direct Government Operation No water committee WASH Liberia 27/05/2011 #2ECC71 6.25 348.75 0.111561 0.825407 -1.211578
22 489544 -0.550333 34.387750 No Low yield|Normally operational 2 1 0 Homa Bay Mbita ... Community Management Per Month Engineering Sciences & Global Development 25/01/2011 #E74C3C -0.75 34.25 9.842547 12.380296 0.130227
23 489510 -0.530250 34.349920 No Low yield|No operation in the dry season 2 1 1 Homa Bay Mbita ... 0 No payment system Engineering Sciences & Global Development 24/01/2011 #E74C3C -0.75 34.25 9.842547 12.380296 0.130227
24 489698 -0.521000 34.275530 No Low yield|Normally operational 2 1 0 Homa Bay Mbita ... Community Management Per Month Engineering Sciences & Global Development 20/01/2011 #E74C3C -0.75 34.25 9.842547 12.380296 0.130227
25 489716 -0.733517 34.190250 Yes No operation in the dry season 2 0 1 Homa Bay Suba ... Community Management Per Month Engineering Sciences & Global Development 02/04/2011 #2ECC71 -0.75 33.75 9.842547 12.380296 0.130227
26 489474 -0.631300 34.459320 Yes In the dry season 2 0 1 Homa Bay NaN ... Community Management Per Month Engineering Sciences & Global Development 24/03/2011 #2ECC71 -0.75 34.25 9.842547 12.380296 -2.929331
27 359736 4.645290 -7.582768 Yes Working but with problems. Insufficient during... 2 0 1 Maryland Pleebo/Sodoken ... Community Management Yes but only in case of breakdown WASH Liberia 07/02/2011 #2ECC71 4.25 352.25 0.141984 0.967595 0.820992
28 489357 -0.650850 34.406880 Yes In the dry season 2 0 1 Homa Bay Ndhiwa ... NaN No payment system Engineering Sciences & Global Development 01/02/2011 #2ECC71 -0.75 34.25 9.842547 12.380296 7.022088
29 489358 -0.642433 34.394130 Yes In the dry season 2 0 1 Homa Bay Ndhiwa ... Direct Government Operation No payment system Engineering Sciences & Global Development 01/02/2011 #2ECC71 -0.75 34.25 9.842547 12.380296 7.022088
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
17981 490728 1.759500 32.375187 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... Community Management No payment system Lifeline 28/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17982 490093 1.820680 32.466522 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... Direct Government Operation No payment system Lifeline 27/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17983 490997 1.848909 32.513946 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... Community Management No payment system Lifeline 27/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17984 490730 2.053498 32.644738 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Kwania ... Community Management No payment system Lifeline 22/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17985 490435 1.733558 32.496951 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... NaN No payment system Lifeline 28/11/2014 #E74C3C 1.25 32.25 16.890741 19.340786 50.986263
17986 490439 1.730313 32.468445 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... NaN No payment system Lifeline 28/11/2014 #E74C3C 1.25 32.25 16.890741 19.340786 50.986263
17987 490440 1.776072 32.458147 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... NaN No payment system Lifeline 28/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17988 490434 1.733324 32.496994 No System is Broken|Seasonal Shortages|Broken Par... 0 1 1 Apac Maruzi ... NaN No payment system Lifeline 28/11/2014 #E74C3C 1.25 32.25 16.890741 19.340786 50.986263
17989 489847 1.866083 32.719903 No System is Broken|Seasonal Shortages|pipes drop... 0 0 1 Apac Kwania ... Community Management No payment system Lifeline 26/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17990 490214 1.997437 32.847491 No System is Broken|Seasonal Shortages|system is ... 0 1 1 Apac Kwania ... Community Management No payment system Lifeline 25/11/2014 #E74C3C 1.75 32.75 4.595492 8.659639 25.297725
17991 489981 1.700627 32.396367 No System is Broken|System is Down for Maintenanc... 0 0 1 Apac Maruzi ... Community Management No payment system Lifeline 27/11/2014 #E74C3C 1.25 32.25 16.890741 19.340786 50.986263
17992 489940 2.088080 32.823175 No System is Broken|System is Down for Maintenanc... 0 1 1 Apac Kwania ... Community Management 1 UGX Per Month Lifeline 24/11/2014 #E74C3C 1.75 32.75 4.595492 8.659639 25.297725
17993 490881 1.912940 32.862696 No System is Broken|System is Down for Maintenanc... 0 0 1 Apac Kwania ... Community Management 1 UGX Per Year Lifeline 26/11/2014 #E74C3C 1.75 32.75 4.595492 8.659639 25.297725
17994 489947 1.919767 32.706646 No System is Broken|System is Down for Maintenanc... 0 0 1 Apac Kwania ... Community Management 1 UGX Per Year Lifeline 25/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17995 490488 2.020943 32.369416 No System is Broken|broken down 0 1 0 Apac Maruzi ... NaN No payment system Lifeline 28/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17996 490747 1.668321 32.541246 No System is Broken|not working|Broken Parts|too ... 0 1 0 Apac Maruzi ... NaN 1 UGX Per Month Lifeline 27/11/2014 #E74C3C 1.25 32.25 16.890741 19.340786 50.986263
17997 490105 2.125865 32.540418 No System is Broken|system not functioning 0 1 0 Apac Maruzi ... NaN No payment system Lifeline 22/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
17998 490814 1.866107 32.900753 No System is Broken|system not working|Broken Par... 0 1 0 Apac Kwania ... NaN 1 UGX Per Month Lifeline 26/11/2014 #E74C3C 1.75 32.75 4.595492 8.659639 25.297725
17999 490288 1.746125 32.419225 No System is Down for Maintenance or Repair|Rati... 0 1 0 Apac Maruzi ... Community Management 5 UGX Per Month Lifeline 27/11/2014 #E74C3C 1.25 32.25 16.890741 19.340786 50.986263
18000 489919 1.851014 32.731565 No System is Down for Maintenance or Repair|Seas... 0 0 1 Apac Kwania ... Community Management No payment system Lifeline 27/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
18001 483831 59.278134 18.015202 missing missing 999 1 0 NaN NaN ... NaN NaN mWater 02/10/2014 #E74C3C 37.25 17.75 0.296950 2.991807 3.041266
18002 483375 8.675615 16.840204 missing missing 999 1 0 NaN NaN ... NaN NaN mWater 15/07/2014 #E74C3C 8.25 16.75 0.101542 13.290039 -18.746097
18003 483695 7.867364 3.961383 missing missing 999 1 0 NaN NaN ... NaN NaN mWater 19/07/2014 #E74C3C 7.75 3.75 1.421210 15.684912 22.464255
18004 484333 7.597172 5.230582 missing missing 999 1 0 NaN NaN ... NaN NaN mWater 16/06/2014 #E74C3C 7.25 4.75 1.287055 15.401714 -0.982665
18005 483321 -2.511946 32.900248 missing missing 999 1 0 NaN NaN ... NaN NaN Mwanza 03/09/2014 #E74C3C -2.75 32.75 7.437896 11.471042 22.979463
18006 483841 7.408475 4.064835 missing missing 999 1 0 NaN NaN ... NaN NaN mWater 18/12/2014 #E74C3C 7.25 3.75 1.287055 15.401714 20.469362
18007 483900 7.408560 4.064823 missing missing 999 1 0 NaN NaN ... NaN NaN mWater 18/12/2014 #E74C3C 7.25 3.75 1.287055 15.401714 20.469362
18008 490128 2.042674 32.745330 No system is down for two years|Broken Parts|pipe... 0 1 0 Apac Kwania ... NaN No payment system Lifeline 23/11/2014 #E74C3C 1.75 32.25 4.595492 8.659639 25.297725
18009 0 0.000000 0.000000 NaN Non-Functional - - Water dried 0 0 1 NaN NaN ... NaN NaN NaN NaN #E74C3C -0.25 0.25 -0.037223 0.818951 -0.652243
18010 0 0.000000 0.000000 NaN Non-Functional - - Water dried\r/\r//\r/- A lo... 0 0 1 NaN NaN ... NaN NaN NaN NaN #E74C3C -0.25 0.25 -0.037223 0.818951 -0.652243

18011 rows × 25 columns

In [10]:
# Is there any relationship between GRACE and the well data?

# 1. Compare g_mean to g_rpt. Make 2 histograms, one for working and one for not. Anything there?


# differences between mean and report time GRACE values
d_mean_rpt_yes = wg[wg['FUNC']=='Yes' ]['grace_at_rpt_date']-wg[wg['FUNC']=='Yes']['grace_mean']
d_mean_rpt_no  = wg[wg['FUNC']=='No'  ]['grace_at_rpt_date']-wg[wg['FUNC']=='No' ]['grace_mean']

bins = np.arange(-40,40)

plt.hist(d_mean_rpt_yes,bins=bins,alpha=0.5,label='Not functioning')
plt.hist(d_mean_rpt_no ,bins=bins,alpha=0.5,label='Functioning')
plt.ylabel('Number of records')
plt.title('GRACE LWE at report dates for African sites minus 2002-2016 mean',size=12)
plt.legend(loc='upper right')
plt.xlabel('Anomaly in liquid water equivalent (cm)')
plt.show()

# 2. Histograms of g_std for working and not. Any difference?

d_std_yes = wg[wg['FUNC']=='Yes' ]['grace_std']
d_std_no  = wg[wg['FUNC']=='No'  ]['grace_std']

bins = np.arange(0,30)

plt.hist(d_std_yes,bins=bins,alpha=0.5,label='Not functioning')
plt.hist(d_std_no ,bins=bins,alpha=0.5,label='Functioning')
plt.ylabel('Number of records')
plt.title('GRACE standard deviation of LWE at African sites, 2002-2016')
plt.legend(loc='upper right')
plt.xlabel('Liquid water equivalent (cm)')
plt.show()