I am having trouble calculating the pixel by pixel correlation coefficient between two datasets. My current code takes in two folders full of rasters and creates two independent raster stacks. These rasters all have the same cellsize and extent. I then try and take the correlation coefficient (spearman in my case) between the column values of those rasters.  library(raster)  r &lt;- raster() raster1 &lt;- list.files(path = "data/List1", pattern = "*.tif$", full.names = T) raster2 &lt;- list.files(path = "data/List2", pattern = "*.tif$", full.names = T) l1 &lt;- stack(raster1) l2 &lt;- stack(raster2)  list1Values &lt;- values(l1) list2Values &lt;- values(l2) corValues &lt;- vector(mode = 'numeric')  for (i in 1:dim(list1Values)[1]){   corValues[i] &lt;- cor(x = list1Values[i,], y = list2Values[i,], method = 'spearman') }  corRaster &lt;- setValues(r, values = corValues)   However at the correlation for loop it gives six error messages saying  In cor(x = list1Values[i, ], y = list2Values[i, ], method = "spearman") :   the standard deviation is zero   Ignoring that and continuing on, the last line errors out saying  Error in setValues(r, values = corValues) :    length(values) is not equal to ncell(x), or to 1   My initial guess was that since the datasets have a lot of NA values (In this case it is a lot of open ocean that there is no data for), that could cause problems, so I added the   use = "complete.obs"   parameter to the cor function. This errors saying  Error in cor(x = list1Values[i, ], y = list2Values[i, ], method = "spearman",  :    no complete element pairs   My guess is that this last error is telling me the matrices it created don't line up, and no non-NA cells match. I have no clue how this is possible because I have been working with these rasters for years and they certainly line up. Other than that, I don't know why this isn't working.    This question is not a duplicate of  in any way. First I am not using point data, but instead raster. I am also interested in a simple pearson and/or spearman correlation coefficient NOT cross correlation nor spatial correlation. I'm also not using two layers, but two sets of hundreds of layers (making the statistical analysis of the correlation valid). Lastly, this code is doing fundamentally different things than that post and I am asking for specific help with the errors arising.
120501
I have a point shapefile concerning asthma instances, a point shapefile concerning rainfall, and a polygon shapefile with population information i.e gender, ages etc.  I would like to research if there is a spatial relationship between the location of asthma and rainfall, and also if there is a spatial relationship between instances of asthma and  the population data. Ultimately I would to be able to say:     There are/are not higher instances of asthma in areas that receive more rainfall and/are there are/are not higher instances of asthma in areas with more men/women.   I have done research on the areas of spatial correlation and spatial autocorrelation which I assumed would be related to this area. However, it seems these only correspond to data within one map layer i.e. you cannot technically look for spatial correlation between two map layers.   Is there any way in R to assess correlation/a relationship between two point layers, or a point layer and a polygon layer?
96348
I have a single band raster image of 30m by 30m resolution. I have already created a polygon fishnet (grid) of 300m by 300m with point labels over that. How can I extract the mean pixel value of the raster image overlaying each feature in the fishnet grid? Is it possible with Arcgis?  Please inform me.
64444
I am trying to use the tool  on a File Geodatabase raster, but I am getting some weird errors.  I've used this tool multiple times before with no issue, but not since I upgraded to 10.1, so I'm wondering if there is an issue with it.  The errors I got back so far have been:   This error happened when I specified a shapefile as the output: The specified geodataset is invalid for the conversion. The syntax help for the tool you are using should indicate the supported types of input data. For example, for the particular conversion, a rasterlayer may be an invalid type of geodataset.  Next I have an error about it needed at least one band in the raster, the raster does has one band. This occurred when I specified a File Geodatabase as the output. Finally I had the good ole Error: 999999 for no apparent reason.   The raster is a DEM derived from LiDAR that was not created by a different company. It's too large to be made available for testing, but here are its specs:  Raster_Information  ColumnsandRows  5119,8241 NumberofBands   1 CellSizeXY      1, 1 UncompressedSize    160.93 MB Format      FGDBR Source_Type     Generic Pixel_Type  floating point Pixel_Depth     32 Bit Colormap    absent
120408
I am running a factor analysis in SPSS and get a "matrix is not positive definite" error from my correlation matrix. I've tried removing correlated variables, but I have to remove all variables down to those with correlations of r = 0.8 before the "positive definite" issue is resolved. This seems weird, as I have read that this problem normally arises where two or more vbls are nearly perfectly collinear (to me, r = 0.8 is a high correlation, but not near-perfect.) Can anyone suggest a way of dealing with this issue that doesn't just involve arbitrary removal of variables until the problem goes away? Or perhaps explain why I am getting the issue even when I seem to be removing all the really highly correlated variables?              [edit by Gottfried Helms:] By comment to an answer the OP says, there are 14 variables on 10 cases per correlation.
120032
I'm trying to incorporate different practices to use when exploring a new data set. Especially, how to examine the association between two variables.  Steps for example (not necessarily by order):   plot a y-by-x scatter plot of the raw data to see the relationship visually. compute summary statistics for each variable (mean and sd) compute correlation coefficient r draw the OLS regression line, compute its slope and intercept etc....   I've come across the "SD line" in Freedman's Statistics book, which is defined as:     "the line that goes through the point of averages and climbs at the rate of one vertical SD for each horizontal SD" Freedman, D., Pisani, R., &amp; Purves, R. (2007). Statistics (4th edn).   Since this book ("Statistics") is a canonical textbook, I consider its choice to discuss the SD line as an indication for the line's importance. However, a simple google search for the term "SD line" doesn't yield as many independent results. Most of them come directly from Freedman's book. This tells me it's not a central concept in bivariate analyses in general.  When comparing the SD with the OLS regression line, it seems like the regression line is more informative (than the SD line) for predicting y from x. Therefore, I'm wondering if bothering to plot the SD line has any benefit or added value that I would not already have when plotting the regression line.  Example using mtcars dataset, focusing on association between weight and mpg  data(mtcars)  ## calculate means mean_wt &lt;- mean(mtcars$wt) mean_mpg &lt;- mean(mtcars$mpg)  ## calculate standard deviations sd_wt &lt;- sd(mtcars$wt) sd_mpg &lt;- sd(mtcars$mpg)  ## scatter plot plot(x = mtcars$wt, y = mtcars$mpg)  ## add the "point of averages" points(mean_wt, mean_mpg, col = "red", cex = 1.5, pch = 16)  ## calculate the slope of the sd line slope &lt;- -1*sd_mpg/sd_wt  ## plot the sd line curve(expr = x*slope + (mean_mpg - slope*mean_wt), add = TRUE, col = 'blue', lwd = 2, type = "l", lty = 2)  ## plot the regression line model &lt;- lm(mpg ~ wt, data = mtcars) abline(model, col = "orange", lwd = 2)  ## legend legend("topright",        legend = c("Regression line", "SD line"),        col = c("orange", "blue"),        lty = c(1, 2),        lwd = c(2, 2))     Thus, my question: how can the SD line increase one's understanding about the relationship between two variables, in a way that's either adding or complementing on what the regression line already tells?
82579
I am trying to do kriging by gstat package through R.   Here is my Code:  krg&lt;-krige(var ~ 1,data,grid,model=vgm,nmax=30)   But each time I recive this warning:  1: In predict.gstat(g, newdata = newdata, block = block,  ... :    Warning: Covariance matrix (nearly) singular at location [1,4501,0]: skipping...   and when I try to plot it by spplot this error appears:  Error in seq.default(zrng[1], zrng[2], length.out = cuts + 2) :    'from' cannot be NA, NaN or infinite In addition: Warning messages: 1: In asMethod(object) :   complete map seems to be NA's -- no selection was made 2: In min(x) : no non-missing arguments to min; returning Inf 3: In max(x) : no non-missing arguments to max; returning -Inf   Here is the full code: You can download my data through this link:  and here is my full code:  library(gstat)                                          library(sp) library(raster)  coordinates(Mydata) &lt;- ~East + North proj4string(Mydata) &lt;- CRS("+proj=utm +zone=10+datum=WGS84") ext &lt;- as(extent(Mydata), "SpatialPolygons") r &lt;- rasterToPoints(raster(ext), spatial = TRUE) proj4string(r) &lt;- proj4string(Mydata)  v1 &lt;- variogram(cpue_female ~ 1, Mydata) vgm_mod_omni &lt;- vgm(psill = 0.0003, model = "Mat", range = 50000, nugget = 33, kappa = 0.5) G1 &lt;- krige(cpue_female ~ 1, Mydata, r, vgm_mod_omni, nmax = 30) gridded(G1) &lt;- TRUE spplot(G1)
51826
I am trying to loop through multiple rasters and extract cells that are above 0.1 The code I am trying to use is:  import arcpy  arcpy.CheckOutExtension('Spatial') arcpy.env.workspace=r'F:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master' #pathway to all rasters in workspace directory rasters = arcpy.ListRasters('*.tif*') for raster in rasters:      #give new file a name      outraster = raster.replace('.tif','_water.tif')      arcpy.gp.RasterCalculator_sa("""raster &gt;= 0.1""",outraster) print('Done Processing')   For some reason I can't copy and paste the error (python actually shuts down when I run this code) but here is a screen shot of part of it.    EDIT: My parameters have changed and I now need everything >=-0.4 and but !=0.00.  I am trying to use:  import arcpy  from arcpy.sa import * arcpy.CheckOutExtension('Spatial') arcpy.env.workspace=r'D:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master' #pathway to all rasters in workspace directory rasters = arcpy.ListRasters('*.tif*') for raster in rasters:     #Save temp raster to disk with new name     ras = Raster(raster)     outraster = Con(ras &gt;= -0.4 &amp; ras !=0.00, ras)     outraster.save(raster.replace('.tif','_water.tif')) print('Done Processing')   but this returns:  ValueError: The truth value of a raster is ambiguous. Invalid use of raster with Boolean operator or function. Check the use of parentheses where applicable.   EDIT:  I think I just needed to change this line:   outraster = Con((ras &gt;= -0.4) &amp; (ras !=0.00), ras)
47695
Similarly to this post    the rasterize function in the raster package is giving me problems with enclosed polygons. In my case, some enclosed polygons are rasterized with cell values equal to the enclosing polygons, as they were somehow 'merged'.   Here's my sample code  require(sp) require(rgdal) require(rgeos) require(raster)  ## set wd wd &lt;- 'M:/foo' setwd(wd)  ## open shapefile sample.poly &lt;- readOGR(dsn=wd, layer='sample')  ## rasterize resol &lt;- 10 r &lt;- raster(res=resol, ext=extent(sample.poly)) sample.r &lt;- rasterize(sample.poly, r, 'Code')   ## export files writeRaster(x=sample.r, filename='sample.tif', format='GTiff') writeOGR(obj=sample.poly, dsn=wd, layer='sample', driver='ESRI Shapefile')   Here's the structure of table of attributes  str(sample.poly@data) 'data.frame':   11 obs. of  2 variables: $ OBJECTID: num  6249 14593 14614 15434 15683 ... $ Code    : int  33200 32110 33300 32410 32210 31210 33300 32110 32210 31210 ..   data$Code are CORINE LCM-like codes. I already tried to set them as characters instead of integers, but the same error occurred.  Here's how the sample shapefile looks like (in QGIS)    and here's the result of rasterize    As you can see, the enclosed polygon marked with yellow arrow is properly rasterized, while the one marked with the red arrow is 'dissolved' in the enclosing polygon.   Any suggestion? Is it a bug? This is just a sample taken from a larger shapefile, with dozens of such errors.
18045
Let's say I test how variable Y depends on variable X under different experimental conditions and obtain the following graph:    The dash lines in the graph above represent linear regression for each data series (experimental setup) and the numbers in the legend denote the Pearson correlation of each data series.  I would like to calculate the "average correlation" (or "mean correlation") between X and Y. May I simply average the r values? What about the "average determination criterion", $R^2$? Should I calculate the average r and than take the square of that value or should I compute the average of individual $R^2$'s?
33403
I want to do something relatively simply but it seems to be very difficult to do, though I am very new to using R for spatial data so apologies if I've missed something obvious.  I have a list containing 161 shapefiles of bird distributions in the world, they are each in the format of a SpatialPolygonsDataFrame. I also have world climate data in a raster stack format for 8 variables. Both sets of files are in the same projection (lat and long WGS84). I want to create new columns for the climate data with 1's where the species polygon overlaps with the climate data and 0 where it doesn't. I tried various functions to do this, quite a few initial options didn't work, the best option appeared to be:  pnt.in.poly(XY coordinates of climate data points, XY coordinates of species polygon)   I changed my raster of climate data into points to get the XY coordinates. I unified the polygons for each species shapefile and then used:  species@polygons[[1]]@Polygons[[1]]@cords   to get the coordinates for my species polygons as a matrix. However, I discovered that:  species@polygons[[1]]@Polygons[[1]]@coords    wasn't actually getting all the coordinates of the polygons only some of them.  I tried using fortify from the ggplot package, and it gives me a dataframe with coordinates but some of them are for holes in the polygons. I ran that with the pnt.in.poly function excluding the holes but it took an hour to run and when I looked at the output maps after it matched for some of the range but not all of it- there were holes in the distribution where there shouldn't be holes.   I also tried this method that I found somewhere else:  lonlat &lt;- lapply(slot(species, "polygons"), function(x) lapply(slot(x,"Polygons"), function(y) slot(y, "coords")))    which gives me a list with 436 items- but I'm still not sure how to get the coordinates into a matrix form that will work in pnt.in.poly.  I'm sure this should be easier to do! Does anyone know a better way?
71444
Since a few days I do not get ahead when trying to compare two Pearson correlation coefficients. Imagine that I've got two datasets where on each I do a correlation between Land Surface Temperature and an urban metric. The datasets are different in their length, so the first one has round about 160.000 observables and the second one has about 2400 observables. For the correlation on the first dataset I get a Pearson of -0.74 and for the second I get -0.885. Now I want to find out whether these coefficients are significantly different from each other. Is there any appropriate method you could suggest?  I already played around with the Fisher-Z-Transformation, but from my point of view with no purposeful results. When I calculate the Fisher-Z on my coefficients, it results in 20.95 (Fisher-Z).  I would be very happy about suitable advices.
60084
I am trying to figure out how to convert a correlation matrix (R) to a covariance matrix (S) for input into a random number generator that only accepts S (rmvnorm("mvtnorm") in R)   library("mvtnorm")   TRUTH= 0.8 # target correlation value between X1 and X2 R &lt;- as.matrix(data.frame(c(1, TRUTH), c(TRUTH, 1))) V &lt;- diag(c(sqrt(1), sqrt(1))) # diagonal matrix of sqrt(variances) S &lt;- V %*% R %*% V cor(rmvnorm(100, sigma=S) )  # repeat this to get an idea of the variance around Pearson's estimator   Instance where variances are not equal to 1  V &lt;- diag(c(sqrt(3), sqrt(2)))  S &lt;- V %*% R %*% V cor(rmvnorm(100, sigma=S) )   This seems to be correct, but I would like expert criticism.
19801
I'm trying to predict a binary outcome using 50 continuous explanatory variables (the range of most of the variables is $-\infty$ to $\infty$). My data set has almost 24,000 rows. When I run glm in R, I get:  Warning messages:   1: glm.fit: algorithm did not converge   2: glm.fit: fitted probabilities numerically 0 or 1 occurred    I've read the other responses that suggest perfect separation might be occurring, but I'm confident that isn't the case in my data (though quasi-complete separation could exist; how can I test to see if that's the case?). If I remove some variables, the "did not converge" error might go away. But that's not always what happens.  I tried using the same variables in a bayesglm function and got the same errors.  What steps would you take to figure out exactly what's going on here? How do you figure out which variables are causing the problems?
103780
I observe a very weird behaviour in the SVD outcome of random data, which I can reproduce in both Matlab and R. It looks like some numerical issue in the LAPACK library; is it?  I draw $n=1000$ samples from the $k=2$ dimensional Gaussian with zero mean and identity covariance: $X\sim \mathcal N (0, \mathbf I)$. I assemble them in a $1000 \times 2$ data matrix $\mathbf X$. (I can optionally center $\mathbf X$ or not, it does not influence the following.) Then I perform singular value decomposition (SVD) to get $\mathbf X=\mathbf{USV}^\top$. Let's take some two particular elements of $\mathbf U$, e.g. $U_{11}$ and $U_{22}$, and ask what is the correlation between them across different draws of $\mathbf X$. I would expect that if the number $N_\mathrm{rep}$ of draws is reasonably big, then all such correlations should be around zero (i.e. population correlations should be zero, and sample correlations will be small).  However, I observe some weirdly strong correlations (around $\pm0.2$) between $U_{11}$, $U_{12}$, $U_{21}$, and $U_{22}$, and only between these elements. All other pairs of elements have correlations around zero, as expected. Here is how the correlation matrix for the $20$ "upper" elements of $\mathbf U$ looks like (first $10$ elements of the first column, then the first $10$ elements of the second column):     Notice strangely high values in the upper-left corners of each quadrant.  It was  that brought this effect to my attention. @whuber argued that the PC1 and PC2 are not independent and presented this strong correlation as an evidence for that. However, my impression is that he accidentally discovered a numerical bug in LAPACK library. What is going on here?  Here is @whuber's R code:  stat &lt;- function(x) {u &lt;- svd(x)$u; c(u[1,1], u[2, 2])}; Sigma &lt;- matrix(c(1,0,0,1), 2); sim &lt;- t(replicate(1e3, stat(MASS::mvrnorm(10, c(0,0), Sigma)))); cor.test(sim[,1], sim[,2]);   Here is my Matlab code:  clear all rng(7)  n = 1000;     %// Number of variables k = 2;        %// Number of observations Nrep = 1000;  %// Number of iterations (draws)  for rep = 1:Nrep     X = randn(n,k);     %// X = bsxfun(@minus, X, mean(X));     [U,S,V] = svd(X,0);      t(rep,:) = [U(1:10,1)' U(1:10,2)']; end  figure imagesc(corr(t), [-.5 .5]) axis square hold on plot(xlim, [10.5 10.5], 'k') plot([10.5 10.5], ylim, 'k')
93254
I'm looking for a method to calculate the area of overlap between two kernel density estimates in R, as a measure of similarity between two samples. To clarify, in the following example, I would need to quantify the area of the purplish overlapping region:  library(ggplot2) set.seed(1234) d &lt;- data.frame(variable=c(rep("a", 50), rep("b", 30)), value=c(rnorm(50), runif(30, 0, 3))) ggplot(d, aes(value, fill=variable)) + geom_density(alpha=.4, color=NA)     A similar question was discussed , the difference being that I need to do this for arbitrary empirical data rather than predefined normal distributions. The overlap package addresses this question, but apparently only for timestamp data, which doesn't work for me. The Bray-Curtis index (as implemented in vegan package's vegdist(method="bray") function) also seems relevant but again for somewhat different data.  I'm interested in both the theoretical approach and the R functions I might employ to implement it.
101676
Context: I want to store the TIFF data to use it for running a linear regression model later on. Ultimate goal: Create a carbon map using field data as a basis. Libraries: library(sp) library(raster) library(tidyverse) library(rgdal)  plotcoord &lt;- SpatialPointsDataFrame(AGCplot[,3:4], AGCplot[,1:2],proj4string=CRS('+proj=longlat'))  pspcoord&lt;-spTransform(plotcoord, CRS('+proj=utm +zone=21')) #Omzetten naar UTM21N (= S1) plot(pspcoord)  s1 &lt;- stack('Some.tif')  colnamen = c(&quot;plotid&quot;,&quot;AGC&quot;,&quot;B2&quot;,&quot;B3&quot;,&quot;B4&quot;,&quot;B5&quot;,&quot;B6&quot;,&quot;B7&quot;,&quot;B8&quot;,&quot;B8A&quot;,&quot;B11&quot;,&quot;B12&quot;,&quot;NDVI&quot;,&quot;NDWI&quot;,&quot;MMRI&quot;,&quot;IRECI&quot;,&quot;MASI&quot;,&quot;VV&quot;,&quot;VH&quot;)  s1_e &lt;- extract(s1,pspcoord,method='bilinear')  Upon extraction I get the following message: Warning message: In proj4string(x) : CRS object has comment, which is lost in output
45331
I have a line shapefile representing a road network. I wish to rasterize this data, with the resulting values in the raster showing the total length of the lines that fall within the raster cell.  The data is in the British National Grid projection so the units will be meters.   Ideally I would like to perform this operation using R, and i'm guessing that the rasterize function from the raster package would play a part in achieving this, I just can't work out what the function applied should be.