WorldWideTelescope / toasty

Break large images into "tile pyramids", with a focus on the all-sky TOAST format.
https://toasty.readthedocs.io/
MIT License
0 stars 4 forks source link

Toasted image doesn't show in wwt #100

Closed keflavich closed 11 months ago

keflavich commented 11 months ago

I'm trying to use toasty for the first time, so I might be doing things wrong, but I can't get a toasted wtml to show up

This: https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_Feather/index_rel.wtml is produced with:

from toasty import study, image as timage, pyramid, builder
from wwt_data_formats import write_xml_doc

tim = timage.Image.from_pil(img)

bui = builder.Builder(pyramid.PyramidIO('/orange/adamginsburg/web/public/ACES/MUSTANG_Feather'))
stud = bui.prepare_study_tiling(tim)
height, width, _ = np.array(img).shape
bui.apply_wcs_info(mustang_fk5, width=width, height=height)
if not os.path.exists('/orange/adamginsburg/web/public/ACES/MUSTANG_Feather/7/44/44_27.png'):
    # to avoid re-running every time I modify something
    bui.execute_study_tiling(tim, stud)
# commented out, but I tried this approach too, which gives the full URL
#bui.imgset.url = "http://data.rc.ufl.edu/web/public/adamginsburg/ACES/MUSTANG_Feather" + bui.imgset.url
fldr = bui.create_wtml_folder()
bui.write_index_rel_wtml()

When I try to load that wtml in WWT https://worldwidetelescope.org/webclient/?wtml=https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_Feather/index_rel.wtml I see nothing, but I see signs of problems in the browser console. It looks like wwt is asking for https://0.0.0.0/0/0_0.png, which looks like a failure to fill '{1}/{3}/{3}_{2}.png' correctly. It's picking the wrong tile, too, maybe?

How do I get my toasted files to show up?

pkgw commented 11 months ago

Hi Adam,

I think you're close. The error that you mention looks like what you might get when the WTML contains a relative URL to the image data, which (unfortunately-ish) isn't allowed in publicly-hosted WTML files.

Pulling down your WTML file right now, it looks like you've updated it? I'm seeing:

Url="http://data.rc.ufl.edu/web/public/adamginsburg/ACES/MUSTANG_Feather/{1}/{3}/{3}_{2}.png"

for the image data URL, which looks more correct. But it's not quite right, since the following doesn't exist:

https://data.rc.ufl.edu/web/public/adamginsburg/ACES/MUSTANG_Feather/0/0/0_0.png

(i.e., if you replace {1} and {2} and {3} with 0, you should get a file.)

Snooping around your webspace, I see two problems. First, the prefix is wrong: it should be ...ufl.edu/pub/adamginsburg/....

Second, it looks like you've only upload the 7/ folder, which contains the highest-resolution tiles. There should be additional folders 6/, 5/, ... 0/ containing downsampled tiles, all the way down to 0/0/0_0.png which definitionally covers the entire image. If those you generated those files locally and just forgot to upload them, no problem.

If you didn't generate them locally, you need to run something like toasty cascade --start 7 . to generate them.

For the initial issue, Toasty has a convention where the index_rel.wtml can contain relative URLs, because we have a command called wwtdatatool preview that starts up a local server with on-the-fly WTML rewriting to allow local testing of your data and WTML. Then when you're happy and you know the final URL at which you'll publish the data, wwtdatatool wtml rewrite-urls can generate a rewritten file (conventionally called index.wtml) with the relative URLs rewritten as absolute ones.

keflavich commented 11 months ago

Thanks @pkgw.

The pub vs public is a typo, that I can fix.

The missing 0 level is the real problem here: when I ran toast, it didn't make that level. I'm using the python module, not the command-line tools; I'll see if I can reconstruct that cascade command.

pkgw commented 11 months ago

Gotcha. It should be pretty straightforward. Here's the CLI code:

https://github.com/WorldWideTelescope/toasty/blob/master/toasty/cli.py#L69

Relevant API docs: https://toasty.readthedocs.io/en/latest/api/toasty.merge.html

keflavich commented 11 months ago

Thanks. I'm very close now, just trying to work out the coordinates:

https://worldwidetelescope.org/webclient/?wtml=https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_Feather/index_rel.wtml#/

my attempt is:

# rotate MUSTANG WCS into fk5
mustang_fk5 = wcs.WCS(naxis=2)
mustang_fk5.wcs.ctype = 'RA---TAN', 'DEC--TAN'
# old version refpos = mustangwcs.pixel_to_world(*mustangwcs.wcs.crpix)
refpos = g00 = SkyCoord(0*u.deg, 0*u.deg, frame='galactic') 
mustang_fk5.wcs.crval = refpos.fk5.ra.deg, refpos.fk5.dec.deg
mustang_fk5.wcs.crpix = mustangwcs.world_to_pixel(g00)

offrefpos = mustangwcs.pixel_to_world(mustang_fk5.wcs.crpix[0] + 5, mustang_fk5.wcs.crpix[1])
# debug posang = -coordinates.position_angle(refpos.fk5.ra, refpos.fk5.dec, offrefpos.fk5.ra, offrefpos.fk5.dec)
# unclear why I'm 90 off...
posang = 90*u.deg - coordinates.position_angle(refpos.fk5.ra, refpos.fk5.dec, offrefpos.fk5.ra, offrefpos.fk5.dec)
print(posang.deg)
pixscale = np.mean(np.abs(mustangwcs.wcs.cdelt)) #<--- this seems to be the problem?
mustang_fk5.wcs.cd = np.array([[np.cos(posang), -np.sin(posang)], [np.sin(posang), np.cos(posang)]]) * pixscale

bui = builder.Builder(pyramid.PyramidIO('/orange/adamginsburg/web/public/ACES/MUSTANG_Feather'))
stud = bui.prepare_study_tiling(tim)
height, width, _ = np.array(img).shape
if not os.path.exists('/orange/adamginsburg/web/public/ACES/MUSTANG_Feather/7/44/44_27.png'):
    bui.execute_study_tiling(tim, stud)
if not os.path.exists('/orange/adamginsburg/web/public/ACES/MUSTANG_Feather/0/0/0_0.png'):
    merge.cascade_images(
        bui.pio, start=7, merger=merge.averaging_merger, cli_progress=False
    )
bui.imgset.url = "https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_Feather/" + bui.imgset.url
bui.apply_wcs_info(mustang_fk5, width=width, height=height)
fldr = bui.create_wtml_folder()
bui.write_index_rel_wtml()
with open(os.path.join(bui.pio._base_dir, "index.wtml"), 'w') as fh:
    write_xml_doc(fldr.to_xml(), dest_stream=fh)

I'm.... pretty sure the WCS itself is right, but it's off-center. Maybe this is one of those fits-vs-png issues....

pkgw commented 11 months ago

If you're manually applying FITS WCS to a derived PNG image, that might be a factor. Here's some more relevant code from the Toasty CLI:

https://github.com/WorldWideTelescope/toasty/blob/master/toasty/cli.py#L614-L673

The Toasty image-handling framework has a utility function ensure_negative_parity() that, well, ensures that a WCS coordinate system has negative parity, which is what you need for tiled images in WWT. See:

https://github.com/WorldWideTelescope/toasty/blob/master/toasty/image.py#L244 — computing parity sign https://github.com/WorldWideTelescope/toasty/blob/master/toasty/image.py#L259 — flipping WCS parity

The way I remember this stuff is that if you have an image that lands on the sky properly (say, a FITS with positive parity), if you flip the WCS parity in the way implemented above and flip the image bitmap top-to-bottom, nothing changes. And since a naive FITS-to-PNG conversion does that top-to-bottom flip, you usually need to apply a WCS parity flip.

keflavich commented 11 months ago

My pixel scale is also set incorrectly, so that must be related to the scales in those snippets.

Is this forcing the pixel scale to 1 deg/pix?

            wcs.wcs.crpix *= scale

            if hasattr(wcs.wcs, "cd"):
                wcs.wcs.cd /= scale
            else:
                wcs.wcs.cdelt /= scale

EDIT: No, it's not.

keflavich commented 11 months ago

So I think this line: https://github.com/WorldWideTelescope/wwt_data_formats/blob/5a95ef5c6209a1170b5f4a4c5ada6e58020af398/wwt_data_formats/imageset.py#L546 is setting the scale incorrectly, but I'm not sure why.

It may just be that I did something dumb.

keflavich commented 11 months ago

yeah, my mistake was using the wrong wcs 🤦

nevertheless, I'm still off by a bit - might be the translation from Galactic to fk5. Any tips on how best to do that? https://worldwidetelescope.org/webclient/?wtml=https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_Feather/index_rel.wtml#/

