Extract daily point precipitation

In AWAPer: Catchment Area Weighted Climate Data Anywhere in Australia

library(AWAPer, warn.conflicts = FALSE)

Make netCDF file

The first step is to create the netCDF files. Here only one netCDF file is created and only for precipitation data (as defined by the following URLS for the source data being set to NA: urlTmin, urlTmax, urlVprp, urlSolarrad) and only between the dates updateTo and updateTo. If the latter two dates were not input then data would be downloaded from 1/1/1900 to yesterday. The netCDF file contains grids of daily rainfall for all of Australia and is used below to extract data at points of interest. Often users run makeNetCDF_file once to build netCDF data files that contain all variables over the entire record length (which requires ~5GB disk storage) and then use the netCDFs grids for multiple projects, rather than re-building the netCDF for each project. Also, if makeNetCDF_file is run with the the netCDF file names pointing to existing files and updateFrom=NA then the netCDF files will be updated to yesterday.

fnames = makeNetCDF_file(ncdfFilename ='AWAP_demo.nc',
                         updateFrom=as.Date("2010-08-01","%Y-%m-%d"),
                         updateTo=as.Date("2010-10-01","%Y-%m-%d"),
                         urlTmin=NA, urlTmax=NA, urlVprp=NA, urlSolarrad=NA)
#> Starting to build both netCDF files.
#> ... Testing downloading of AWAP precip. grid
#> ... Getting grid gemoetry from file.
#> ... Deleting /home/timjp/Documents/SRC/AWAPer/vignettes/precip.20000101.grid.gz
#> ... Building AWAP netcdf file.
#> ... NetCDF data will be  extracted from  2010-08-01  to  2010-10-01
#> ... Starting to add data AWAP netcdf file.
#> Working on grid time point: 2010-08-01
#> Working on grid time point: 2010-08-02
#> Working on grid time point: 2010-08-03
#> Working on grid time point: 2010-08-04
#> Working on grid time point: 2010-08-05
#> Working on grid time point: 2010-08-06
#> Working on grid time point: 2010-08-07
#> Working on grid time point: 2010-08-08
#> Working on grid time point: 2010-08-09
#> Working on grid time point: 2010-08-10
#> Working on grid time point: 2010-08-11
#> Working on grid time point: 2010-08-12
#> Working on grid time point: 2010-08-13
#> Working on grid time point: 2010-08-14
#> Working on grid time point: 2010-08-15
#> Working on grid time point: 2010-08-16
#> Working on grid time point: 2010-08-17
#> Working on grid time point: 2010-08-18
#> Working on grid time point: 2010-08-19
#> Working on grid time point: 2010-08-20
#> Working on grid time point: 2010-08-21
#> Working on grid time point: 2010-08-22
#> Working on grid time point: 2010-08-23
#> Working on grid time point: 2010-08-24
#> Working on grid time point: 2010-08-25
#> Working on grid time point: 2010-08-26
#> Working on grid time point: 2010-08-27
#> Working on grid time point: 2010-08-28
#> Working on grid time point: 2010-08-29
#> Working on grid time point: 2010-08-30
#> Working on grid time point: 2010-08-31
#> Working on grid time point: 2010-09-01
#> Working on grid time point: 2010-09-02
#> Working on grid time point: 2010-09-03
#> Working on grid time point: 2010-09-04
#> Working on grid time point: 2010-09-05
#> Working on grid time point: 2010-09-06
#> Working on grid time point: 2010-09-07
#> Working on grid time point: 2010-09-08
#> Working on grid time point: 2010-09-09
#> Working on grid time point: 2010-09-10
#> Working on grid time point: 2010-09-11
#> Working on grid time point: 2010-09-12
#> Working on grid time point: 2010-09-13
#> Working on grid time point: 2010-09-14
#> Working on grid time point: 2010-09-15
#> Working on grid time point: 2010-09-16
#> Working on grid time point: 2010-09-17
#> Working on grid time point: 2010-09-18
#> Working on grid time point: 2010-09-19
#> Working on grid time point: 2010-09-20
#> Working on grid time point: 2010-09-21
#> Working on grid time point: 2010-09-22
#> Working on grid time point: 2010-09-23
#> Working on grid time point: 2010-09-24
#> Working on grid time point: 2010-09-25
#> Working on grid time point: 2010-09-26
#> Working on grid time point: 2010-09-27
#> Working on grid time point: 2010-09-28
#> Working on grid time point: 2010-09-29
#> Working on grid time point: 2010-09-30
#> Working on grid time point: 2010-10-01
#> Warning in makeNetCDF_file(ncdfFilename = "AWAP_demo.nc", updateFrom = as.Date("2010-08-01", : Thesolar radiation data netCDF file will not be built or updated.
#> Data construction FINISHED..

