# load the library, marvel at the many dependencies :D
library(SpatialVx)
# plot options for this notebook
options( repr.plot.width=10, repr.plot.height=6, repr.plot.res=150 )
Lade nötiges Paket: spatstat
Lade nötiges Paket: spatstat.data
Lade nötiges Paket: spatstat.geom
spatstat.geom 3.1-0
Lade nötiges Paket: spatstat.random
spatstat.random 3.1-4
Lade nötiges Paket: spatstat.explore
Lade nötiges Paket: nlme
spatstat.explore 3.1-0
Lade nötiges Paket: spatstat.model
Lade nötiges Paket: rpart
spatstat.model 3.2-1
Lade nötiges Paket: spatstat.linnet
spatstat.linnet 3.0-6
spatstat 3.0-3
For an introduction to spatstat, type ‘beginner’
Lade nötiges Paket: fields
Lade nötiges Paket: spam
Spam version 2.9-1 (2022-08-07) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attache Paket: ‘spam’
Die folgenden Objekte sind maskiert von ‘package:base’:
backsolve, forwardsolve
Lade nötiges Paket: viridis
Lade nötiges Paket: viridisLite
Try help(fields) to get started.
Lade nötiges Paket: smoothie
Lade nötiges Paket: smatr
Lade nötiges Paket: turboEM
Lade nötiges Paket: doParallel
Lade nötiges Paket: foreach
Lade nötiges Paket: iterators
Lade nötiges Paket: parallel
Lade nötiges Paket: numDeriv
Lade nötiges Paket: quantreg
Lade nötiges Paket: SparseM
Attache Paket: ‘SparseM’
Das folgende Objekt ist maskiert ‘package:spam’:
backsolve
Das folgende Objekt ist maskiert ‘package:base’:
backsolve
Attache Paket: ‘turboEM’
Die folgenden Objekte sind maskiert von ‘package:numDeriv’:
grad, hessian
R 101: use help() or ? to access the documentation. In this Notebook, you can also click on a bit of code and hit shift + tab.
As an example, try ?help to learn about the help function
help(package=SpatialVx)
Dokumentation für Paket ‘SpatialVx’
Information für Paket ‘SpatialVx’
Description:
Package: SpatialVx
Version: 1.0-1
Date: 2022-10-25
Title: Spatial Forecast Verification
Authors@R: c(person("Eric", "Gilleland", role = c("aut", "cre"
), email = "EricG@ucar.edu" ), person("Kim",
"Elmore", role = "ctb", email =
"kim.elmore@noaa.gov" ), person("Caren", "Marzban",
role = "ctb", email = "marzban@stat.washington.edu"
), person("Matt", "Pocernich", role = "ctb", email
= "matt_pocernich@hotmail.com" ), person("Gregor",
"Skok", role = "ctb", email =
"Gregor.Skok@fmf.uni-lj.si" ) )
Author: Eric Gilleland [aut, cre], Kim Elmore [ctb], Caren
Marzban [ctb], Matt Pocernich [ctb], Gregor Skok
[ctb]
Maintainer: Eric Gilleland <EricG@ucar.edu>
Depends: R (>= 4.1.0), spatstat (>= 2.0-0), fields (>= 6.8),
smoothie, smatr, turboEM
Imports: spatstat.geom, spatstat.linnet, spatstat.model,
distillery, maps, methods, boot, CircStats,
fastcluster, waveslim
Description: Spatial forecast verification refers to verifying
weather forecasts when the verification set
(forecast and observations) is on a spatial field,
usually a high-resolution gridded spatial field.
Most of the functions here require the forecast and
observed fields to be gridded and on the same grid.
For a thorough review of most of the methods in
this package, please see Gilleland et al. (2009)
<doi: 10.1175/2009WAF2222269.1> and for a tutorial
on some of the main functions available here, see
Gilleland (2022) <doi: 10.5065/4px3-5a05>.
License: GPL (>= 2)
URL: https://projects.ral.ucar.edu/icp/SpatialVx/
BugReports: https://staff.ral.ucar.edu/ericg/
NeedsCompilation: no
Packaged: 2022-11-08 14:06:47 UTC; gille
Repository: CRAN
Date/Publication: 2022-11-08 15:50:03 UTC
Built: R 4.2.2; ; 2022-11-11 07:40:10 UTC; unix
Index:
Aindex Area Index
CSIsamples Forecast Verification with Cluster Analysis:
The Variation
Cindex Connectivity Index
CircleHistogram Circular Histogram
EBS Elmore, Baldwin and Schultz Method for Field
Significance for Spatial Bias Errors
ExampleSpatialVxSet Simulated Spatial Verification Set
FQI Forecast Quality Index
FeatureAxis Major and Minor Axes of a Feature
FeatureFinder Threshold-based Feature Finder
FeatureMatchAnalyzer Analyze Features of a Verification Set
FeatureProps Single Feature Properties
FeatureTable Feature-Based Contingency Table
Fint2d 2-d Interpolation
GFSNAMfcstEx Example Verification Set
Gbeta Spatial-Alignment Summary Measures
GeoBoxPlot Geographic Box Plot
LocSig Temporal Block Bootstrap Keeping Locations in
Space Constant
MCdof Monte Carlo Degrees of Freedom
MergeForce Force Merges in Matched Feature Objects
Mij Raw Image Moments.
OF Optical Flow Verification
S1 S1 Score, Anomaly Correlation
Sindex Shape Index
SpatialVx-package Spatial Forecast Verification
TheBigG The Spatial Alignment Summary Measure Called G
UKobs6 Example Precipitation Rate Verification Set
(NIMROD)
abserrloss Loss functions for the spatial prediction
comparison test (SPCT)
bearing Bearing from One Spatial Location to Another
binarizer Create Binary Fields
calculate_FSSvector_from_binary_fields
Calculate FSS Values
calculate_FSSwind Calculate FSSwind Score Value
calculate_dFSS Calculate dFSS Score Value
censqdelta Centered on Square Domain Baddeley's Delta
Metric
centdist Centroid Distance Between Two Identified
Objects
clusterer Cluster Analysis Verification
combiner Combine Features/Matched Objects
compositer Create Composite Features
deltamm Merge and/or Match Identified Features Within
Two Fields
disjointer Identify Disjoint Sets of Connected Components
expvg Exponential Variogram
expvgram Exponential Variogram
fss2dPlot Create Several Graphics for List Objects
Returned from hoods2d
fss2dfun Various Verification Statistics on Possibly
Neighborhood-Smoothed Fields.
gmm2d 2-d Gaussian Mixture Models Verification
griddedVgram Variograms for a Gridded Verification Set
hiw Spatial Forecast Verification Shape Analysis
hoods2d Neighborhood Verification Statistics for a
Gridded Verification Set.
hoods2dPlot Quilt Plot and a Matrix Plot.
hump Simulated Forecast and Verification Fields
imomenter Image Moments
interester Feature Interest
iwarper Image Warping By Hand
locmeasures2d Binary Image Measures
locperf Localization Performance Measures
lossdiff Test for Equal Predictive Ability on Average
Over a Regularly Gridded Space
make.SpatialVx Spatial Verification Sets - SpatialVx Object
metrV Binary Location Metric Proposed in Zhu et al.
(2011)
minboundmatch Minimum Boundary Separation Feature Matching
obs0601 Spatial Forecast Verification Methods
Inter-Comparison Project (ICP) Test Cases and
other example verification sets
optflow Optical Flow
pphindcast2d Practically Perfect Hindcast Neighborhood
Verification Method
rigider Rigid Transformation
saller Feature-based Analysis of a Field (Image)
spatbiasFS Field Significance Method of Elmore et al.
(2006)
spct Spatial Prediction Comparison Test
structurogram Structure Function for Non-Gridded Spatial
Fields.
structurogram.matrix Structure Function for Gridded Fields
surrogater2d Create Surrogate Fields
thresholder Apply a Threshold to a Field
upscale2d Upscaling Neighborhood Verification on a 2-d
Verification Set
variographier Variography Score
vxstats Some Common Traditional Forecast Verification
Statistics.
warper Image Warp
waveIS Intensity Scale (IS) Verification
wavePurifyVx Apply Traditional Forecast Verification After
Wavelet Denoising
waverify2d High-Resolution Gridded Forecast Verification
Using Discrete Wavelet Decomposition
You can find the full documentation here http://cran.r-project.org/web/packages/SpatialVx/SpatialVx.pdf. There is also an extended Tutorial by Eric Gilleland, which can be found here http://dx.doi.org/10.5065/4px3-5a05 .
R objects can be saved with save() and loaded with load(). Let's load some forecasts and observations!
load("MVdat.rdata", verbose=TRUE)
Loading objects: MVdat
# dimensions of the data: nx x ny x n_datasets x n_timesteps
print( dim( MVdat ) )
print( dimnames( MVdat ) )
[1] 209 193 7 8 [[1]] NULL [[2]] NULL [[3]] [1] "VERA" "BOLAM007+1d" "MOL00225+1d" "BOLAM007+2d" "MOL00225+2d" [6] "BOLAM007+3d" "MOL00225+3d" [[4]] [1] "2007-06-22 03:00" "2007-06-22 19:00" "2007-07-20 11:00" "2007-09-27 07:00" [5] "2007-09-27 21:00" "2007-08-08 01:00" "2007-08-08 21:00" "2007-07-10 13:00"
# basic plot methods in R: use par() to set parameters, use functions like plot(), image.plot(), ... to make plots
par( mfrow=c(2,2), mar=c(2,2,3,4) )
image.plot( MVdat[ , , 1, 3 ], main="obs" )
image.plot( MVdat[ , , 2, 3 ], main="forecast 1" )
image.plot( MVdat[ , , 3, 3 ], main="forecast 2" )
image.plot( MVdat[ , , 4, 3 ], main="forecast 3" )
To try out the spatial verification methods, first put the data into a SpatialVx object:
# ?make.SpatialVx
obj <- make.SpatialVx( X = MVdat[ , , 1, ], # observation
Xhat = list( MVdat[ , , 2, ], MVdat[ , , 3, ], MVdat[ , , 4, ], MVdat[ , , 5, ], MVdat[ , , 6, ], MVdat[ , , 7, ] ), # competing forecasts
thresholds = c(0.1, 0.5,1,5,10), # precipitation thresholds of interest
obs.name = "VERA",
model.name = dimnames( MVdat )[[3]][2:7]
)
# plot method
plot(obj, model=2, time.point=3)
Try out the neighbourhood verification, for example FSS. The thresholds were already set when we created the object, now just select the neighbourhood sizes (smoothing levels) and check out the results.
fss <- hoods2d(obj, which.methods = c("fss"), levels=c(1,5,11,15,25,35), model=2, time.point=3, verbose=TRUE )
Looping through thresholds. Setting up binary objects for threshold 1 Looping through levels. Neighborhood length = 1 Neighborhood length = 5 Neighborhood length = 11 Neighborhood length = 15 Neighborhood length = 25 Neighborhood length = 35 Setting up binary objects for threshold 2 Looping through levels. Neighborhood length = 1 Neighborhood length = 5 Neighborhood length = 11 Neighborhood length = 15 Neighborhood length = 25 Neighborhood length = 35 Setting up binary objects for threshold 3 Looping through levels. Neighborhood length = 1 Neighborhood length = 5 Neighborhood length = 11 Neighborhood length = 15 Neighborhood length = 25 Neighborhood length = 35 Setting up binary objects for threshold 4 Looping through levels. Neighborhood length = 1 Neighborhood length = 5 Neighborhood length = 11 Neighborhood length = 15 Neighborhood length = 25 Neighborhood length = 35 Setting up binary objects for threshold 5 Looping through levels. Neighborhood length = 1 Neighborhood length = 5 Neighborhood length = 11 Neighborhood length = 15 Neighborhood length = 25 Neighborhood length = 35 Time difference of 1.294549 secs
plot(fss)
Dashed lines show the scores for a uniform forecast, dotted lines correspond to a random forecast.
Try out other models, or time steps, or neighbourhood sizes if you want.
https://doi.org/10.1175/2008MWR2415.1 , https://doi.org/10.1175/2009WAF2222271.1
First, let us try to find some features:
feat <- FeatureFinder( obj, model=2, time.point=3,
thresh=1, min.size=5, smoothpar=1 )
summary(feat)
plot(feat)
Verification field ( ) feature properties:
centroidX centroidY area OrientationAngle AspectRatio Intensity0.25
[1,] 94.37893 118.5739 636 141.4859 0.3797866 1.5575
[2,] 24.97767 145.9495 1030 128.8125 0.6870459 1.3400
Intensity0.9
[1,] 6.740
[2,] 4.211
Forecast field ( ) feature properties:
centroidX centroidY area OrientationAngle AspectRatio Intensity0.25
[1,] 6.50000 10.00000 6 NA NA 1.1996462
[2,] 18.16667 23.83333 6 135.00000 0.8000000 0.7239479
[3,] 87.60784 86.19608 51 29.33890 0.3338269 1.2836358
[4,] 86.73913 106.69565 23 49.41227 0.5000290 0.8597148
[5,] 73.63158 107.94737 19 28.56521 0.8179922 1.2710418
[6,] 23.33333 117.77778 9 115.85537 0.4475109 1.4393665
[7,] 15.57230 140.90411 657 146.08604 0.7176647 1.2918923
[8,] 60.50000 123.00000 6 NA NA 1.3964723
[9,] 60.53005 143.97541 366 111.46316 0.4054323 1.3629926
[10,] 31.40000 152.70000 10 136.51346 0.4078405 1.1410787
Intensity0.9
[1,] 1.912341
[2,] 4.308489
[3,] 3.855635
[4,] 3.273141
[5,] 2.512367
[6,] 2.477558
[7,] 5.694697
[8,] 2.272766
[9,] 8.978880
[10,] 1.366515
Features are detected by
min.size grid points You can play with the threshold, min.size, and the smoothing strength (kernel size) smoothpar. Next, let us compute the SAL scores.
SAL <- saller(feat)
summary(SAL)
[1] "" [1] "" "VERA" "MOL00225+1d" Structure Component (S): -0.9818771 Amplitude Component (A): -0.3585697 Location Component (L): 0.1956943
S and A range from $-2$ to $2$, the location score is in $[0,2]$. The scores indicate that the forecast is far too small-scaled or peaked ($S<<0$) and somewhat underestimates the amplitude. The location is apparently quite good, but this score is often not very informative for large domains with several discrete objects.
Try verifying a different forecast for this day, or look at another day.
IS1 <- waveIS( obj, model=2, time.point=3, verbose=TRUE )
Looping through thresholds = Threshold 1 Threshold 2 Threshold 3 Threshold 4 Threshold 5 Finished computing DWTs for each threshold.
plot(IS1)
Careful though, for this score it makes a big difference whether your data was $2^n \times 2^n$ or not! Let's try it out.
# a simple function to pad our fields with zeros to a square of size N x N
pad <- function( x, N){
nx <- nrow(x)
ny <- ncol(x)
px <- 1:nx + ceiling( (N - nx)/2 )
py <- 1:ny + ceiling( (N - ny)/2 )
res <- array( dim=c(N,N), data=0 )
res[ px, py ] <- x
return( res )
}
dyadic_dat <- make.SpatialVx( pad( MVdat[ ,,1,3 ], 256 ), pad( MVdat[ ,,3,3 ], 256 ), thresholds = c(0.1, 0.5,1,5,10) )
plot(dyadic_dat)
IS2 <- waveIS( dyadic_dat, model=1, time.point=1 )
plot(IS2)
Find a vector field that transports the predicted precipitation to the observed locations. Here, we use my implementation of the Pyramid matching algorithm used by the Displacement and Amplitude score (https://doi.org/10.1175/2009WAF2222247.1) .
# not implemented in SpatialVx, still interesting though
source( "KCflow.r" )
flow <- KCflow( obs=MVdat[ ,,1,3 ], forc=MVdat[ ,,3,3 ],
F = 5, # maximum displacement scale
f = 1, # minimum displacement scale
d = 1, # maximum displacement to test at each scale
sigma = 2, # smoothing applied to the fields
smflow =3 # smoothing applied to the flow vectors
)
plot(flow)
Apart from some artifacts, the algorithm has smoothly shifted some of the misplaced precipitation to the observed locations. Let us look at the output of the function:
str(flow)
List of 5 $ obs : num [1:209, 1:193] 0.46 0.45 0.45 0.44 0.43 ... $ forc : num [1:209, 1:193] 0 0 0 0 0 0 0 0 0 0 ... $ dx : num [1:209, 1:193] 0.758 1.083 1.432 1.798 2.171 ... $ dy : num [1:209, 1:193] -7.91 -7.94 -7.96 -7.98 -7.99 ... $ warped: num [1:209, 1:193] 0 0 0 0 0 ... - attr(*, "class")= chr "oflow"
How much did each pixel move?
disp <- sqrt( flow$dx**2 + flow$dy**2 )
disp[MVdat[ ,,3,3 ] <0.1] <- NA
image.plot( disp )
# plot histogram of displacement lengths
hist(disp)
We could also look at the reverse problem of deforming the observation to the predicted locations:
inverse_flow <- KCflow( obs=MVdat[ ,,3,3 ], forc=MVdat[ ,,1,3 ] )
plot(inverse_flow)
disp2 <- sqrt( inverse_flow$dx**2 + inverse_flow$dy**2 )
disp2[MVdat[ ,,1,3 ] <0.1] <- NA
image.plot( disp2 )
hist( disp2 )
For example https://doi.org/10.1175/WAF-D-10-05061.1
# compute distance maps with functions from the spatstat library
dm_obs = distmap( solutionset( im( MVdat[ ,,1,3 ] )>1 ) )$v
dm_for = distmap( solutionset( im( MVdat[ ,,3,3 ] )>1 ) )$v
par(mfrow=c(1,2))
image.plot(dm_obs)
contour( MVdat[ ,,1,3 ], add=TRUE, col="white" , levels=c(1,10))
image.plot(dm_for)
contour( MVdat[ ,,3,3 ], add=TRUE, col="white" , levels=c(1,10))
Scores can be computed based on the difference in these distance maps:
image.plot(dm_obs-dm_for)
# distance measure scores:
dist_scores <- locmeasures2d( obj, model=2, time.point=3,
k=c(4, 0.975), alpha=c(0.1,0.9)
)
summary(dist_scores)
[1] ""
Comparison for:
[1] "" "VERA" "MOL00225+1d"
Threshold(s) is (are):
[1] "0.1" "0.5" "1" "5" "10"
Baddeley Delta Metric with c = Inf
0.1 0.5 1 5 10
p = 2; NA NA NA NA NA
Hausdorff distance
0.1 0.5 1 5 10
[1,] 90.14214 101.0122 87.91169 121.6518 85.24264
Quantile (if k in (0,1) or k-th highest (if k = 1, 2, ...) difference in distance maps.
0.1 0.5 1 5 10
k = 4; 90.14214 99.84062 87.91169 121.6518 85.24264
k = 0.975; 61.83574 85.56854 50.40833 104.6646 74.35534
Mean error distance (Miss)
0.1 0.5 1 5 10
[1,] 6.924713 6.379885 4.3018 11.72145 31.66671
Mean error distance (FalseAlarm)
0.1 0.5 1 5 10
[1,] 3.181524 4.707234 6.342921 16.02486 53.21314
Mean square error distance (Miss)
0.1 0.5 1 5 10
[1,] 226.5624 238.7637 50.74211 194.9927 1004.406
Mean square error distance (FalseAlarm)
0.1 0.5 1 5 10
[1,] 53.19609 153.4358 140.2165 402.3831 3288.283
Partial Hausdorff distance
0.1 0.5 1 5 10
k = 4; NA NA NA NA NA
k = 0.975; NA NA NA NA NA
Pratt's figure of merit (FOM)
0.1 0.5 1 5 10
alpha = 0.1; 0.5792179 0.6197010 0.6099739 0.18868909 0.0062005025
alpha = 0.9; 0.4260656 0.4572163 0.4093034 0.09028741 0.0006951138
Careful, results could be sensitive to small features in otherwise empty regions.