Togo Models for Temperatures and Rainfall (proposed)

Andy Long
Spring, 2018

Setting Up:

Importing Data

datahome="/home/andy/classes/2018Spring/MathModeling/projects/Final/Model/";
datahome="http://ceadserv1.nku.edu/longa/classes/mat375/finalproject/Data/v2/" ; (* web *)
(* I had to remove missing values from the V2 file *)
file=StringJoin[datahome,"All_Clean_v2_missing_removed.xlsx"];
list=Import[file];
city="Togo"

list[[1,1,All]]

rainfall=Rest[Flatten[list[[1,All,4]]]];
location=Rest[Flatten[list[[1,All,1]]]];
longitude=Rest[Flatten[list[[1,All,6]]]];
latitude=Rest[Flatten[list[[1,All,7]]]];
elevation=Rest[Flatten[list[[1,All,8]]]];
decYear=Rest[Flatten[list[[1,All,10]]]];
(* Note that we're shifting time to 0 at the beginning of year 1961, the start of our study. This changes constant terms in models, essentially, but generally makes linear regression better conditioned: *)
decYear=decYear-1961;
minTemp=Rest[Flatten[list[[1,All,11]]]];
maxTemp=Rest[Flatten[list[[1,All,5]]]];
nvals=Length[rainfall]

These are the unique coordinates for the cities (in alphabetical order):

lats=DeleteDuplicates[latitude]
longs=DeleteDuplicates[longitude]
elevs=DeleteDuplicates[elevation]
meanElevs=Mean[elevs];

realdatamin=Table[{decYear[[i]],minTemp[[i]]},{i,1,nvals}];
realdatamax=Table[{decYear[[i]],maxTemp[[i]]},{i,1,nvals}];realdatarain=Table[{decYear[[i]],rainfall[[i]]},{i,1,nvals}];

A couple of functions we need:

OurSquared replaces R-squared when we move to a model without an intercept

OurSquared[real_,pred_]:= 
Module[{mean=Mean[real],nvals=Length[real]},
1-((∑_(i=1)^nvals)(((real[[i]]-pred[[i]]))^2))/((∑_(i=1)^nvals)(((real[[i]]-mean))^2))
]

TestResiduals does just that -- for normality.

TestResiduals[resids_]:=
Module[
{s=StandardDeviation[resids],
m=Mean[resids],𝒹,lp,table,show},
𝒹=NormalDistribution[m,s];
lp=ListPlot[resids,PlotLabel->city];
DistributionFitTest[resids];
ℋ=DistributionFitTest[resids,Automatic,"HypothesisTestData"];
show=Show[
Histogram[resids,Automatic,"ProbabilityDensity"],Plot[PDF[ℋ["FittedDistribution"],x],{x,-5,5},PlotStyle->Thick]
];
table=DistributionFitTest[resids,𝒹,{"TestDataTable","AndersonDarling"}];
GraphicsGrid[{{lp,show},{table,
ProbabilityPlot[resids,𝒹]}}]
]

Global Togo Model for Minimum Temperature

We generate the Minimum Temperature Model:

mindata=Transpose[{latitude,longitude,elevation,decYear,minTemp}];

finalmin=LinearModelFit[mindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev
},{lat,long,elev,x}]
finalmin["ParameterTable"]
finalmin["AdjustedRSquared"]

The minimum residuals are tested against a normal, and not normal:

minResiduals=finalmin["FitResiduals"];
TestResiduals[minResiduals]

minparams=finalmin["BestFitParameters"];
minTempFn[long_,lat_,elev_,x_]:=minparams.{
1,
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev}

Let’s take a look at the “long/lat trend surface”: I fix an elevation, and a year, and then let long and lat vary. I also include the coordinates of the cities, and their values (using their own particular elevations, so that one sees the impact of elevation on the model at their locations):

yearSeason=1.75;
maxTrend=Plot3D[minTempFn[x,y,meanElevs,yearSeason],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]}];
cityPoints=Table[Graphics3D[{Black,PointSize[0.05],Point[{longs[[i]],lats[[i]],minTempFn[longs[[i]],lats[[i]],elevs[[i]],yearSeason]}]}],
{i,1,10}];
Show[maxTrend,cityPoints,PlotRange->All]

Global Togo Model for Maximum Temperature

We generate the Maximum Temperature Model:

