<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/analysing-soil-moisture-data-in-netcdf-format-with-the-ncdf4-library/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://r.prevos.net/netcdf-format/"> The Devil is in the Data</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>The <a rel="nofollow" href="http://www.unidata.ucar.edu/software/netcdf/" target="_blank">netCDF</a> format is popular in sciences that analyse sequential spatial data. It is a self-describing, machine-independent data format for creating, accessing and sharing array-oriented information. The netCDF format provides spatial time-series such as meteorological or environmental data. This article shows how to visualise and analyse this data format by reviewing soil moisture data published by the Australian Bureau of Statistics.</p>
<h2>Soil Moisture data</h2>
<p>The <a rel="nofollow" href="http://www.bom.gov.au/water/landscape" target="_blank">Australian Bureau of Meteorology</a> publishes hydrological data in both a simple map grid and in the NetCDF format. The map grid consists of a flat text file that requires a bit of data jujitsu before it can be used. The NetCDF format is much easier to deploy as it provides a three-dimensional matrix of spatial data over time.</p>
<p>We are looking at the possible relationship between sewer main blockages and deep soil moisture levels. You will need to manually download this dataset from the Bureau of Meteorology website. I have not been able to scrape the website automatically. For this analysis, I use the actual <a rel="nofollow" href="http://www.bom.gov.au/water/landscape/#/sd_pct/Actual/Day/-39.00/130.40/4/Point/Separate/0/0/2017/8/21" target="_blank">deep soil moisture </a>level, aggregated monthly in NetCDF 4 format.</p>
<p><a rel="nofollow" href="http://www.bom.gov.au/water/landscape/#/sd_pct/Actual/Day/-39.00/130.40/4/Point/Separate/0/0/2017/8/21" target="_blank"><img alt="Soil moisture data from the Australian Bureau of meteorology in netCDF format" src="https://i1.wp.com/www.bom.gov.au/water/landscape/static/images/sd_avg.png?resize=584%2C503" data-recalc-dims="1"></a></p>
<h2>Reading, Extracting and Transforming the netCDF format</h2>
<p>The ncdf4 library, developed by David W. Pierce, provides the necessary functionality to manage this data. The first step is to load the data, extract the relevant information and transform the data for visualisation and analysis. When the data is read, it essentially forms a complex list that contains the metadata and the measurements.</p>
<p>The <code>ncvar_get</code> function extracts the data from the list. The lon, lat and dates variables are the dimensions of the moisture data. The time data is stored as the number of days since 1 January 1900. The spatial coordinates are stored in decimal degrees with 0.05-decimal degree intervals. The moisture data is a three-dimensional matrix with longitue, latitude and time as dimensions. Storing this data in this way will make it very easy to use.</p>
<pre>
library(ncdf4)
# Load data
bom &lt;- nc_open(&quot;Hydroinformatics/SoilMoisture/sd_pct_Actual_month.nc&quot;)
print(bom) # Inspect the data

# Extract data
lon &lt;- ncvar_get(bom, &quot;longitude&quot;)
lat &lt;- ncvar_get(bom, &quot;latitude&quot;)
dates &lt;- as.Date(&quot;1900-01-01&quot;) + ncvar_get(bom, &quot;time&quot;)
moisture &lt;- ncvar_get(bom, &quot;sd_pct&quot;)
dimnames(moisture) &lt;- list(lon, lat, dates)
</pre>
<h2>Visualising the data</h2>
<p>The first step is to check the overall data. This first code snippet extracts a matrix from the cube for 31 July 2017 and plots it. This code pipe extracts the date for the end of July 2017 and creates a data frame which is passed to ggplot for visualisation. Although I use the Tidyverse, I still need reshape2 because the gather function does not like matrices.</p>
<pre>
library(tidyverse)
library(RColorBrewer)
library(reshape2)

d &lt;- &quot;2017-07-31&quot;
m &lt;- moisture[, , which(dates == d)] %>%
       melt(varnames = c(&quot;lon&quot;, &quot;lat&quot;)) %>%
       subset(!is.na(value))

ggplot(m, aes(x = lon, y = lat, fill = value)) + borders(&quot;world&quot;) + 
    geom_tile() + 
    scale_fill_gradientn(colors = brewer.pal(9, &quot;Blues&quot;)) + 
    labs(title = &quot;Total moisture in deep soil layer (100-500 cm)&quot;,
    subtitle = format(as.Date(d), &quot;%d %B %Y&quot;)) + 
    xlim(range(lon)) + ylim(range(lat)) + coord_fixed()