pkgw commented 11 months ago

For reference, the line that you pointed out has to do with the way that the WWT data formats interpret the same field in different ways depending on the image mode (docs).

What's your definition of "off by a bit"? The default DSS background map has an astrometric offset at the 1 arcsec level (which I would dearly love to fix one day, sigh).

I could see how a galactic/equatorial translation could induce issues, though. I'm pretty sure that the only way to avoid issues there would be to resample the original data onto a TAN equatorial system, and then tile that resampled image. The advantage of this approach is that if you can do that resampling, then you can probably reuse your current code more-or-less as-is.

Alternatively, you could tile the image as-is into a TOAST pixelization, rather than a tiled TAN. The TOAST format is mainly used for all-sky maps, but there's no reason it can't be used for relatively small images — toasty has infrastructure to filter the processing down to only the areas where there's actual image data. Not only does TOAST render more precisely than TAN over large-ish angular scales (degree-ish and larger), but the way that the tiling works, any kind of wacky WCS in the input image will be handled correctly (basically, instead of pre-resampling as above, the resampling happens on the fly). This approach should give the best results, but would require a bit of reworking of your code. I have a sample script that I use for DASCH that you could build off of, though.

(BTW, if you have any interest in optical photometry over 100-year timescales, I just refreshed the DASCH website and will be making it a lot better in the coming months! And we have a new API coming that I am going to ask you about integrating into astroquery (i.e. I'm up for writing all of the code, but want to do it according to your design preferences).)

keflavich commented 11 months ago

What's your definition of "off by a bit"? The default DSS background map has an astrometric offset at the 1 arcsec level (which I would dearly love to fix one day, sigh).

image

Sgr B2, The Brick, and Sgr A* are all up and to the right by ~a few arcminutes. This looks like a linear shift in Galactic lat, which is surprising.

Reprojecting into TAN is something I very strongly want to avoid. It requires duplicating data and resampling into a much more sparsely populated projection. I'm treating that as a last resort, and possibly even after last - I'm going through all this effort to figure out how I can show Galactic images in WWT more generally.

If you have a good example of how to project into TOAST, I'd be happy to take that approach - I'm only following the current approach because it's what I could figure out.

(I saw your blog post about DASCH on my RSS feed - that's cool. Post something on astroquery issues about that)

pkgw commented 11 months ago

Do you have a sample source file that I could play around with? I think that would help me try to get a handle on whether a reprojection can be avoided, and/or test out TOASTing code.

keflavich commented 11 months ago

https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_12m_feather_noaxes.png has embedded AVM

keflavich commented 11 months ago

OK, after further investigation, it's a problem with my translation from Galactic -> fk5.

This is something that Should Be Easy but clearly isn't. Or maybe it's just me.

keflavich commented 11 months ago

A solution that works is to just fit the WCS:

from astropy.wcs.utils import fit_wcs_from_points

points = np.array([(0,0), (0, rslt.data.shape[0]), (rslt.data.shape[1], 0), (rslt.data.shape[1], rslt.data.shape[0])])
wpoints = ACESwcs.pixel_to_world(points[:,0], points[:,1])
ACES_fk5 = fit_wcs_from_points((points[:,0], points[:,1]), wpoints.fk5)

However, when I try this, the parity ends up wrong.... so I can force it to be inverted:

# really dumb, definitely wrong hack to make the builder not complain about parity
ACES_fk5.wcs.cd[:,1] *= -1

This works better, but it's still offset.

https://worldwidetelescope.org/webclient/?wtml=https://data.rc.ufl.edu/pub/adamginsburg/ACES/MUSTANG_Feather/index_rel.wtml#/

image
keflavich commented 11 months ago

OK, I solved the offset: it was the classic FITS vs python indexing issue.

Anyone coming to this thread late: the only toast issue was the first thing Peter reported, which was the missing 0/0_0/0.png caused by running the wrong tiler.

pkgw commented 11 months ago

Ah, yes, I'm familiar with the feeling of slogging through a series of WCS issues manifesting as smaller and smaller coordinate offsets ...

Would it potentially make sense to take what you've worked out and adapt it into some utility function to help people with the galactic → equatorial translation? Clearly it's non-trivial to get the details right.

Other than that, is it time to close this issue?

keflavich commented 11 months ago

from astropy.wcs.utils import fit_wcs_from_points is a straightforward enough solution, it's just hacky.

But yes, closing.