Skip to content

Starting to add more (con)text #220

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Starting to add more (con)text #220

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
143 changes: 117 additions & 26 deletions Measurement/UndersampledMoments.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,16 +8,16 @@
"source": [
"# Biased moments of undersampled sources\n",
"<br>Owner(s): **Andrew Bradshaw** ([@andrewkbradshaw](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@andrewkbradshaw))\n",
"<br>Last Verified to Run: **2019-08-09**\n",
"<br>Verified Stack Release: **18.0**\n",
"<br>Last Verified to Run: **2019-08-23**\n",
"<br>Verified Stack Release: **18.1.0**\n",
"\n",
"This notebook demonstrates how the measurement of centroid, size, and orientation angles of small sources (relative to the pixel grid) are biased due to undersampling and/or assumptions in the method of analysis. A possible improvement, in the form of a modified weight function, is shown to have a positive effect.\n",
"This notebook performs measurements of the centroid, size, and orientation angles of small sources. It goes on to show how these measurements are biased (relative to the pixel grid) due to undersampling and/or assumptions in the method of analysis. A possible improvement, in the form of a modified weight function, is shown to have a positive effect.\n",
"\n",
"### Learning Objectives:\n",
"\n",
"After working through this tutorial you should be able to: \n",
"1. Simulate stars and galaxies with a variety of properties, including sub-pixel centroid, size, and orientation angle.\n",
"2. Measure the centroid, size, and orientation angle using a variety of algorithms.\n",
"2. Measure the centroid, size, and orientation angle using a variety of algorithms (including the Stack default).\n",
"\n",
"### References:\n",
"* [galsim.Shear](http://galsim-developers.github.io/GalSim/_build/html/shear.html) - Documentation on how GalSim defines shear\n",
Expand Down Expand Up @@ -65,7 +65,7 @@
"\n",
"We start by setting up a class that uses [GalSim](http://galsim-developers.github.io/GalSim/_build/html/index.html) to generate simulated images and sources. GalSim can provide various levels of realism (e.g., readout noise, sky background, silicon sensor effects, etc.), and these are configurable during the creation of an instance. We also implement a function to use GalSim to generate a Gaussian source with a given centroid, ellipticity, and orientation.\n",
"\n",
"We break up the creation of the `undersampled_galsim` class so that the various components can be explored in turn, rather than as one monolithic block. To do this, we subclass the existing class and add new methods."
"We break up the creation of the `UndersampledGalSim` class so that the various components can be explored in turn, rather than as one monolithic block. To do this, we subclass the existing class and add new methods."
]
},
{
Expand Down Expand Up @@ -93,7 +93,26 @@
" diffusion_factor=0,stamp_width=31,pixel_scale=0.2,gain=1,\n",
" rngseed=8675309,do_hsm_moments_truefalse=False,sim_sensor_truefalse=False,\n",
" default_gsparams=True):\n",
" \"\"\"Initialize the imager with some basic settings\"\"\"\n",
" \"\"\"Initialize the imager with some basic settings\n",
" \n",
" Parameters\n",
" ----------\n",
" method_of_imaging: how to generate the image ('phot','fft','no_pixel')\n",
" read_noise: Gaussian std dev. of CCD readout noise\n",
" sky_level_pixel: the background noise level (subtracted off)\n",
" diffusion_factor: silicon sensor diffusion factor\n",
" stamp_width: width of each postage stamp (default: 31)\n",
" pixel_scale: pixel scale in arcsec/pixel (LSST = 0.2)\n",
" gain: CCD gain (default: 1)\n",
" rngseed: random number generator seed\n",
" do_hsm_moments_truefalse: do hsm moments measurement\n",
" sim_sensor_truefalse: whether or not to simulate the sensor & diffusion\n",
" default_gsparams: use default GalSim parameters\n",
" \n",
" Returns\n",
" -------\n",
" None\n",
" \"\"\"\n",
" self.method = method_of_imaging # default: phot, fft; no_pixel' might take forever?\n",
" self.sensor = sim_sensor_truefalse # whether or not to simulate the sensor & diffusion\n",
" self.pixel_scale=pixel_scale # pixel scale in arcsec/pixel, LSST = 0.2\n",
Expand Down Expand Up @@ -123,7 +142,20 @@
" def draw_gaussian(self,offset=(0.0,0.0),gal_sigma_pix=1.5,gal_ellip=0,\n",
" gal_angle_deg=0,gal_flux=1e5):\n",
" \"\"\"Draw a Gauassian on the image, using the provided centroid offsets, widths, ellipticity\n",
" angle, and flux. Returns the galsim object and HSM moments (if do_hsm is set to True)\"\"\"\n",
" angle, and flux. Returns the galsim object and HSM moments (if do_hsm is set to True)\n",
" \n",
" Parameters\n",
" ----------\n",
" offset : coordinate offset (preserves sub-pixel offset in even-sized images)\n",
" gal_sigma_pix : standard deviation of Gaussian (unit: pixels)\n",
" gal_ellip : source ellipticity, e; see GalSim.Shear\n",
" gal_angle_deg : source position angle, beta; see GalSim.shear (unit: deg)\n",
" gal_flux : source flux\n",
" \n",
" Returns\n",
" -------\n",
" final,hsm : image array, hsm fit parameters\n",
" \"\"\"\n",
" gal = galsim.Gaussian(sigma = gal_sigma_pix*self.pixel_scale, gsparams=self.gsparams)\n",
" gal_shape = galsim.Shear(e=gal_ellip, beta=gal_angle_deg * galsim.degrees)\n",
" gal = gal.shear(gal_shape)\n",
Expand All @@ -150,7 +182,9 @@
"toc-hr-collapsed": true
},
"source": [
"## Functions to compute first and second moments on images"
"## Custom calculations of first and second moments on images\n",
"\n",
"We define two custom methods to calculate first and second moments on images. These are non-adaptive methods that radially weights by the standard deviation of the input Gaussian multiplied by some chosen pre-factor.\n"
]
},
{
Expand All @@ -162,7 +196,20 @@
"def calc_firstmoms(img,coord,winsize,thresh,gwin_sig):\n",
" #print(img[:5,:5],coord,winsize,thresh,gwin_sig)\n",
" \"\"\"Calculates the x and y first moments of the image windowed by winsize at\n",
" coord (x,y), weighted with a Gaussian of width gwin_sig\"\"\"\n",
" coord (x,y), weighted with a Gaussian of width gwin_sig\n",
" \n",
" Parameters\n",
" ----------\n",
" img : input image array\n",
" coord : centroid coordinate (pixels)\n",
" winsize : window size (pixels)\n",
" thresh : threshold value for defining pixels in the source\n",
" gwin_sig : width of the Gaussian weighting function (pixels)\n",
" \n",
" Returns\n",
" -------\n",
" xmom, ymom : the x and y first moments of the source\n",
" \"\"\"\n",
" xc,yc=coord\n",
" xmom,ymom=0.,0.\n",
" thesum=0.\n",
Expand All @@ -186,7 +233,20 @@
"def calc_secmoms_gwin(img,coord,winsize,thresh,gwin_sig):\n",
" \"\"\"Calculates the x, y, xy, and radial second moments of the image windowed by\n",
" winsize at coord (x,y) and thresholded at thresh, weighted by a Gaussian of \n",
" width gwin_sig\"\"\"\n",
" width gwin_sig\n",
" \n",
" Parameters\n",
" ----------\n",
" img : input image array\n",
" coord : centroid coordinate (pixels)\n",
" winsize : window size (pixels)\n",
" thresh : threshold value for defining pixels in the source\n",
" gwin_sig : width of the Gaussian weighting function (pixels)\n",
"\n",
" Returns\n",
" -------\n",
" rmom, xmom, ymom, xymom : radial, x, y, and xy second moments of the source\n",
" \"\"\"\n",
" xc,yc=coord\n",
" r2mom,x2mom,y2mom,xy2mom=0.,0.,0.,0.\n",
" thesum=0.\n",
Expand Down Expand Up @@ -215,6 +275,13 @@
" return rmom, xmom, ymom, xymom\n",
"\n",
"def calc_gauss_diff(i,j,xcen,ycen,sigx,sigy):\n",
" \"\"\" \n",
" Parameters\n",
" ----------\n",
" \n",
" Returns\n",
" -------\n",
" \"\"\"\n",
" ydiff=1./24.*(((i-xcen)**2/sigx**2-1.)/sigx**2 + ((j-ycen)**2/sigy**2-1.)/sigy**2)\n",
" ygauss=np.exp(-.5*((i-xcen)**2/sigx**2))*np.exp(-.5*((j-ycen)**2/sigy**2))\n",
" return ydiff*ygauss\n",
Expand Down Expand Up @@ -272,6 +339,17 @@
"outputs": [],
"source": [
"def calc_and_print_moments(stamp,hsm):\n",
" \"\"\"Calculate and pring the first and second moments for a postage stamp with predetermined hsm measurements\n",
" \n",
" Parameters\n",
" ----------\n",
" stamp : postage stamp array\n",
" hsm : moments calculated from hsm\n",
" \n",
" Returns\n",
" -------\n",
" None\n",
" \"\"\"\n",
" x_hsm,y_hsm,sig_hsm=hsm.moments_centroid.x,hsm.moments_centroid.y,hsm.moments_sigma\n",
" # calculate my moments\n",
" thresh=0\n",
Expand Down Expand Up @@ -330,7 +408,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Making mosaics of objects"
"## Making mosaics of objects\n",
"\n",
"Now that we have the ability to generate and measure one object, we would like to expand our analysis to a statistical set of objects. We add functionality to make a mosaic of sources on a single image."
]
},
{
Expand Down Expand Up @@ -435,7 +515,9 @@
"toc-hr-collapsed": true
},
"source": [
"## Making a mosaic of small spots with a range of sub-pixel centroids"
"## Making a mosaic of small spots with a range of sub-pixel centroids\n",
"\n",
"We create a mosaic of small spots distributed over a range of centroids. The centroids are given sub-pixel offsets to determine the biases of our various fitting methods."
]
},
{
Expand All @@ -456,7 +538,7 @@
" angle_minmax=(0,0),flux_minmax=(1e5,1e5))\n",
"\n",
"mosaic,cat_in=ugal.make_mosaic(cat_in)\n",
"print(time.time()-tstart,\" to simulate \",str(len(cat_in)),\" galaxies\")\n",
"print(\"It took %.2fs to simulate %s spots\"%(time.time()-tstart,str(len(cat_in))))\n",
"\n",
"plt.figure(figsize=(16,16))\n",
"plt.imshow(mosaic,origin='lower',vmax=100),plt.colorbar()"
Expand Down Expand Up @@ -567,7 +649,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"### Merge in stack measurements to catalog"
"### Merge in stack measurements to catalog\n",
"\n",
"Now that we've run the Stack measurements, we would like to merge them into the same catalog that contains the hsm measurements performed when generating the simulated sources."
]
},
{
Expand Down Expand Up @@ -628,7 +712,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Measurement using sextractor wrapper"
"## Measurement using SExtractor wrapper\n",
"\n",
"Below we include some code that can be used to install `sep`, a python wrapper around Source Extractor, and calculate the SExtractor moments."
]
},
{
Expand Down Expand Up @@ -668,7 +754,9 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Measure moments using a non-adaptive method, radially weighting by the input Gaussian $\\sigma$ multiplied by a factor"
"## Custom moments measurement \n",
"\n",
"We now perform a first moments calculation using the method that we defined above in `calc_firstmoms`. This is a non-adaptive method that radially weights by the standard deviation of the input Gaussian multiplied by some chosen pre-factor (we choose this factor to be 2 and 5 below)."
]
},
{
Expand Down Expand Up @@ -701,7 +789,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plotting the centroid bias"
"### Plot centroid bias"
]
},
{
Expand Down Expand Up @@ -741,7 +829,10 @@
},
"source": [
"# Size bias\n",
"Making a mosaic of spots with a range of input sizes $\\sigma$, takes about a minute for 100 objects with 5 million photons"
"\n",
"Make a mosaic of spots with a range of input sizes, $\\sigma$. When generating the mosaic, we also perform a simple moments calculation, which we store in the `cat_in` variable. \n",
"\n",
"This cell takes about **1 minute** to run for 100 objects with 5 million photons."
]
},
{
Expand Down Expand Up @@ -773,6 +864,13 @@
"fits.writeto(mo_dir+'mosaic_sizes.fits',mosaic,overwrite=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we run the stack on the spot mosaic. We add the Stack moment measurements (i.e., the quantities starting with `base_SdssShape`) to the catalog containing the simple moments measurements that we performed when generating the mosaic."
]
},
{
"cell_type": "code",
"execution_count": null,
Expand Down Expand Up @@ -880,7 +978,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"## plotting size bias"
"### Plot size bias"
]
},
{
Expand Down Expand Up @@ -923,13 +1021,6 @@
"plt.tight_layout(),plt.subplots_adjust(wspace=.4)\n",
"plt.savefig(mo_dir+'/biased_moments-size.png',dpi=150)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down