maxdata=Transpose[{latitude,longitude,elevation,decYear,maxTemp}];
maxdata[[1,All]]

finalmax=LinearModelFit[maxdata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev
},{lat,long,elev,x}]
finalmax["ParameterTable"]
finalmax["AdjustedRSquared"]

The maximum residuals are tested against a normal, and are not bad!

maxResiduals=finalmax["FitResiduals"];
TestResiduals[maxResiduals]

maxparams=finalmax["BestFitParameters"];
maxTempFn[long_,lat_,elev_,x_]:=maxparams.{1,x,Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev}

Let’s take a look at the “long/lat trend surface”: I fix an elevation, and a year, and then let long and lat vary. I also include the coordinates of the cities, and their values (using their own particular elevations, so that one sees the impact of elevation on the model at their locations):

yearSeason=1.75
maxTrend=Plot3D[maxTempFn[x,y,meanElevs,yearSeason],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]}];
cityPoints=Table[Graphics3D[{Black,PointSize[0.05],Point[{longs[[i]],lats[[i]],maxTempFn[longs[[i]],lats[[i]],elevs[[i]],yearSeason]}]}],
{i,1,10}];
Show[maxTrend,cityPoints,PlotRange->All]

Global Togo Model for Rainfall Data

Time now for the rainfall: let’s see how it depends on max and min temperatures:

raindata=Transpose[{decYear,latitude,longitude,elevation,minTemp,maxTemp,rainfall}];

finalrain=LinearModelFit[raindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt
},{x,lat,long,elev,mint,maxt}]

Since the constant term is not significant, it seems wise to drop it. I haven’t done so just yet, however. There are only modest changes to the other terms:

finalrainNoConstant=LinearModelFit[raindata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt
},{x,lat,long,elev,mint,maxt},IncludeConstantBasis->False]
finalrain["ParameterTable"]
finalrainNoConstant["ParameterTable"]

Note that when we drop the constant term, the R^2 rockets up. We don’t really believe that, however, so we consult “Our^2”, and it shows that they’re really about the same:

finalrain["AdjustedRSquared"]
finalrainNoConstant["AdjustedRSquared"]
(* Note the ugly but easy way to compute the estimates *)
OurSquared[rainfall,rainfall-finalrain["FitResiduals"]]
OurSquared[rainfall,rainfall-finalrainNoConstant["FitResiduals"]]