Set the points for data extraction

Set coordinates to the location of one groundwater bore and one rainfall gauge and then convert the points to a spatial object and set projection to GDA94.

coordinates.data = data.frame( ID =c('Bore-10084446','Rain-63005'),
                               Longitude = c(153.551875, 149.5559),
                               Latitude =  c(-28.517974,-33.4289))

sp::coordinates(coordinates.data) <- ~Longitude + Latitude

sp::proj4string(coordinates.data) = '+proj=longlat +ellps=GRS80 +no_defs'
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition

Extract daily precipitation data

Extract the daily precipitation at the two locations. The data is extracted from the netCDF file ncdfFilename and between the dates extractFrom and extractTo. The other AWAPer variables are not extracted because getTmin, getTmax, getVprp, getSolarrad and getET are set to F. Note the netCDF file ncdfFilename must be in the working directory or the full file path must be given.

climateData.data = extractCatchmentData(ncdfFilename='AWAP_demo.nc',
                                        extractFrom=as.Date("2010-08-01","%Y-%m-%d"),
                                        extractTo=as.Date("2010-10-01","%Y-%m-%d"),
                                        catchments=coordinates.data,
                                        getTmin=F, getTmax=F, getVprp=F, getSolarrad=F, getET=F)
#> Warning in sp::proj4string(catchments): CRS object has comment, which is lost in output

#> Warning in sp::proj4string(catchments): CRS object has comment, which is lost in output
#> WARNING: The projection string of the catchment boundaries does not appear to be +proj=longlat +ellps=GRS80. Attempting to transform coordinates...
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
#> Extraction data summary:
#>     NetCDF non-solar radiation climate data exists from  1900-01-01  to  2010-10-01
#>     Data will be extracted from  2010-08-01  to  2010-10-01  at  2  catchments
#> Starting data extraction:
#> ... Building catchment weights:
#> ... Starting to extract data across all catchments:
#> ... Calculating catchment weighted daily data.
#> Data extraction FINISHED..

Plot the daily precipitation at each site

Plot time series of the etracted daily precipitation.

for (i in 1:nrow(coordinates.data)){
  filt = climateData.data$CatchmentID.ID == coordinates.data$ID[i]

  data2plot = climateData.data$precip_mm[filt]
  dates = ISOdate(climateData.data$year[filt],climateData.data$month[filt],climateData.data$day[filt])

  if (i==1){
    plot(dates, data2plot, main='Extracted precip.', ylab='Precip [mm/d]', xlab='Date [day/month]', type='l',lty = i+1, cex=0.2)
  } else {
    lines(dates, data2plot, lty = i+1)
  }
}
legend('topright',legend=c('Bore-10084446','Rain-63005'),lty=2:3)

Check AWAPer data against precipitation gauge date

Use the following observed precipitation (from gauge_63005) to check the AWAPer data. Note, the observed data is hardcoded in and sourced from the Australian Bureau of Meteorology.

obsPrecip <- data.frame(
  year= c(2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010,
          2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010,
          2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010),
  month = c(8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9,
            9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10),
  day = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 1, 2, 3, 4, 5,
          6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 1),
  precip_mm = c(0.6, 5.2, 0.8, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 15.8, 15.6, 7.6, 0.7, 0.4, 1.4, 1.0, 0.0, 0.0, 30.4, 1.0, 0.0,
                0.0, 5.0, 2.2, 0.3, 0.8, 13.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 13.8, 15.2, 0.4, 0.2, 0.0, 0.0, 12.4, 0.9,
                0.0, 0.0, 0.1, 13.0, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3, 0.0, 0.0, 0.0, 0.0))

# Plot the observed vs AWAPer rainfall at gauge ID 63005
filt2 = climateData.data$CatchmentID.ID=='Rain-63005'
plot(obsPrecip$precip_mm,climateData.data$precip_mm[filt2],
     xlim = c(0,35),ylim = c(0,35),
     main='Obs. vs. AWAPer precip.',
     xlab='Obs. [mm/d]', ylab='AWAPer [mm/d]', cex=0.2)
abline(0,1, col='grey', lty=2)

# Plot the cumulative observed vs AWAPer rainfall at gauge ID 63005
plot(cumsum(obsPrecip$precip_mm),cumsum(climateData.data$precip_mm[filt2]),
     xlim = c(0,175),ylim = c(0,175),
     main='Cum. obs. vs. AWAPer precip.',
     xlab='Obs. [mm]', ylab='AWAPer [mm]', type='l', cex=0.2)
abline(0,1, col='grey', lty=2)

AWAPer documentation built on Feb. 18, 2021, 1:06 a.m.