What I think would be the best way to do what you want in ITensor is to change the local basis used for your DMRG calculation. So instead of computing your ground state in a real-space basis as you are doing now, and then desiring to transform your wavefunction as a secondary step into a y-momentum basis, I think it would be both much easier and more efficient to work in a y-momentum (and of course x-real-space) basis from the very beginning, at the level of defining your site indices and Hamiltonian.
So the broad steps to do this would be as follows:
1.a. define a new site type that carries y-momentum quantum numbers, such as QN("Ky",0), QN("Ky",1), QN("Ky",-1), QN("Ky",2). Because the values of quantum numbers in ITensor are required to be integers, these integers are to be understood as QN("Ky",n) where the true momentum value is ky=(2pin)/Ny. So maybe a better name for these QN's would actually be QN("ny",0), QN("ny",1), etc. The name is really up to you.
1.b. to actually define such a new site type in the C++ version of ITensor unfortunately requires a bit of work and expertise with how to define a C++ class and its constructor and other similar things. If you are comfortable with this, and understand the general pattern of the existing site types included in ITensor such as in itensor/mps/sites/spinhalf.h, itensor/mps/sites/electron.h, etc. then it should not be too hard to make a new site type (call it KySpinHalf, say) and use it in your code instead of the SpinHalf sites. A key step here is that you'll have to take advantage of the fact that the named Args passed to the site type constructors gets a value "SiteNumber" passed to it, and you'll have to also pass "Ny" to it, so that these can be combined to work out what QN("ny",n) value that site should have based on where it's located around the cylinder. I can help you with this as I know this is probably the trickiest and least clear part.
the next step is to transform your Hamiltonian to a hybrid y-momentum / x-real-space basis. This is just a standard excersise in changing variables to spin operators which are plane-wave linear combinations of real-space spin operators by substituting them into your real-space Hamiltonian. I'm not saying "standard" to imply it's obvious or easy: just that this part is a thing which is purely mathematical and outside of ITensor.
You can use the newly created site types along with AutoMPO to input your momentum-space Hamiltonian and get an MPO for it. It might be a slightly bigger MPO than a purely real-space one because now the spin interactions will be non-local in the y-direction around each ring or rung of the cylinder. But the extra savings coming from the block-sparsity of the momentum-conserving ITensors should easily make up for this extra cost.
Finally, and we may have to discuss the particulars of this more, you can SVD your ground state MPS after it's optimized to get the entanglement spectrum, except now the bond indices of your MPS will automatically be labeled by the QN("ny",n) quantum numbers. So by analyzing the spectrum along with these quantum number labels, you can work out which sector the different singular values correspond to. Basically this could be as easy as SVDing part of your MPS and just printing out the "S" tensor of singular values. But we can explore this in some follow-up questions if it turns out to be more complicated in practice for some reason.
Hope that helps you to get started. One of the postdocs at Flatiron CCQ already implemented many of the steps above and it worked well, so if you email me I could put you in touch with him to see if he's willing to share some tips or perhaps even some code examples. We should ultimately put this into the library or at least the Code Formulas page of the website.