The residuals are not normal :(                :

rainResiduals=finalrain["FitResiduals"];
TestResiduals[rainResiduals]

The following was my kludgy way of showing the same thing:

s=StandardDeviation[rainResiduals]
m=Mean[rainResiduals]
rainHisto=Histogram[(rainResiduals-m)/s,Automatic,"ProbabilityDensity"];
rainNormal=Plot[PDF[NormalDistribution[0,1],x],{x,Min[rainResiduals],Max[rainResiduals]}];
Show[rainHisto,rainNormal]

This manipulate command allows us to examine the fit of the rainfall model across the data set, and across time. Notice that there may be some discontinuities when we switch cities (the data is all 10 cities pasted together, in alphabetic order).

Manipulate[
Module[
{
end=t+n*12,
rainmin,
rainmax,
rainavg,
yearmin,
yearmax,
p1, p2
},
yearmin=Min[Part[decYear,t;;end]];
yearmax=Max[Part[decYear,t;;end]];
rainmin=Min[Part[rainfall,t;;end]];
rainmax=Max[Part[rainfall,t;;end]];
rainavg=Mean[{rainmin,rainmax}];
rainmin=rainavg-1.5(rainavg-rainmin);
rainmax=rainavg+1.5(rainavg-rainmin);
p1=ListPlot[Table[{decYear[[i]],finalrain[decYear[[i]],latitude[[i]],longitude[[i]],elevation[[i]],minTemp[[i]],maxTemp[[i]]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{-50,500}},Joined->True,PlotLabel->location[[t]],AxesLabel->{"Year","rainerature"}];
p2=ListPlot[Table[{decYear[[i]],rainfall[[i]]},{i,t,end,1}],PlotRange->{{yearmin,yearmax},{rainmin,rainmax}},PlotStyle->PointSize[Medium],PlotLabel->location[[t]],AxesLabel->{"Year","rainerature"}];
Show[p1,p2]
],
{{t,1},1,nvals-n*12,12},
{{n,1},1,5,1}
]

Now we create a function that computes the modeled rainfall:

rainparams=finalrain["BestFitParameters"];
rainFn[long_,lat_,elev_,mint_,maxt_,x_]:=rainparams.{
1,
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev,
mint,
maxt
}

Let’s take a look at the “long/lat trend surface” for rainfall: I fix an elevation, and a year, and then let long and lat vary. I also include the coordinates of the cities, and their values (using their own particular elevations, so that one sees the impact of elevation on the model at their locations):

rainTrend=Plot3D[rainFn[x,y,meanElevs,23,28,1],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]}];
cityPoints=Table[Graphics3D[{Black,PointSize[0.05],Point[{longs[[i]],lats[[i]],rainFn[longs[[i]],lats[[i]],elevs[[i]],23,28,1]}]}],
{i,1,10}];
Show[rainTrend,cityPoints,PlotRange->All]

Now we create a function that computes the modeled rainfall, given only the long/lat coordinates of a point, and the year, using modeled min and max temperatures:

It turns out that Mathematica has a function that computes the elevation on Earth for particular latitudes and longitudes. Very cool, but very intense to compute. So if you try to plot too many points, it takes a long time....

In the end, however, all I need is the (long,lat,year) and I can predict temperatures and rainfall.

rainFnWithModel[long_,lat_,x_]:=
Module[{elevation=GeoElevationData[GeoPosition[{lat,long}],UnitSystem->"Metric"]//QuantityMagnitude},
rainparams.{
1,
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elevation,
minTempFn[long,lat,elevation,x],
maxTempFn[long,lat,elevation,x]
}
]

By choosing MaxRecursion as 0, and PlotPoints at 10, we can compute the surface fairly quickly. You might try increasing those, but it may take some time. It depends on how good your computer is....

yearChosen=39.8 (* fall, 2000: *)
rains=Table[rainFnWithModel[longs[[i]],lats[[i]],yearChosen],{i,1,10}]
cityPoints=Table[Graphics3D[{Black,PointSize[0.05],Point[{longs[[i]],lats[[i]],rains[[i]]}]}],
{i,1,10}];
rainplot=Plot3D[rainFnWithModel[x,y,yearChosen],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]},PlotPoints->10,Mesh->All,MaxRecursion->0];

Show[rainplot,cityPoints,PlotRange->All]

I was worried that Sokode wasn’t showing up, but it just kept falling below the surface.... You can see that by making it red.

Sokode=Table[Graphics3D[{Red,PointSize[0.05],Point[{longs[[i]],lats[[i]],rains[[i]]}]}],
{i,8,8}];
Show[rainplot,cityPoints,Sokode,PlotRange->All]

This manipulate command allows me to look at the modeled rainfall throughout time:

Manipulate[
Module[
{rains,cityPoints},
rains=Table[rainFnWithModel[longs[[i]],lats[[i]],yearChosen],{i,1,10}];cityPoints=Table[Graphics3D[{Black,PointSize[0.05],Point[{longs[[i]],lats[[i]],rains[[i]]}]}],
{i,1,10}];
Show[
Plot3D[rainFnWithModel[x,y,yearChosen],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]},PlotPoints->10,Mesh->All,MaxRecursion->0],
cityPoints,PlotRange->All]
],
{{yearChosen,39.5},0,139,.085}
]

How will rainfall change with time? Let’s look at rainfall at one place, but many years (100 years) apart:

p1962=Plot3D[rainFnWithModel[x,y,1],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]},PlotPoints->10,Mesh->All,MaxRecursion->0];p2062=Plot3D[rainFnWithModel[x,y,101],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]},PlotPoints->10,Mesh->All,MaxRecursion->0];
Show[p1962,p2062,PlotRange->All]

Here I look at a town’s rainfall at times 100 years apart, as we move through 100 years of time:

(* All towns look the same: *)
yearsApart=100;
indexTown=2; (* Atakpame *)
indexTown=700; (* Dapaong *)
indexTown=1300; (* Kara *)
indexTown=2000; (* Kouma *)
indexTown=2600; (* Lome *)
indexTown=3300; (* Mango *)
indexTown=3900; (* Niamtougou *)
indexTown=4500; (* Sokode *)
indexTown=5100; (* Sotouboua *)
indexTown=5800; (* Tabligbo *)
Plot[
rainFnWithModel[longitude[[indexTown]],latitude[[indexTown]],x+yearsApart]
-
rainFnWithModel[longitude[[indexTown]],latitude[[indexTown]],x]
,{x,0,100}
,PlotLabel->location[[indexTown]]
]

