We present the first direct evaluation of Delta I = 3/2K -> pi pi matrix elements with the aim of determining all the low-energy constants at NLO in the chiral expansion. Our numerical investigation demonstrates that it is indeed possible to determine the K -> pi pi matrix elements directly for the masses and momenta used in the simulation with good precision. In this range however, we find that the matrix elements do not satisfy the predictions of NLO chiral perturbation theory. For the chiral extrapolation we therefore use a hybrid procedure which combines the observed polynomial behavior in masses and momenta of our lattice results, with NLO chiral perturbation theory at lower masses. In this way we find stable results for the quenched matrix elements of the electroweak penguin operators ((1=2) <pi pi vertical bar O(8)vertical bar K-0> = (0.68 +/- 0.09) GeV3 and (1=2) <pi pi vertical bar O(7)vertical bar K-0> = (0.12 +/- 0.02) GeV3 in the NDR-(MS) over bar scheme at the scale 2 GeV), but not for the matrix elements Of 04 (for which there are too many low-energy constants at NLO for a reliable extrapolation). For all three operators we find that the effect of including the NLO corrections is significant (typically about 30%). We present a detailed discussion of the status of the prospects for the reduction of the systematic uncertainties.