</pre>
<p><a rel="nofollow" href="http://r.prevos.net/?attachment_id=986" target="_blank"><img data-orig-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Moisture_Australia.png?fit=1248%2C869" data-orig-size="1248,869" data-medium-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Moisture_Australia.png?fit=300%2C209" data-attachment-id="986" data-image-meta="{" data-image-description="" src="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Moisture_Australia.png?resize=584%2C407" data-permalink="http://r.prevos.net/netcdf-format/moisture_australia/" data-large-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Moisture_Australia.png?fit=584%2C407" alt="Deep soil moisture: Source Bureau of Meteorology, Australia" data-image-title="Moisture_Australia" sizes="(max-width: 584px) 100vw, 584px" data-comments-opened="1" class="wp-image-986" data-recalc-dims="1"></a></p>
<p>With the ggmap package we can create a nice map of a local area.</p>
<pre>
library(ggmap)
loc &lt;- round(geocode(&quot;Bendigo&quot;) / 0.05) * 0.05 
map_tile &lt;- get_map(loc, zoom = 12, color = &quot;bw&quot;) %>% 
    ggmap()

map_tile + 
    geom_tile(data = m, aes(x = lon, y = lat, fill = value), alpha = 0.8) + 
    scale_fill_gradientn(colors = brewer.pal(9, &quot;Blues&quot;)) + 
    labs(title = &quot;Total moisture in deep soil layer (100-500 cm)&quot;,
        subtitle = format(as.Date(d), &quot;%d %B %Y&quot;))
</pre>
<p><a rel="nofollow" href="http://r.prevos.net/?attachment_id=987" target="_blank"><img data-orig-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Bendigo_July2017.png?fit=1299%2C869" data-orig-size="1299,869" data-medium-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Bendigo_July2017.png?fit=300%2C201" data-attachment-id="987" data-image-meta="{" data-image-description="" src="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Bendigo_July2017.png?resize=584%2C391" data-permalink="http://r.prevos.net/netcdf-format/bendigo_july2017/" data-large-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/08/Bendigo_July2017.png?fit=584%2C391" alt="" data-image-title="Bendigo_July2017" sizes="(max-width: 584px) 100vw, 584px" data-comments-opened="1" class="wp-image-987" data-recalc-dims="1"></a></p>
<h2>Analysing the data</h2>
<p>For my analysis, I am interested in the time series of moisture data for a specific point on the map. The previous code slices the data horizontally over time. To create a time series we can pierce through the data for a specific coordinate. The purpose of this time series is to investigate the relationship between sewer main blockages and deep soil data, which can be a topic for a future post.</p>
<pre>
mt &lt;- data.frame(date = dates, 
                 dp = moisture[as.character(loc$lon), as.character(loc$lat), ])
ggplot(mt, aes(x = date, y = dp)) + geom_line() + 
    labs(x = &quot;Month&quot;,
         y = &quot;Moisture&quot;,
         title = &quot;Total moisture in deep soil layer (100-500 cm)&quot;,
         subtitle = paste(as.character(loc), collapse = &quot;, &quot;))
</pre>
<p>The latest version of this code is available on my <a rel="nofollow" href="https://github.com/pprevos/r.prevos.net/tree/master/Hydroinformatics/SoilMoisture" target="_blank">GitHub repository</a>.</p>
<p><a rel="nofollow" href="http://r.prevos.net/?attachment_id=994" target="_blank"><img data-orig-file="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/09/Bendigo_TimeSeries.png?fit=2271%2C1168" data-orig-size="2271,1168" data-medium-file="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/09/Bendigo_TimeSeries.png?fit=300%2C154" data-attachment-id="994" data-image-meta="{" data-image-description="" src="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/09/Bendigo_TimeSeries.png?resize=584%2C300" data-permalink="http://r.prevos.net/netcdf-format/bendigo_timeseries/" data-large-file="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/09/Bendigo_TimeSeries.png?fit=584%2C301" alt="Deep soil moisture time series." data-image-title="Bendigo_TimeSeries" sizes="(max-width: 584px) 100vw, 584px" data-comments-opened="1" class="wp-image-994" data-recalc-dims="1"></a></p>
<p>The post <a rel="nofollow" href="http://r.prevos.net/netcdf-format/" target="_blank">Analysing soil moisture data in NetCDF format with the ncdf4 library</a> appeared first on <a rel="nofollow" href="http://r.prevos.net/" target="_blank">The Devil is in the Data</a>.</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://r.prevos.net/netcdf-format/"> The Devil is in the Data</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/y_kAUpjfwo8" height="1"></a></div>