Here I look at a town’s rainfall value compared to its value an integer number of years into the future as we move through time:

(* All towns look the same: *)
yearSeason=0.25;
indexTown=2; (* Atakpame *)
indexTown=700; (* Dapaong *)
indexTown=1300; (* Kara *)
indexTown=2000; (* Kouma *)
indexTown=2600; (* Lome *)
indexTown=3300; (* Mango *)
indexTown=4500; (* Sokode *)
indexTown=5800; (* Tabligbo *)
indexTown=5100; (* Sotouboua *)
indexTown=3900; (* Niamtougou *)
ListPlot[Table[{1961+i+yearSeason,
rainFnWithModel[longitude[[indexTown]],latitude[[indexTown]],i+yearSeason]
-
rainFnWithModel[longitude[[indexTown]],latitude[[indexTown]],yearSeason]}
,{i,0,139}]
,PlotLabel->location[[indexTown]]
]

Here I look at a town’s difference in maxTemp and minTemp values as we move through time: as Jacob pointed out, this is a measure of variability in temperature

(* All towns look the same: *)
yearSeason=0.25;
indexTown=3900; (* Niamtougou *)
indexTown=2; (* Atakpame *)
ListPlot[Table[{1961+i+yearSeason,
maxTempFn[longitude[[indexTown]],latitude[[indexTown]],elevation[[indexTown]],i+yearSeason]
-
minTempFn[longitude[[indexTown]],latitude[[indexTown]],elevation[[indexTown]],yearSeason]}
,{i,0,139}]
,PlotLabel->location[[indexTown]]
]

Let’s take a look at Diurnal temperature ranges (DTR): differences of max and min. 
The Daily Temperature Range.

https://digitalcommons.unl.edu/cgi/viewcontent.cgi?referer=https://www.google.com/&httpsredir=1&article=1187&context=natrespapers
fulltext.pdf
Suggests that diurnal temperatures are decreasing -- variability is declining. We show the opposite, if we control for the usual suspects.... A linear regression against time suggests a decrease:

diurnalm=LinearModelFit[Transpose[{decYear,maxTemp-minTemp}],{x},{x}]

But if we add in the things we’ve been using to model, we get a positive increase:

diurnaldata=Transpose[{latitude,longitude,elevation,decYear,maxTemp-minTemp}];
diurnaldata[[1,All]]

finaldiurnal=LinearModelFit[maxdata,{
x,
Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)], 
lat,
lat*long,
(lat^2),
long,
(long^2),
elev
},{lat,long,elev,x}]
finaldiurnal["ParameterTable"]
finaldiurnal["AdjustedRSquared"]

diurnalResiduals=finaldiurnal["FitResiduals"];
TestResiduals[diurnalResiduals]

diurnalparams=finaldiurnal["BestFitParameters"];
diurnalTempFn[long_,lat_,elev_,x_]:=diurnalparams.{1,x,Sin[2*π*x],Cos[2*π*x],
Sin[2*π*x/(1/2)],Cos[2*π*x/(1/2)],
Sin[2*π*x/(1/3)],Cos[2*π*x/(1/3)],
Sin[2*π*x/(1/4)],Cos[2*π*x/(1/4)],
Sin[2*π*x/(1/5)],Cos[2*π*x/(1/5)],
Sin[2*π*x/(6)],Cos[2*π*x/(6)],
Sin[2*π*x/(12)],Cos[2*π*x/(12)],
lat,
lat*long,
(lat^2),
long,
(long^2),
elev}

Let’s take a look at the “long/lat trend surface”: I fix an elevation, and a year, and then let long and lat vary. I also include the coordinates of the cities, and their values (using their own particular elevations, so that one sees the impact of elevation on the model at their locations):

yearSeason=1.75
diurnalTrend=Plot3D[diurnalTempFn[x,y,meanElevs,yearSeason],{x,Min[longitude],Max[longitude]},{y,Min[latitude],Max[latitude]}];
cityPoints=Table[Graphics3D[{Black,PointSize[0.05],Point[{longs[[i]],lats[[i]],diurnalTempFn[longs[[i]],lats[[i]],elevs[[i]],yearSeason]}]}],
{i,1,10}];
Show[diurnalTrend,cityPoints,PlotRange->All]