Timings

ClearAll[timings]
timings[m_Integer] := Module[{t1, t2, t3, t4, n = 1000, a, b, c, b2a, b2b, b2c, b2d, sa}, 
   SeedRandom[1];
   {a, b} = RandomReal[1, {2, n, n}];
   c = RandomSample[Tuples[Range[n], {2}], m];
   b2a = b2b = b2c = b2d = b;
   t1 = First[AbsoluteTiming[sa = SparseArray[c -> 1, {n, n}];  
   b2a = (1 - sa) b2a + sa a;]];
   t2 = First[AbsoluteTiming[(b2b[[##]] = a[[##]]) & @@@ c;]];
   t3 = First[AbsoluteTiming[(Part[b2c, #[[1]], #[[2]]] = Part[a, #[[1]], #[[2]]])&/@c;]];
   t4 = First[AbsoluteTiming[b2d = ReplacePart[b2d, Thread[c -> Extract[a, c]]];]]; 
   {m, t1,  t2, t3, t4, b2a == b2b == b2c == b2d}]

tr = timings /@ {100, 500, 5000, 50000, 100000, 300000};
backgroundrule = # -> Item[#, Background -> Yellow] & /@ Min /@ tr[[All, 2 ;; 5]];
headers = {{"Length@desired", "method", SpanFromLeft, SpanFromLeft, 
  SpanFromLeft, " b2a==b2b==b2c==b2d "}, {SpanFromAbove,   
  " SparseArray ", " ApplySetPart ", " MapSetPart ", " ReplacePart ", SpanFromAbove}};
Grid[Join[headers, timingresults /. backgroundrule], 
  Dividers -> All,  Alignment -> {Center, Center}] 

SeedRandom[1]
n = 10;
array = ImageData[RandomImage[1, {n, n}]];
array2 = ImageData[RandomImage[1, {n, n}]];


sa = SparseArray[desired -> 1, Dimensions@array];
newarray = (1 - sa) array2 + sa array;

Showing only the desired parts of the three matrices:

Row[Labeled[MatrixPlot[MapAt[Black &, #, Complement[Tuples[Range[n], {2}], desired]], 
     ImageSize -> 200], #2, Top] & @@@ 
  Transpose[{{array, array2, array2a}, {"array", "array2", "array2a"}}]]

All four methods suggested in answers and comments so far give the same result:

array2a = array2b = array2c = array2d = array2;
array2a = (1 - sa) array2 + sa array;
(Part[array2b, #[[1]], #[[2]]] = Part[array, #[[1]], #[[2]]]) & /@ desired; (* user6014 *)
(array2c[[##]] = array[[##]]) & @@@ desired; (* Albert Retey's comment *)
array2d = ReplacePart[array2d, 
  Thread[desired -> (array[[##]] & @@@ desired)]]; (* J.M.'s comment *)

 array2a == array2b == array2c == array2d

True

SeedRandom[1]
n = 10;
array = ImageData[RandomImage[1, {n, n}]];
array2 = ImageData[RandomImage[1, {n, n}]];


sa = SparseArray[desired -> 1, Dimensions@array];
newarray = (1 - sa) array2 + sa array;

Showing only the desired parts of the three matrices:

Row[Labeled[MatrixPlot[MapAt[Black &, #, Complement[Tuples[Range[n], {2}], desired]], 
     ImageSize -> 200], #2, Top] & @@@ 
  Transpose[{{array, array2, array2a}, {"array", "array2", "array2a"}}]]

All four methods suggested in answers and comments so far give the same result:

array2a = array2b = array2c = array2d = array2;
array2a = (1 - sa) array2 + sa array;
(Part[array2b, #[[1]], #[[2]]] = Part[array, #[[1]], #[[2]]]) & /@ desired; (* user6014 *)
(array2c[[##]] = array[[##]]) & @@@ desired; (* Albert Retey's comment *)
array2d = ReplacePart[array2d, 
  Thread[desired -> Extract[array, desired]]]; (* J.M.'s comment *)

 array2a == array2b == array2c == array2d

True