Radio Galaxy Zoo Talk

I'm thinking of making the 'FIRST contour overlay' Python code I wrote available to any zooite ...

  • JeanTate by JeanTate

    The thread Suggestions for RGZ Objects to show with detailed FIRST contours overlaid on SDSS images, started on July 5 2014 by WizardHowl, is quite popular. And in the past year+ I have created thousands of these overlays, many of which I've posted here in RGZ Talk.

    The boilerplate footer I created, which I use for the overlays I post here, is something like "Boilerplate: SDSS image per http://skyservice.pha.jhu.edu/DR12/ImgCutout/getjpeg.aspx, FIRST (red) and NVSS (cyan) contours derived from FITS files produced using SkyView with Python code described in this RGZ Talk thread. Image center per the ARG image ({link}; J2000)."

    As you will have noticed, I'm falling increasingly behind in creating, and posting, the overlays (based on requests made in that thread).

    Perhaps it's time that I made my method, and code, available to any zooite? Perhaps there are zooites who would like to learn how to create such overlays, and would like to start to learn by using my code?

    What do you think?

    Posted

  • sisifolibre by sisifolibre

    I would like to learn how to make your overlaids, but I haven't any knowledge about programing, etc... Any case I'm sure that if you make your method available to any zoonite someone else will learn how to make it.

    Posted

  • JeanTate by JeanTate in response to sisifolibre's comment.

    Thanks sisifolibre, and those who PMed me on this topic.

    I plan to do this:

    • in each new post in this thread, I'll introduce one step in the process I have used
    • I'll aim to start with how to do a single contour overlay image; later I'll talk about how to do many in one go (batch mode)
    • When I've finished, I'll bundle everything into a single, new thread, and ask someone to "sticky" it
    • I do hope that everyone who's interested will post questions, comments, etc, as we go along! 😃

    Comments?

    Next post: an overview of the steps, before starting on the first one.

    Posted

  • Dolorous_Edd by Dolorous_Edd

    I also would be interested

    I think you should add a section about Python setup, that is a topic by itself

    Posted

  • JeanTate by JeanTate

    Next post: an overview of the steps, before starting on the first one.

    To produce a contour overlay image - FIRST on SDSS - you need to do these things:

    1. produce the SDSS 'canvas' that the contours will be overlaid onto
    2. obtain the FIRST data, from which the contours will be made
    3. create the contours
    4. overlay the contours onto the canvass
    5. produce the contour overlay image
    6. save it, publish it, etc

    Steps 1 and 2 can be done independently; you can produce the canvass first, or you can get the FIRST data first. 😉

    The Python code I developed will do steps 3 to 5.

    I also developed code to automate steps 1 and 2, or at least semi-automate them. The 'semi-automated' code for step 2 does not use Python. And I also semi-automated steps 3 to 5 in my Python code, so I could produce a lot of images, in 'batch mode'. However, I found that there will always be cases where things fail. For a variety of reasons. So you may need to do some intensive manual tweaking, in perhaps as many as 5% of cases.

    In these posts, I will start by going through how to produce a single FIRST contour overlay image on SDSS, 'manually'. Later, I will describe how I semi-automated the processes.

    In addition to Python (and what's needed for semi-automation of step 2), I use(d) a browser, a spreadsheet, and an image processing app. My choice? Firefox, OpenOffice, and GIMP, all of which are free (I have a Windows 7 machine). I'm pretty sure you can do everything I'll be describing with any other (sensible!) choice of browser, spreadsheet, and image processing app.

    Next: how to manually produce an SDSS canvass.


    And some background etc.

    I should stress that what I'll be outlining in these posts is the way I developed a method for producing contour overlay images. It is almost certainly not how a professional would have done it (astronomer, person with computer science training, etc). And if you decide to work on your own, you may well develop something different again. In a later post I'll say something about other approaches/methods I've come across.

    Also, while I'll be describing 'FIRST contour overlays on SDSS', my code will also produce NVSS overlays on SDSS (and both FIRST and NVSS). And it can be tweaked to use any canvass, or any WISE band contours, or ... I'll describe these things at the end.


    And it goes without saying: comments, questions, suggestions, ... feedback! is always welcome. 😄

    Posted

  • JeanTate by JeanTate

    Next: how to manually produce an SDSS canvass.

    I'll illustrate with an example, from the ARG0001ly0 thread: a field centered on SDSS J093139.03+320400.1.

    Here's how it appears in that thread:

    enter image description here

    And here's the URL for that image, sans the http:// part:

    skyservice.pha.jhu.edu/DR12/ImgCutout/getjpeg.aspx?ra=142.912666110626&dec=32.0667050436208&scale=0.2&width=200&height=200&opt=G

    It's easy to play with this, to produce different images, e.g. without the grid you simply leave off the "&opt=G" part:

    enter image description here

    And to make it 600x600 pix instead of 200x200, just change the "width=200&height=200" to "width=600&height=600":

    enter image description here

    There are three other parts of the URL that we need, "ra=", "dec=", and "scale=":

    1. "ra" and "dec" are, as you probably already know, Right Ascension and Declination, the celestial coordinates of the center of the field; they correspond to the Earthly longitude and latitude (Wikipedia on this). In making an SDSS canvass, you need only four significant digits in each (and sometimes three will do, at high Dec); here's the original with "ra=142.9127&dec=32.0667", can you see the difference?

    enter image description here

    1. "scale" is indeed a scale, in units of arcsecs per pixel, so smaller values give a zoomed-in image, larger a zoomed-out one; here are "scale=0.1" and scale=0.3":

    enter image description here enter image description here

    I settled on 600x600 pix images, with a scale of 0.24, for most of my FIRST contour overlay images. Why? Because:

    • they look pretty
    • most FIRST radio emission associated with a host of interest is smaller than ~2.4' across
    • 2.4' (2.4 arcmins) = 600 x 0.24" (0.24 arcsecs), so converting in my head was easy.

    So, how to capture/create a 600x600 pix SDSS canvass, centered on (RA, Dec), and 2.4'x2.4' in 'sky size'?

    • I used the SDSS DR12 Navigate tool, plugged in my RA and Dec values, and pressed the Search button
    • If I've done this correctly (i.e. the host I want is in the center), I then click on the Explore link, in the right-hand panel
    • With my mouse over the Explore image, I right-click, and select "Copy Image Location" (may be different if you use something other than Firefox)
    • I head on over to RGZ Talk, open a new post, post the image, but with this 'tail string': "scale=0.24&width=600&height=600"
    • after checking that this is what I want, I right-click, and select "Save Image as ..."
    • give it a (file) name, making sure to add ".jpg"
    • delete the post (no need to clutter up RGZ Talk)
    • Done! 😃

    Next: How to obtain the FIRST data, from which the contours will be produced


    Some background etc:

    In a later post I'll go over some hard-won Tips and Traps, especially with regard to semi-automating this; it is really time consuming to produce even 20 canvasses this way, let alone 200 (say).

    The JPG file you've saved is dumb, it contains no information on its size, orientation, location, etc. You need to keep/record this information, in order to ensure that the FIRST contours match the size, orientation, and location of the SDSS canvass. I'll be explaining how I did this in later posts.

    The SDSS DR12 JPG images, from the SDSS ImgCutout service on its SkyServer website, are "Luptonized". They combine data from the g, r, and i band sensors on the SDSS camera (the other two bands, u and z, are rarely used, but the data is available if you ever want to use it). The mapping is chromatic, it preserves wavelength order: g is mapped to JPGs' B, r to G, and i to R. The gri bands are not even close to the blue, green, and red cones in our eyes; for starters, i is in the infrared. So the colors are artificial.

    "Luptonize" refers to a method for presenting (exaggerating) color in astronomical images (the paper is Lupton+ (2004) "Preparing Red-Green-Blue Images from CCD Data", and is available from ADS, here). It involves transforming the range of intensities to ensure as much detail is 'visible' in the 256 levels in a JPG (CCDs produce 12 or 14 bit intensity outputs; a JPG image has 8 bits), and turning color contrasts waaaay up (without this last trick, most astronomical images would be mostly shades of grey, with perhaps some pale pastel colors).


    And it goes without saying: comments, questions, suggestions, ... feedback! is always welcome. 😄

    Posted

  • JeanTate by JeanTate

    Next: How to obtain the FIRST data, from which the contours will be produced

    This is, perhaps, the most straight-forward step. If you want the FIRST data for just one object.

    Go to the SkyView Query Form. And fill it in. It doesn't matter which order you enter the inputs:

    • "Coordinates or Source": I always use "ddd.dddd, ±dd.dddd" format, so for the object in the last post, "142.9127, +32.0667" (the "+" is optional)
    • "SkyView Surveys": click on "Radio"-> "VLA FIRST (1.4 GHz)"; make sure this is the ONLY one checked
    • "Common Options" -> "Image size (pixels)": enter "600"
    • "Common Options" -> "Image size (degrees)": enter "0.04"; this is another place where choosing 600 pixels, and scale 0.24"/pix as the size helps, dividing 2.4 by 600 is easily done in your head!
    • NOTE: if your SDSS canvass was not 2.4'x2.4' and 600x600 pix, you'll need to enter the values appropriate to that canvass
    • Click on the button "Submit Request"
    • a new tab (or window) will open
    • scroll down to "Download FITS"
    • right-click on "FITS"
    • select "Save As ...": you should have a new window, with the choice of where to save the file, its name, and file type ("fits File" for me)
    • make any changes you want, but make sure the file type is "FITS" (or "fits"), and note the filename and location
    • Save
    • close the SkyView windows (optional)
    • Done! 😃

    Next: how to create the contours.


    Some background etc:

    In a later post I'll describe how you can (semi-) automate this process. As with creating the SDSS canvass, it gets tedious if you have to collect the FIRST data for 20 sources, let alone 200. Not to mention error-prone.

    FITS stands for Flexible Image Transport System, and is widely used by astronomers. For our purposes, a key feature is that the location of the center of the image, size (in pix), scale, and orientation are all nicely stored in the file 'header'. And that the keywords for each are standard.

    It's a quite rich file format, and can be used for things other than images (and one file can contain several, linked, images).

    The NASA page A Primer on the FITS Data Format is a good place to explore further. There's also a NASA page called "FITS Image Software Packages for image viewing, analysis, and format conversion" (link). You can spend many hours exploring there; whether you end up thinking those hours were happily spent, or not, ... well, be warned that you may be in for some frustration (a lot of the packages do not have 'ordinary zooite friendly' user manuals! 😃


    And it goes without saying: comments, questions, suggestions, ... feedback! is always welcome. 😄

    Posted

  • ivywong by ivywong scientist, admin

    This is a great idea @JeanTate!. Perhaps you might like to start an iPython notebook ? Don't worry if you do not have the time, it's a passing thought.

    Posted

  • JeanTate by JeanTate in response to ivywong's comment.

    Thanks Ivy.

    Once I'm ~done with going through the basics I'll take a look into what an iPython notebook is (I'd never heard of such a thing until I read your post 😮), and maybe start one 😃

    Posted

  • JeanTate by JeanTate in response to ivywong's comment.

    Next: how to create the contours.

    This is the first step which requires Python. I am particularly unqualified to explain how to install Python, and the basics of how to code in it, etc, so I'm not going to try. Perhaps we could get a discussion going on this topic, in a separate thread? For the record, I got started on Python with Ivezić et al.'s "Statistics, Data Mining, and Machine Learning in Astronomy: A Practical Python Guide for the Analysis of Survey Data" (link), a book I can highly recommend (though it's not cheap).

    In this post I'll include code snippets, which are not something you could get to work by simply copy/pasting; the full code will come in a later post

    The FIRST FITS file saved in the last step is just a 2D array of numbers (well, the Data part is; the Header is a quite different animal). And just before you saved it, you saw a representation of that data, in the SkyView page, in the form of a JPG (or GIF?). It looked pretty boring, something like this:

    enter image description here

    What's involved in 'creating contours' is:

    • setting a threshold, so 'noise' can be ignored
    • deciding what the 'spacing' between the contours will be, the ratio of intensity between two adjacent contours: it's commonly SQRT(2) or 2
    • how to 'smooth' the contours: you don't want lots of squarish contours, but equally not all the details smoothed out either

    But before any of that, you have to get the data into a form Python can work with!

    For that, I use the package "astropy.io", and within it "fits" (# is used for comments):

    from astropy.io import fits
    Ffield = fits.open(myFITSfile.fits) # "myFITSfile.fits" is the name of the FITS file I want to open
    Ffieldd = Ffield[0].data # puts the image data into memory, "[0]" indicates the the first "extension"; there's only one
    Ffield.close() # we've got the data, no need to keep the file open
    

    Why use "F"? It stands for "FIRST", and my most often used code has corresponding "N" things ("NVSS")!

    A convenient, and common, choice of 'the noise floor' or 'blank sky' is the median, so:

    import numpy as np
    Ffieldmax = np.amax(Ffieldd) # used for contours
    Fmed = np.median(Ffieldd) # median flux (intensity) value, whole image
    Ffieldsky = Fmed # this is the 'sky' (see RGZ Talk thread)
    Ffieldd0 = Ffieldd - Ffieldsky # sky=0, I once thought it makes things easier, but probably doesn't, so this can be omitted
    

    "RGZ Talk thread" is this RGZ Talk thread

    What about the "0-th contour"? A natural choice is the standard deviation of the noise, but a better one involves the MAD statistic, "median absolute deviation":

    from pylab import * # I don't know if this is necessary, anyone?
    Fmad = np.median(np.abs(Ffieldd-Fmed)) # median absolute deviation
    Ffloor = 1.4826 * Fmad # convert MAD to estimate of sigma, the standard deviation
    

    That same RGZ Talk thread explains what 1.4826 is. This sets the 0-th contour at 1-sigma; to set it at twice sigma, simply add "2" to the equation:

    Ffloor = 2.0 * 1.4826 * Fmad
    

    Suppose we want the contour spacing to be 2, and some generality in the code (so making tweaks will be less tedious and error prone). Here's what I did:

    import scipy.ndimage as ndimage
    Fsf = 2.0 # contour scale factor
    Fsmooth = 4.0 # sigma for gaussian smoothing, about right for 2.4'x2.4' FIRST images
    Flsf = np.log(Fsf) # you'll see why soon
    Ffielddr = (Ffieldmax-Ffieldsky)/Ffloor # dynamic range, above noise
    Fncont = np.int(np.log(Ffielddr)/Flsf) # number of contours (above 0-th), handy for when things go wrong
    Flc0 = np.log(Ffloor)/Flsf # 0-th contour, as a log, base scale factor
    Flmaxcont = Fncont + Flc0 # flux (intensity) of 'highest' contour
    Fcontlev = np.logspace(Flc0,Flmaxcont,num=Fncont+1,base=Fsf) # array, contour levels
    Fhfieldd = ndimage.gaussian_filter(Ffieldd0, sigma=Fsmooth, order=0)
    

    And that's it! Well, things like the contour color, line style, thickness have yet to be set, but those can wait until a later step.

    Lots of routines (?) to understand, some easy (e.g. "median", "log", "abs", "int"), some not (e.g. "logspace", "gaussian_filter"). And some of the math(s) you may not have covered in school/university.

    Next: how to overlay the contours onto the canvass


    Some background etc

    A lot of the background to this post (and more) can be found in the How to decide the 'zero point' for radio contours? thread, here in RGZ Talk. There are quite a few, very helpful, posts by several RGZ Science Team members (although not all are relevant to this post, on how to create contours). As you'll see from that thread, it took me quite a while to figure out how to set the 'floor', create contours, etc!

    In later posts I'll go over how I modified this code to run in 'batch mode', how I handled (some rather too common) 'exceptions', etc. I'll also mention a quite different way to create contours (etc), which I discovered long after I'd developed this code; it's likely easier to use 😦 so I wish I'd learned of it earlier!


    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    I expect there will be plenty of questions this time!

    Posted

  • sisifolibre by sisifolibre

    great work and explanation Jean, I hope to see overlays of other zoonites soon.

    When I have more time I will try to follow your steps 😃

    Posted

  • JeanTate by JeanTate in response to sisifolibre's comment.

    Thanks for the feedback and support, sisifolibre! 😃

    Due to some, um, high-level IRL interrupts, it will be several days before the next installment.

    Posted

  • JeanTate by JeanTate

    Thanks to my fellow zooite who pointed out an omission in my last installment; I added the missing line of code.

    Next: how to overlay the contours onto the canvass

    The core code for this step isn't all that long, but it's far from simple (at least, to me). And the way I coded this step is most definitely not the only way!

    This post also covers some of step 5 (produce the contour overlay image).

    import matplotlib.pyplot as plt
    from matplotlib.pyplot import figure
    ifil = "MySDSSfileS24.jpg" # input SDSS image file
    SDimg = plt.imread(ifil)
    

    matplotlib is incredibly rich, and quite challenging to understand how it works, what bits to use, etc. All the above snippet does is read the SDSS JPG file, which will be used as the canvass, and give it a name.

    Next comes the set-up; establishing the parameters, parameter values, etc to produce an output file that's what you want (or expect)

    mydpi = 100 # "dots per inch"; funny how a digital image connects with the physical world
    xpix = Ffieldh['NAXIS1'] # size of image, in pix
    ypix = Ffieldh['NAXIS2']
    

    NAXIS1 and NAXIS2 are keywords in the FIRST FITS Header. Of course, you could hard code this (e.g. xpix = 600), but since the FIRST data contains these parameters and their values, why not use it? Besides, you don't have to change this every time you use a different-sized image.

    import numpy as np
    w = float(xpix)/mydpi # Ha! without this, figure sizes will be calc using INT division!!!
    h = float(ypix)/mydpi  
    

    The width of the image, and its height. Yes, my 'real' code does contain that comment!

    f = plt.figure(figsize=(w, h), dpi = mydpi, frameon=False)
    ax = plt.Axes(f, [0.,0.,1.,1.])
    ax.set_axis_off()
    f.add_axes(ax)
    ax.imshow(np.flipud(SDimg), origin='lower')
    

    We don't want our image to have axes, but we do need to make sure that the canvass has the same orientation as the contours; the code above does that. And no, I'm not 100% certain I could explain exactly why it works ... but it does.

    Fclinthick = 0.4 # thickness of FIRST contour lines
    Fclincolor = 'r' # color of FIRST contour lines; red in this case
    plt.contour(Fhfieldd, Fcontlev, colors=Fclincolor, linewidths=Fclinthick)
    

    And that's it! 😃

    Except to say that, for plt.contour, the kwarg is "linewidths", with an "s". And "colors", with an "s". Etc. Elsewhere in Python, the kwargs are "linewidth" and "color". {insert roll-eyes smilie here}

    Next: how to produce the image (things not included in this post), and save it.


    Some background etc

    Long after I wrote my code, I learned of an entirely different way to do this step. One that's perhaps a lot more straight-forward. Anyway, I'll talk about this, and other possible ways of overlaying the contours onto the canvass, in a (much) later post.

    In most of my contour overlay images I have text - such as "coords (RA, Dec): ", "scale 2.4'x2.4' (600x600pix)" - a scale bar, and a pair of lines (indicating orientation). In a later post I'll go over the code I used to put such things in the image.


    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    Posted

  • JeanTate by JeanTate

    Next: how to produce the image (things not included in this post), and save it.

    The Python part of this post is really short:

    import matplotlib.pyplot as plt
    from matplotlib.pyplot import figure
    ojpg = "myOverlay.jpg" # name of (output) jpg file with FIRST contours overlaid on an SDSS canvass
    plt.axis('off') # remove the axes, or at least prevent them from showing
    f.savefig(ojpg,dpi=mydpi)
    plt.close(f)
    

    Now you've got a FIRST contour overlay image on your PC/whatever; how to post it here, into RGZ Talk?

    These days I use Imgur: it's easy to use, and to add an image in your account - in an RGZ Talk post, for example - the process is extremely straight-forward:

    • create an Imgur account
    • click on the "upload images" button
    • upload your image
    • once uploaded, open the image (click on its thumbnail)
    • select and copy the "Direct Link (email & IM)" link
    • in the RGZ Talk post in which you want to add your image, click on the "image" icon (to the left of the "code" one)
    • paste the Imgur link
    • click "OK"
    • Done! 😃

    Next: pulling all the code snippets together into a single, working, program

    Also: a post on what I plan to write about afterwards


    Some background, etc

    Thanks to RGZ (and GZ) Science Team member Kyle Willett for the Imgur recommendation (in GZ Talk, here), and to c_cld for drawing my attention to it, and Dolorous_Edd for independently recommending it (see p13 of the Suggestions for RGZ Objects to show with detailed FIRST contours overlaid on SDSS images thread).

    Originally, I used some version or other (possibly more than one, I don't remember any more) of Google's Photo service (or whatever it is/was called), but some time ago they changed their policies, and it became difficult to impossible to use any more. There are a lot, and I mean a LOT, of image (photo) sharing sites and services out there on the internet! Of the ones I tried myself, few are any good for publishing images in RGZ Talk, and some are apparently so poorly designed/monitored that you run a real risk of having your PC infected by adware (or worse, malware).

    It's very easy to save your contour overlay file in a different format; perhaps PNG is the most robust, because this format is 'native' to matplotlib.


    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    Posted

  • antikodon by antikodon

    Great work Jean!

    here test image
    ARG0000u17 SDSS J084712.54+454033.1 enter image description here

    I use Anaconda , which include all package and Spyder as IDE.
    About NAXIS1-2 should I call up header before 'x/ypix = Ffieldh..'?

    Posted

  • JeanTate by JeanTate in response to antikodon's comment.

    Very cool, antikodon! Well done! 😃

    I use Anaconda , which include all package and Spyder as IDE

    When we get to discussing Python installation, etc - perhaps in a separate thread - I hope you'll share your experience.

    About NAXIS1-2 should I call up header before 'x/ypix = Ffieldh..'?

    Yes. The code snippets I've been introducing were not intended to indicate the order in which they appear within a working program. Here's a snippet of the actual program (I'll go over the whole thing in my 'next' post); ignore the fact that some things may have slightly different names than in the snippets above:

    Ffield = fits.open(Fifits)
    Ffieldh = Ffield[0].header 
    Ffieldd = Ffield[0].data
    Ffield.close()
    # extract metadata from header, calculate values
    RA = Ffieldh['CRVAL1'] # I have not discussed this yet
    Dec = Ffieldh['CRVAL2']
    xpix = Ffieldh['NAXIS1'] # size of image, in pix
    ypix = Ffieldh['NAXIS2']
    

    Hope this helps.

    Posted

  • JeanTate by JeanTate in response to JeanTate's comment.

    Also: a post on what I plan to write about afterwards

    Somewhat in order:

    • adding text etc ("annotations")
    • creating FIRST FITS in a semi-automated way
    • ditto, for SDSS canvasses
    • handling exceptions; what to do when things don't work out right
    • adding a scale bar
    • creating an overlay image with both FIRST and NVSS contours
    • creating more than one (2, 20, 200, ...) overlay images in one go; semi-automating the whole process
    • keeping track of everything
    • how to tweak the code to do other overlays (e.g. WISE), use other canvasses (e.g. POSS), etc
    • how to make your own Luptonized images
    • other packages, and options, I've come across
    • more (I'll edit this post to add other items, as I think of them, or my colleagues suggest etc them!)

    As I wrote earlier - "I am particularly unqualified to explain how to install Python, and the basics of how to code in it, etc, so I'm not going to try. Perhaps we could get a discussion going on this topic, in a separate thread?" - I'd be interesting in joining in a discussion on how to set up Python, get the packages you need, etc. However, other than perhaps starting a thread on that topic, I cannot 'lead' any such discussion.

    Posted

  • JeanTate by JeanTate

    Next: pulling all the code snippets together into a single, working, program

    Here it is, with only a couple of changes from what's in the snippets (see below).

    # code for first demo
    # see RGZ Talk thread:
    # I'm thinking of making the 'FIRST contour overlay' Python code I wrote available to any zooite ...
    #
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.pyplot import figure
    import scipy.ndimage as ndimage
    from astropy.io import fits
    from pylab import *
    
    # parameters, and values, fixed for this program
    Fnsignoise = 2.5 # FIRST noise floor, in sigmas
    Fsf = 2.0 # FIRST contour scale factor
    Fsmooth = 4.0 # sigma for gaussian smoothing (FIRST contours)
    
    Fclinthick = 0.4 # thickness of FIRST contour lines
    Fclincolor = 'r' # color of FIRST contour lines
    
    Fifits = "myFITSfile.fits" # input FIRST FITS file
    ifil = "MySDSSfileS24.jpg" # input SDSS image file; the canvass
    ojpg = "myOverlay.jpg" # output contour overlay image file
    
    # open files, obtain headers and data, close file
    Ffield = fits.open(Fifits)
    Ffieldh = Ffield[0].header 
    Ffieldd = Ffield[0].data
    Ffield.close()
    
    # extract metadata from header, calculate values
    xpix = Ffieldh['NAXIS1'] # size of image, in pix
    ypix = Ffieldh['NAXIS2']
    
    # noise floor; 'true sky' is median flux
    Ffieldmax = np.amax(Ffieldd) # used for contours
    Fmed = np.median(Ffieldd) # median flux value, whole image
    Fmad = np.median(np.abs(Ffieldd-Fmed)) # median absolute deviation; see RGZ Talk
    Ffloor = Fnsignoise * 1.4826 * Fmad # convert to est sigma, set min as nsignoise sigma
    Ffieldsky = Fmed # this is the 'sky' (see RGZ Talk)
    Ffieldd0 = Ffieldd - Ffieldsky # sky=0, makes calcs, contours, etc easier (or not)
      
    # contours
    Flsf = np.log(Fsf)
    Ffielddr = (Ffieldmax-Ffieldsky)/Ffloor # dynamic range, above noise floor
    Fncont = np.int(np.log(Ffielddr)/Flsf) # number of contours (above 0-th)
    Flc0 = np.log(Ffloor)/Flsf # 0-th contour, as a log, base sf
    Flmaxcont = Fncont + Flc0
    Fcontlev = np.logspace(Flc0,Flmaxcont,num=Fncont+1,base=Fsf) # array, contour levels
    Fhfieldd = ndimage.gaussian_filter(Ffieldd0, sigma=Fsmooth, order=0) # 'hidden'
    
    # time to plot!
    SDimg = plt.imread(ifil)
    
    # setup
    mydpi = 100
    w = float(xpix)/mydpi # Ha! without this, figure sizes will be calc 
    h = float(ypix)/mydpi # using INT division!!!
    f = plt.figure(figsize=(w, h), dpi = mydpi,frameon=False)
    ax = plt.Axes(f, [0.,0.,1.,1.])
    ax.set_axis_off()
    f.add_axes(ax)
    ax.imshow(np.flipud(SDimg), origin='lower')
              
    # proper contours
    plt.contour(Fhfieldd, Fcontlev, colors=Fclincolor, linewidths=Fclinthick)
           
    # finishing off
    plt.axis('off') # remove the axes, or at least prevent them from showing
    f.savefig(ojpg,dpi=mydpi)
    plt.close(f)
    

    Changes? In an earlier snippet:

    Ffloor = 1.4826 * Fmad # convert MAD to estimate of sigma, the standard deviation
    

    "This sets the 0-th contour at 1-sigma; to set it at twice sigma, simply add "2" to the equation:"

    Ffloor = 2.0 * 1.4826 * Fmad
    

    This is now:

    Fnsignoise = 2.5 # FIRST noise floor, in sigmas. This is in the "parameters" section
    Ffloor = Fnsignoise * 1.4826 * Fmad 
    

    Does the code work?

    Here's myOverlay.jpg, produced with the SDSS canvass I created in an earlier post (centered on SDSS J093139.03+320400.1, (RA, Dec)=(142.9127, 32.0667)), and the FITS file I described earlier too.:

    enter image description here

    😃


    Some background, etc

    I haven't said much about comments yet; in the code you can see I have them both on their own lines and at the end of 'active' lines (recap: a "comment", a piece of code which is not 'interpreted', begins with "#" in Python). There are blank lines too; I could have coded them as an 'empty comment' (a line beginning with "#" and containing nothing else). Blank lines are sometimes required!

    I know Python can work 'interactively', i.e. you can start a program, and have it wait for inputs from your keyboard (for example). This code is not interactive, and in fact I'm not sure I would know how to make it so! 😮

    For example, instead of hard-coding the input files, presumably there's a way to get them interactively, from what you enter on your keyboard once the program is running. Anyone know how to do this?


    Next: adding text etc ("annotations")

    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    Posted

  • JeanTate by JeanTate in response to JeanTate's comment.

    Next: adding text etc ("annotations")

    A contour overlay image is certainly nice, but to have scientific merit it needs to some context. Perhaps the most important context is:

    • position: where, on the sky, it is; this is best done with some reference or other to (RA, Dec)
    • image scale, or size; there are quite a few ways to do this, perhaps the most common is a scale bar
    • orientation: kinda like a compass, pointing which direction is N, and which E (say)

    I'll give the code I used for all three; however, I'll illustrate the second by a method other than a scale bar. I've also indicated in which section of the code, in my post yesterday, these code snippets go.

    # where does this code go? In the "extract metadata from header, calculate values" section
    RA = Ffieldh['CRVAL1']
    Dec = Ffieldh['CRVAL2']
    xscal = abs(Ffieldh['CDELT1']) # scale, in deg/pix; abs 'cause of orientation
    yscal = abs(Ffieldh['CDELT2'])
    xsiz = xscal * xpix # CDELT1 is -ve because E faces left?
    ysiz = yscal * ypix
    
    # where does this code go? After the "proper contours" section
    # coords, as annotation
    RAs = '%8.4f' % (RA)
    Decs = '%8.4f' % (Dec)
    coords = 'coords (RA, Dec): '+RAs+' '+Decs+' (center)'
    ax.annotate(coords, xy=(0.57,0.97),xycoords='axes fraction', color='w', fontsize=8)
    
    # scale, as text, not scale bar; this is in the same section as "coords"
    scal = 'scale: '+str(xsiz*60)+"' x "+str(ysiz*60)+"' ("+str(xpix)+'x'+str(ypix)+'pix)'
    ax.annotate(scal, xy=(0.02,0.97), xycoords='axes fraction', color='w', fontsize=8)
    
    # orientation, just a pair of lines, with "E" and "N"; this is in the same section as "coords"
    plt.vlines(60,10,50,colors='y')
    plt.hlines(10,20,60,colors='y')
    ax.annotate('E', xy=(12,8), xycoords='axes pixels', color='y', fontsize=8)
    ax.annotate('N', xy=(60,58), xycoords='axes pixels', color='y', fontsize=8)
    

    Two comments:

    • there are (at least) two different ways to specify where, on your image, you want the text to go, "axes fraction" or "axes pixels"; you can see I've used both
    • concatenating bits of text and parameter values, to produce the output text string you want, can be tricky; lots of trial and error!

    What's the result? Here:

    enter image description here


    Some background, etc

    It is not necessary to provide any of position, scale, and orientation on the image itself. For example, you can add axes to provide position information, and/or the scale and orientation information in the accompanying text (in you paper, or post, or whatever). However, putting it on the image reduces the chances for error.

    You can see that I used two new sets of parameters, obtained from the FITS Header: ('CRVAL1', CRVAL2') and ('CDELT1', 'CDELT2'). This is one of the great thing about FITS: the key information about the data (image) is in the Header! 😃 All the more reason to be extremely careful that you keep good track of the key metadata for the 'dumb' SDSS JPG canvass! 😮

    I left the scale bar to a later post. Why? Because my implementation is extremely kludgy, clumsy, etc. So I would definitely not recommend it to anyone! However, it works, and perhaps illustrates what you sometimes need to do if your sole goal is "it works".


    Next: creating FIRST FITS in a semi-automated way

    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    Posted

  • JeanTate by JeanTate

    Next: creating FIRST FITS in a semi-automated way

    SkyView has a nice JAR tool that I have used. It's available from this site, and the documentation is easy to follow, and short.

    Here's code for obtaining the FIRST FITS file I've been using, as an output file called myFITSfile.fits:

    java -jar skyview.jar survey=FIRST position='142.9127, 32.0667' size=0.04 Pixels=600 output=myFITSfile
    

    Looks pretty simple, right? All the parameters are obvious, and the values easy to understand, right? Just remember to enclose the position in (single) quote marks. And put spaces where you need to, and none where you don't.

    To run this, I open a "Command Prompt" window (you can find it under "All Programs", from the "Start" button, for machines running Windows), then navigate to the directory where I want the output file to go. You could then simply type the command (line of code above), but that's highly error prone, and useless if you want, say, 20 FITS files.

    So, I open a Notepad file, type the code into it, then save it as a "Windows Batch File", e.g. MyJAR.bat (the ".bat" extension means the Windows OS will treat it as a Windows Batch File). In the same directory as where I want the output file to go.

    Then, at the Command Prompt, type "MyJar.bat" (actually, I think simply "MyJAR" will work), without the quotes, of course.

    A few seconds later, Done! 😃

    And to create two FIRST FITS files from the one .bat file? Easy! Here's such a .bat file, with the second FIRST position (194.1800, 11.3000), and saved as myFITSfile2.fits:

     java -jar skyview.jar survey=FIRST position='142.9127, 32.0667' size=0.04 Pixels=600 output=myFITSfile
     java -jar skyview.jar survey=FIRST position='194.1800, 11.3000' size=0.04 Pixels=600 output=myFITSfile2
    

    It can become tedious typing so many, very similar, lines. Error-prone too. So why not semi-automate the process?

    I use a spreadsheet, which I also use to keep track of my work ("keeping track of everything" is the topic of a later post). In that I have columns for:

    • RA
    • Dec
    • image size, in degrees
    • an object name

    I could also have a column for the number of pixels, but since I almost always use 600, I hard code that.

    To create a JAR line, I use Open Office's CONCATENATE function; for example (I have to split this across two lines):

    =CONCATENATE("java -jar skyview.jar survey=FIRST position='";G170;",";H170;"' size=";O170;
    "Pixels=600 output=";C170;"F";INT(600*O170))
    

    Cells G170 and H170 contain the (RA, Dec) coordinates of my target, whose name is in cell C170. Cell O170 is 0.04, for an image of size 2.4'x2.4'. If the target name is "MyFITS", the JAR output filename will be "MyFITSF24.fits"

    To create a .bat file, simple copy/paste into a Notepad file, and save as MyJAR.bat

    Done! 😃


    Some background, etc

    If you don't know how to navigate to a directory, in the Command Prompt window, holler.

    There's something sneaky about the RA and Dec coordinates that are used as inputs to the SkyView JAR tool: they must contain no more than four significant digits! 😮 Why? I have no idea, just that if you have five or more the application will choke (not always, strangely, but often enough).

    What if there's no FIRST data for the field around (RA, Dec) you entered? I'll cover that in a later post; suffice it to say that the app won't alert you to that fact.

    Make sure your output FITS filenames are unique, a different one for each target; if not, you'll lose some work, as there's no warning about overwriting/replacing.


    Next: creating SDSS canvasses in a semi-automated way

    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    Posted

  • JeanTate by JeanTate

    Next: creating SDSS canvasses in a semi-automated way

    The Python code I use for this task introduces three new things, cf what's in my earlier Python code posts:

    • csv: a tool to read a CSV (comma-separated values) file
    • urllib: a tool to retrieve the contents of a URL (among other things)
    • a loop, implemented by "try" ... "finally"

    You might guess that this Python code takes as input a CSV file, and outputs a JPG one; you'd be right! 😃

    The input CSV file is simplicity itself: just an ID (name of the target), and (RA, Dec). As described in an earlier post, I open a blank Notepad file, and copy/paste the cells in my spreadsheet, containing the relevant data. Here's the line in the spreadsheet which does this:

    =CONCATENATE(C170;", ";ROUND(G170;5);", ";ROUND(H170;5))
    

    Cell C170 contains the target's name ("myTarget", say), and (G170, H170) its (RA, Dec).

    Save the Notepad file, as a CSV; e.g. "SDSS01.csv"

    Here's the Python code:

    # retrieve SDSS .jpg image
    # using single input file, in csv format
    # scale and image size fixed in tailstr
    # file suffix in the urllib line
    
    import csv
    import urllib
    import numpy as np
    
    headstr = "http://skyservice.pha.jhu.edu/DR12/ImgCutout/getjpeg.aspx?ra="
    tailstr = "&scale=0.24&width=600&height=600"
    
    f = open('SDSS01.csv')
    
    try:
        reader = csv.reader(f)
        for row in reader:
            ID = row[0]
            ra = row[1]
            dec = row[2]
            url = headstr+ra+"&dec="+dec+tailstr
            urllib.urlretrieve(url, ID+"S24.jpg")
            print ID, ra, dec
    
    finally:
        f.close()
    

    If there are three targets, say, with names "myTarget01", "myTarget02", and "myTarget03", this code will produce three JPG files, 600x600 pix in size, with a scale of 0.24"/pix; they will have these names: "myTarget01S24.jpg", "myTarget02S24.jpg", and "myTarget03S24.jpg" (respectively).

    Done! 😄


    Some background, etc

    There are several ways to implement a loop in Python; I don't think "try" ... "finally" has any particular advantages, or disadvantages, over other ways.

    It's easy to generalize this code so that the image size, in pixels, and scale are inputs, rather than hard-coded. And to change the last three characters of the output filename accordingly.

    I use "S24" (for an image derived from SDSS with a scale of 0.24"/pix) as an easy way to keep track, and reduce errors. For example, in creating a FITS contour overlay image 2.4'x2.4' and 600x600 pix in size, if the FITS input filename ends with "F24" and the SDSS canvass "S24", I can exploit that in the Python code, so that if either is missing, or misnamed, it'll create an exception. Of course, this requires that the two input files were created correctly ...

    I'll discuss this further in my "handling exceptions; what to do when things don't work out right" and "keeping track of everything" future posts.


    Next: adding a scale bar

    And it goes without saying: comments, questions, suggestions, ... feedback is always welcome. 😄

    Posted

  • JeanTate by JeanTate

    I'm going to take a break from writing these, so the next one in the series will be at least a week from now.

    In the meantime, please let me know if you've been able to follow along, what the results of trying to create your own contour overlay images have been, etc, etc, etc.

    As I often say, happy hunting! 😃

    Posted

  • ChrisMolloy by ChrisMolloy

    Just a quick update on producing single First contour overlay images with Python code. Here's my first attempt below, ARG0001ly0, SDSS J093139.03+320400.1.

    enter image description here

    I'm using Python 3.5.1 so I had to make slight amendments to the code Jean provided, which I believe is Python 2.7. Also I'm using a mac osx system, which has its own quirks and variables which are different from a windows operating system.

    Firstly when running the python code I kept getting an error stating no " fitsfile found", so I had to specify the location directly, and make this applicable for the other files. The modifications I made to the code are as follows:

    Fifits = "/Users/chris/Desktop/RGZPythonimages/fitsfile.fits" # input FIRST FITS file
    ifil = "/Users/chris/Desktop/RGZPythonimages/test.jpeg" # input SDSS image file; the canvass
    ojpg = "/Users/chris/Desktop/RGZPythonimages/myOverlay.jpeg" # output contour overlay image file
    

    Once I got past this error, I then had to make another alteration, as I was getting an error relating to the Finishing off section. Here's the code that worked for me, again, specifying the file location directly, but also changing the file from a jpeg to a png.

    f.savefig("/Users/chris/Desktop/RGZPythonimages/myOverlay.png")
    

    Hope this post helps for anyone thinking of producing First contour overlay images with Python code. The process outlined above by Jean is relatively straight forward and easy to follow.

    Posted

  • JeanTate by JeanTate in response to ChrisMolloy's comment.

    @ChrisMolloy Congrats, well done! 😃

    Posted

  • JeanTate by JeanTate

    A small update: how to deal with NaN?

    What is "NaN"? It stands for "not a number", and you'll come across it in FITS files of fields which include regions that are not "imaged". If you look at FIRST images, via the cutout service or SkyView, you'll see these are blank/black parts of the image.

    In the code I've posted, if there's even one NaN array element in the data block, no contours will be created.

    Here's one way round this: replace NaN with 0, using isnan. For example, in a data block called Ffieldd, this line of code will replace all instances of NaN with 0.0:

    Ffieldd[np.isnan(Ffieldd)] = 0.0
    

    As you might guess, isnan is part of numpy (np in the code above), so you need to make sure you have numpy imported ... which you almost certainly will have, as there are plenty of other lines of code which use it! 😃

    Posted

  • JeanTate by JeanTate

    There's a new, cool, section in Zooniverse Talk, Data processing.

    I wrote a few words about this RGZ Talk thread, in a Discussion in that section, which I titled One volunteer's experience with Python (perhaps not directly relevant?).

    One fellow zooite/volunteer suggested that I put the code up on (in?) GitHub, a site which I knew very little about. I'm a bit reluctant, other than perhaps creating an account of my own and putting the code there.

    What do you think?

    Posted