Calculate a Latitudinal Gradient
Calculate a Latitudinal Gradient
Jose Lopez-Collado
April, 2016
From time to time, it is necessary to compute a statistic (mean, total, min, etc.) along a gradient in a raster file. By gradient, I mean a profile of a statistic across a list of vector values. Here, I present a simple example how to do just that to compute a latitudinal gradient by rows in a matrix.
The asc file (as1fin.asc) contain the mean annual temperatures for the state of Veracruz.
Notice: The example file is in the File Box page
The steps to calculate the gradient are:
1. Write the function to compute the statistic, taking into account null values.
2. Load the asc file, strip the header rows. Before, check what the null value and the starting coordinate values are.
3. Apply the function.
4. Generate the coordinate values.
5. Join the gradient with the coordinate values.
6. Display.
Here we go:
Write the Function
In this example, we want to compute the mean of the temperatures; to be more generic though, we will define a custom function as a mask, so we can change the function without modifying the module that explores the row/column of values:
(*
Define a custom function, the mean in this case.
If you need to compute another function, based on a vector of values, then assign that function to myF, for example, myF[v_]:= Min[v] would compute the minimum value
*)
myF[v_] := Mean[v];
(*
Write a module to select valid values from a vector
The arguments are the vector of values x_, the null value nullv_ and a value to signal that the vector is empty: flag_
Within the module, a list of valid values is selected, then list is checked to see if is not empty, so the statistic is computed; if empty, set the flag value. Check that we use N[] to force Mathematica to do the computations, otherwise it will display symbolic values if present.
*)
gFunction[x_, nullv_, flag_] := Module[{v1, avg},
v1 = Select[x, # > nullv &];
If[Dimensions[v1][[1]] > 0, N[myF[v1]], flag]
]
Read the ASC File
(* Read the file and drop the header lines, the file should be in the working directory *)
matrix = Drop[Import["asc1fin.asc", "Table"], 6];
Apply the Function
(*
Apply the function using Map
This mapping applies the gFunction to each of the rows within the matrix, if you need to compute the values by column, then you need to transpose the matrix, Transpose[matrix]
*)
gradient = gFunction[#, -9999, -1] & /@ matrix;
Generate the coordinate values
(*
Here, we establish the cell size, obtained from the asc file, and assign it to delta, then compute the latitude values by summing the interval at each row, the number of rows is also obtained from the original data matrix
*)
(* get the number of rows, 641 in this case*)
Dimensions[matrix]
{641, 610}
(*Compute a list of integers to use to compute the latitude interval, delta is the cell size *)
k = Table[i, {i, 0, 640}];
delta = 0.0083333337679505;
(* Now compute the vector of latitude values; 17.1 is the starting latitude value obtained from the original asc file *)
v = 17.133336679271 + k*delta;
(*
Join the Gradient with the Coordinate Values
Use thread to make the list of latitude and gradient pairs, reverse the list to start from the bottom value
*)
xy = Thread[{v, Reverse[gradient]}];
Display the Gradient
(*
This is the rewarding step, where we see the analysis
*)
ListLinePlot[xy, BaseStyle -> {FontSize -> 16, FontFamily -> "Helvetica"}, Filling -> Axis, AxesOrigin -> {17.15, 16.2}, PlotStyle -> {Thick, Black}, FillingStyle -> {LightOrange}, AxesLabel -> {"LATITUDE", "TEMPERATURE"}]
The drop in temperature in the middle is because at these latitudes a cold zone is present in Veracruz, where the Orizaba volcano is present.
A nice alternate plot is a smoothed histogram:
SmoothHistogram3D[xy, ColorFunction -> "NeonColors", AxesLabel -> {Style["LATITUDE", 16, Bold], Style["TEMPERATURE", 16, Bold], Style["DENSITY", 16, Bold]}, BaseStyle -> {FontFamily -> "Arial", FontSize -> 14}, ImageSize -> 650]
The SmoothHistogram3D generates a 2D kernel distribution; in other words, it represents the mathematical function of the frequency of counts. Therefore, the peaks represent the higher occurrence, here, one near 17 and the other close to 21 °C.
(* You can export the gradient list as csv for further analysis *)
Export["gradient.csv", xy]
Key words: latitudinal gradient, Mathematica, species distribution, raster file