I've noticed that GSObject.draw can be called without the image keyword. The documentation says that an automatically sized image is returned if image = None. How exactly is this automatic size determined? Does it make sure to retain some fraction of the galaxy's flux within the image, and if so what fraction?
Yes. The summary is that it tries to get at least 99.5% of the flux.
Each atomic object has an internal attribute, stepk, which is the grid size to use in the fourier image. This is based on a real-space image including at least 99.5% of the flux to avoid aliasing in the FFT.
Then as you combine objects (with Add, Convolve, etc.) or transform them (shear, dliate, etc.) the value of stepk is kept up to date for the new profile, sometimes with heuristics when the precise formula would be too unwieldy, but we usually tried to err on the side of being conservative here.
Then when you draw the final object with image=None, it uses the final stepk value to calculate a stamp size with N = 2pi * wmult / (scale * stepk), where scale is the pixel scale of the image, which you would normally specify explicitly (im = obj.draw(scale=pixel_scale)), and wmult is an optional parameter whose job is pretty much just this -- to make larger images than it would do by default.
Then this N is rounded up to either 2^k or 3 * 2^k. This is to help make subsequent FFTs more efficient.
If you want to change the 99.5% value, this is also possible with the GSParams class. It is the parameter alias_threshold. Well, actually 1-alias_threshold. The default alias_threshold is 5e-3, but you can lower that if you want by doing something like the following:
gsparams = galsim.GSParams(alias_threshold=1.e-3)
gal = galsim.Sersic(..., gsparams=gsparams)
psf = galsim.Kolmogorov(..., gsparams=gsparams)
pix = galsim.Pixel(..., gsparams=gsparams)
conv = galsim.Convolve([gal,psf,pix])
im = conv.draw(scale=pixel_scale)