How to rearrange 2D data to get given correlation? I have the following simple dataset with two continuous variables; i.e.:
d = data.frame(x=runif(100,0,100),y = runif(100,0,100))
plot(d$x,d$y)
abline(lm(y~x,d), col="red")
cor(d$x,d$y) # = 0.2135273


I need to rearrange the data in the way to have correlation between variables to be ~0.6. I need to keep means and other descriptive statistics (sd,min,max,etc.) of both variables constant. 
I know it is possible to make almost any correlation with the given data i.e.:
d2 = with(d,data.frame(x=sort(x),y=sort(y)))
plot(d2$x,d2$y)
abline(lm(y~x,d2), col="red")
cor(d2$x,d2$y) # i.e. 0.9965585


If I try to use sample function for this task:
cor.results = c()
for(i in 1:1000){
    set.seed(i)
    d3 = with(d,data.frame(x=sample(x),y=sample(y)))
    cor.results =  c(cor.results,cor(d3$x,d3$y))
}

I get quite wide range of correlations:
> summary(cor.results)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-0.281600 -0.038330 -0.002498 -0.001506  0.034380  0.288800

but this range depends on number of rows in data frame and decreasing with increase of size.
> d = data.frame(x=runif(1000,0,100),y = runif(1000,0,100))
> cor.results = c()
> for(i in 1:1000){
+ set.seed(i)
+ d3 = with(d,data.frame(x=sample(x),y=sample(y)))
+ cor.results =  c(cor.results,cor(d3$x,d3$y))
+ }
> summary(cor.results)
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.1030000 -0.0231300 -0.0005248 -0.0005547  0.0207000  0.1095000

My question is:
How to rearrange such dataset to get given correlation (i.e. 0.7)?
(It will be also good if method will remove dependence on dataset size) 
 A: Here is one way to rearrange the data that is based on generating additional random numbers.
We draw samples from a bivariate normal distribution with specified correlation. Next, we compute the ranks of the $x$ and $y$ values we obtain. These ranks are used to order the original values. For this approach, we have top sort both the original $x$ and $y$ values.
First, we create the actual data set (like in your example).
set.seed(1)
d <- data.frame(x = runif(100, 0, 100), y = runif(100, 0, 100))

cor(d$x, d$y)
# [1] 0.01703215

Now, we specify a correlation matrix.
corr <- 0.7  # target correlation
corr_mat <- matrix(corr, ncol = 2, nrow = 2)
diag(corr_mat) <- 1
corr_mat
#      [,1] [,2]
# [1,]  1.0  0.7
# [2,]  0.7  1.0

We generate random data following a bivariate normal distribution with $\mu = 0$, $\sigma = 1$ (for both variables) and the specified correlation. In R, this can be done with the mvrnorm function from the MASS package. We use empirical = TRUE to indicate that the correlation is the empirical correlation (not the population correlation).
library(MASS)
mvdat <- mvrnorm(n = nrow(d), mu = c(0, 0), Sigma = corr_mat, empirical = TRUE)

cor(mvdat)
#      [,1] [,2]
# [1,]  1.0  0.7
# [2,]  0.7  1.0

The random data perfectly matches the specified correlation.
Next, we compute the ranks of the random data.
rx <- rank(mvdat[ , 1], ties.method = "first")
ry <- rank(mvdat[ , 2], ties.method = "first")

To use the ranks for the original data in d, we have to sort the original data.
dx_sorted <- sort(d$x)
dy_sorted <- sort(d$y)

Now, we can use the ranks to specify the order of the sorted data.
cor(dx_sorted[rx], dy_sorted[ry])
# [1] 0.6868986

The obtained correlation does not perfectly match the specified one, but the difference is relatively small.
Here, dx_sorted[rx] and dy_sorted[ry] are resampled versions of the original data in d.
A: To generate two uniform distributions with a specified correlation, the Ruscio & Kaczetow (2008) algorithm will work. They provide R code.  You can then transform with a simple linear function to get your target min, max, mean, and SD.
Ruscio & Kaczetow Algorithm
I'll summarize the bivariate case, but it can also work with multivariate problems.  Uncorrelated $X_o$ and $Y_o$ are generated with any shape (e.g., uniform). Then, $X_1$ and $Y_1$ are generated as bivariate normal with an intermediate correlation. $X_1$ and $Y_1$ are replaced by $X_0$ and $Y_0$ in a rank-preserving fashion. Adjust the intermediate correlation to be higher or lower depending on whether the r($X_1,Y_1$) is too low or too high. $X_2$ and $Y_2$ are generated as bivariate normal with the new intermediate correlation.  Repeat.
Notice that this is very similar to @Sven Hohenstein's solution, except that it's iterative, so the intermediate correlation will get closer and closer to the target correlation until they are indistinguishable.  Also, note that this algorithm can be used to generate a large population (e.g., N=1 million) from which to draw smaller samples - that is useful if you need to have sampling error.
For a related post: Correlation and non-normal distributions
Preserving Descriptive Statistics
There is no guarantee that the algorithm will produce the exact same descriptives.  However, because a uniform distribution's mean and SD are determined by its min and max, you can simply adjust the min and max to fix everything.  
Let $X_g$ and $Y_g$ be your generated variables from the last iteration of the Ruscio & Kaczetow algorithm, $X_f$ and $Y_f$ be your final variables that you hope to have (with target descriptives), and $X$ and $Y$ be your the original variables in your dataset.  
Calculate  $X_f=(X_g - min(X))*(max(X)-min(x))/(max(X_g)-min(X_g))$
Do the same for $Y_f$
Reference:
Ruscio, J., & Kaczetow, W. (2008). Simulating multivariate nonnormal data using an iterative algorithm. Multivariate Behavioral Research, 43, 355–381. doi:10.1080/00273170802285693
A: I'm guessing that when you say "resample" you mean "simulate," which is more general.  The following is the most concise way I know to simulate normal, bivariate data with a specified correlation.  Substitute your own desired values for r and n.
r = .6
n = 1000
x = rnorm(n) 
z = rnorm(n) 
y = (r/(1-r^2)^.5)*x + z

cor(x,y)
plot(x,y)
abline(lm(y~x), col="red")

