I have two sets of 2d-data points which can be transformed in each other by using a certain transformation function (see also this question).

The problem is that in both data sets there are points which do not have associated points in the other list.

data1 = 
{
{33.9168, 631.006}, {48.8067, 673.198}, {3.59394, 671.167},
{64.1931, 632.506}, {58.7559, 613.401}, {5.45129, 635.602}, 
{40, 500}, {55.6619, 651.298}, {40, 850}, {18.1513, 671.949},
{54.6781, 598.251}, {23.5348, 608.289}, {65.0549, 531.442},
{74.4132, 479.425}, {32.9808, 671.931}, {46.4516, 750.192},
{26.9262, 650.35}, {27.1816, 413.334}, {20.3858, 633.391}, 
{50.9284, 770.49}, {64.1628, 670.801}, {13.1805, 652.588}, 
{41.4876, 650.752}, {82.9996, 514.631}, {36.0045, 612.007},
{26.4914, 548.723}, {58.3295, 458.015}, {21.557, 801.607}, 
{5.84689, 800.425}
};

data2 = 
{
{1532.93, 536.587}, {1514.13, 789.}, {1530.22, 596.423}, 
{1520.66, 640.844}, {1540.5, 660.237}, {1530.03, 790.2}, 
{1559.17, 758.9}, {1556.15, 661.154}, {1580.39, 467.111}, 
{1525.63, 660.167}, {1571.44, 620.556}, {1512.62, 623.985}, 
{1520, 500}, {1533.79, 638.607}, {1526.88, 621.69}, 
{1560.9, 586.053}, {1572.13, 658.656}, {1548.37, 638.933}, 
{1532.8, 400.935}, {1540.44, 618.794}, {1590.15, 501.882}, 
{1554.5, 738.5}, {1564.73, 445.615}, {1543.06, 600.093}, 
{1565.69, 601.532}, {1562.55, 639.132}, {1511.34, 659.395}, 
{1580, 400}, {1585, 700}, {1571.9, 519.25}
};

In the upper plots I have marked these particular points.

Question

How can I remove these marked non-associated points and then sort the remaining points in such a way that the first new point of the corrected data1 list (data1corrected) is corresponding to the new first point of the corrected data2 list (data2corrected) and so on.

I wish then to receive e.g. for the corrected lists:

data1corrected=
{
{21.557, 801.607}, {5.84689, 800.425}, {50.9284, 770.49}, 
{46.4516, 750.192}, {32.9808, 671.931}, {48.8067, 673.198}, 
{3.59394, 671.167}, {18.1513, 671.949}, {64.1628, 670.801}, 
{13.1805, 652.588}, {55.6619, 651.298}, {26.9262, 650.35}, 
{41.4876, 650.752}, {5.45129, 635.602}, {20.3858, 633.391}, 
{64.1931, 632.506}, {33.9168, 631.006}, {58.7559, 613.401}, 
{36.0045, 612.007}, {23.5348, 608.289}, {54.6781, 598.251}, 
{26.4914, 548.723}, {65.0549, 531.442}, {82.9996, 514.631}, 
{74.4132, 479.425}, {58.3295, 458.015}, {27.1816, 413.334}
};

data2corrected=
{
{1530.03, 790.2}, {1514.13, 789.}, {1559.17, 758.9}, 
{1554.5, 738.5}, {1540.5, 660.237}, {1556.15, 661.154}, 
{1511.34, 659.395}, {1525.63, 660.167}, {1572.13, 658.656}, 
{1520.66, 640.844}, {1562.55, 639.132}, {1533.79, 638.607},     
{1548.37, 638.933}, {1512.62, 623.985}, {1526.88, 621.69},   
{1571.44, 620.556}, {1540.44, 618.794}, {1565.69, 601.532}, 
{1543.06, 600.093}, {1530.22, 596.423}, {1560.9, 586.053}, 
{1532.93, 536.587}, {1571.9, 519.25}, {1590.15, 501.882}, 
{1580.39, 467.111}, {1564.73, 445.615}, {1532.8, 400.935}
};

Then I could use the following transformation:

trafo = 
  FindGeometricTransform[data2corrected, data1corrected, TransformationClass -> "Affine"];

which gives for each point {x, y} the coordinate transormation:

trafo[[2]][{x, y}]

{1502.56 + 1.00443 x + 0.0069496 y, -13.1299 - 0.00592912 x + 1.00225 y}

Applying this transformation to data1corrected results in:

data1correctedtransformed= trafo[[2]][data1corrected];

And finally I can plot the transormed coordinates of data1correctedtransformed together with the original points data2corrected:

ListPlot[{data2corrected, data1correctedtransformed}, 
 PlotRange -> {{1510, 1600}, {300, 900}}, PlotStyle -> {Blue, Green}, 
 Frame -> True, FrameLabel -> {{"y", ""}, {"x", 
    "data1correctedtransformed (Green), data2corrected (Blue)"}}, 
 BaseStyle -> {FontWeight -> "Bold", FontSize -> 15, 
   FontFamily -> "Calibri"}, ImageSize -> Large]