Skip to content
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

Jump to bottom

DM-20566: Replace afwGeom with geom and fix some resource warnings #62

Merged
merged 2 commits into from
Jul 24, 2019
Merged

DM-20566: Replace afwGeom with geom and fix some resource warnings #62

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e6476af
Replace afwGeom with geom where appropriate
timj Jul 23, 2019
a0da2a2
Fix some ResourceWarnings
timj Jul 23, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
78 changes: 39 additions & 39 deletions python/lsst/obs/sdss/convertasTrans.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,9 +25,9 @@
import numpy as np

import lsst.afw.image as afwImage
import lsst.afw.geom as afwGeom
from lsst.afw.geom import makeSkyWcs
import lsst.afw.table as afwTable
import lsst.geom as geom
import lsst.meas.astrom.sip as sip

deg2rad = np.pi / 180.
Expand Down Expand Up @@ -134,8 +134,8 @@ def createWcs(x, y, mapper, order=4, cOffset=1.0):
src.set(centroidKey.getY(), y[i])

cat = catTable.makeRecord()
cat.set(catTable.getCoordKey().getRa(), afwGeom.Angle(ra_rad[i], afwGeom.radians))
cat.set(catTable.getCoordKey().getDec(), afwGeom.Angle(dec_rad[i], afwGeom.radians))
cat.set(catTable.getCoordKey().getRa(), geom.Angle(ra_rad[i], geom.radians))
cat.set(catTable.getCoordKey().getDec(), geom.Angle(dec_rad[i], geom.radians))

mat = afwTable.ReferenceMatch(cat, src, 0.0)
matches.append(mat)
Expand All @@ -144,11 +144,11 @@ def createWcs(x, y, mapper, order=4, cOffset=1.0):

# CRPIX1 = Column Pixel Coordinate of Ref. Pixel
# CRPIX2 = Row Pixel Coordinate of Ref. Pixel
crpix = afwGeom.Point2D(x[0] + cOffset, y[0] + cOffset)
crpix = geom.Point2D(x[0] + cOffset, y[0] + cOffset)

# CRVAL1 = RA at Reference Pixel
# CRVAL2 = DEC at Reference Pixel
crval = afwGeom.SpherePoint(ra_rad[0], dec_rad[0], afwGeom.radians)
crval = geom.SpherePoint(ra_rad[0], dec_rad[0], geom.radians)

# CD1_1 = RA degrees per column pixel
# CD1_2 = RA degrees per row pixel
Expand All @@ -158,9 +158,9 @@ def createWcs(x, y, mapper, order=4, cOffset=1.0):
ULl = mapper.xyToRaDec(0., 1.)
LRl = mapper.xyToRaDec(1., 0.)

LLc = afwGeom.SpherePoint(LLl[0], LLl[1], afwGeom.radians)
ULc = afwGeom.SpherePoint(ULl[0], ULl[1], afwGeom.radians)
LRc = afwGeom.SpherePoint(LRl[0], LRl[1], afwGeom.radians)
LLc = geom.SpherePoint(LLl[0], LLl[1], geom.radians)
ULc = geom.SpherePoint(ULl[0], ULl[1], geom.radians)
LRc = geom.SpherePoint(LRl[0], LRl[1], geom.radians)

cdN_1 = LLc.getTangentPlaneOffset(LRc)
cdN_2 = LLc.getTangentPlaneOffset(ULc)
Expand All @@ -180,7 +180,7 @@ def validate(xs, ys, mapper, wcs):
dists = []
for i in range(len(xs)):
tuple1 = mapper.xyToRaDec(xs[i], ys[i])
coord1 = afwGeom.SpherePoint(tuple1[0], tuple1[1], afwGeom.radians)
coord1 = geom.SpherePoint(tuple1[0], tuple1[1], geom.radians)
coord2 = wcs.pixelToSky(xs[i], ys[i])
dist = coord1.separation(coord2).asArcseconds()
dists.append(dist)
Expand All @@ -189,36 +189,36 @@ def validate(xs, ys, mapper, wcs):


def convertasTrans(infile, filt, camcol, field, stepSize=50, doValidate=False):
hdulist = fits.open(infile)
t0 = hdulist[0].header['ccdarray']
if t0 != 'photo':
raise RuntimeError('*** Cannot support ccdarray: %s' % (t0,))

camcols = hdulist[0].header['camcols']
filters = hdulist[0].header['filters']
node_deg = hdulist[0].header['node']
incl_deg = hdulist[0].header['incl']
node_rad = node_deg * deg2rad
incl_rad = incl_deg * deg2rad

cList = [int(cc) for cc in camcols.split()]
fList = filters.split()

try:
cIdx = cList.index(camcol)
except Exception:
print("Cannot extract data for camcol %s" % (camcol))
return None

try:
fIdx = fList.index(filt)
except Exception:
print("Cannot extract data for filter %s" % (filt))
return None

ext = cIdx * len(fList) + (fIdx + 1)
ehdr = hdulist[ext].header
edat = hdulist[ext].data
with fits.open(infile) as hdulist:
t0 = hdulist[0].header['ccdarray']
if t0 != 'photo':
raise RuntimeError('*** Cannot support ccdarray: %s' % (t0,))

camcols = hdulist[0].header['camcols']
filters = hdulist[0].header['filters']
node_deg = hdulist[0].header['node']
incl_deg = hdulist[0].header['incl']
node_rad = node_deg * deg2rad
incl_rad = incl_deg * deg2rad

cList = [int(cc) for cc in camcols.split()]
fList = filters.split()

try:
cIdx = cList.index(camcol)
except Exception:
print("Cannot extract data for camcol %s" % (camcol))
return None

try:
fIdx = fList.index(filt)
except Exception:
print("Cannot extract data for filter %s" % (filt))
return None

ext = cIdx * len(fList) + (fIdx + 1)
ehdr = hdulist[ext].header
edat = hdulist[ext].data

if ehdr['CAMCOL'] != camcol or ehdr['FILTER'] != filt:
print("Extracted incorrect header; fix me")
Expand Down
100 changes: 50 additions & 50 deletions python/lsst/obs/sdss/convertfpM.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
from astropy.io import fits

import lsst.afw.image as afwImage
import lsst.afw.geom as afwGeom
import lsst.geom as geom


class Span(object):
Expand Down Expand Up @@ -102,55 +102,55 @@ def setMask(self, mask):


def convertfpM(infile, allPlanes=False):
hdr = fits.open(infile)
hdr[0].header['RUN']
hdr[0].header['CAMCOL']
hdr[0].header['FIELD']
nRows = hdr[0].header['MASKROWS']
nCols = hdr[0].header['MASKCOLS']
hdr[0].header['NPLANE']

names = hdr[-1].data.names
if ("attributeName" not in names) or ("Value" not in names):
raise LookupError("Missing data in fpM header")

planes = hdr[-1].data.field("attributeName").tolist()
mask = afwImage.Mask(afwGeom.ExtentI(nCols, nRows))

# Minimal sets of mask planes needed for LSST
interpPlane = planes.index("S_MASK_INTERP") + 1
satPlane = planes.index("S_MASK_SATUR") + 1
crPlane = planes.index("S_MASK_CR") + 1

interpBitMask = afwImage.Mask.getPlaneBitMask("INTRP")
satBitMask = afwImage.Mask.getPlaneBitMask("SAT")
crBitMask = afwImage.Mask.getPlaneBitMask("CR")

listToSet = [(interpPlane, interpBitMask),
(satPlane, satBitMask),
(crPlane, crBitMask)]

# Add the rest of the SDSS planes
if allPlanes:
for plane in ['S_MASK_NOTCHECKED', 'S_MASK_OBJECT', 'S_MASK_BRIGHTOBJECT',
'S_MASK_BINOBJECT', 'S_MASK_CATOBJECT', 'S_MASK_SUBTRACTED', 'S_MASK_GHOST']:
idx = planes.index(plane) + 1
planeName = re.sub("S_MASK_", "", plane)
mask.addMaskPlane(planeName)
planeBitMask = afwImage.Mask.getPlaneBitMask(planeName)
listToSet.append((idx, planeBitMask))

for plane, bitmask in listToSet:
if len(hdr) < plane:
continue

if hdr[plane].data is None:
continue

nmask = len(hdr[plane].data)
for i in range(nmask):
frow = hdr[plane].data[i]
Objmask(frow, bitmask).setMask(mask)
with fits.open(infile) as hdr:
hdr[0].header['RUN']
hdr[0].header['CAMCOL']
hdr[0].header['FIELD']
nRows = hdr[0].header['MASKROWS']
nCols = hdr[0].header['MASKCOLS']
hdr[0].header['NPLANE']

names = hdr[-1].data.names
if ("attributeName" not in names) or ("Value" not in names):
raise LookupError("Missing data in fpM header")

planes = hdr[-1].data.field("attributeName").tolist()
mask = afwImage.Mask(geom.ExtentI(nCols, nRows))

# Minimal sets of mask planes needed for LSST
interpPlane = planes.index("S_MASK_INTERP") + 1
satPlane = planes.index("S_MASK_SATUR") + 1
crPlane = planes.index("S_MASK_CR") + 1

interpBitMask = afwImage.Mask.getPlaneBitMask("INTRP")
satBitMask = afwImage.Mask.getPlaneBitMask("SAT")
crBitMask = afwImage.Mask.getPlaneBitMask("CR")

listToSet = [(interpPlane, interpBitMask),
(satPlane, satBitMask),
(crPlane, crBitMask)]

# Add the rest of the SDSS planes
if allPlanes:
for plane in ['S_MASK_NOTCHECKED', 'S_MASK_OBJECT', 'S_MASK_BRIGHTOBJECT',
'S_MASK_BINOBJECT', 'S_MASK_CATOBJECT', 'S_MASK_SUBTRACTED', 'S_MASK_GHOST']:
idx = planes.index(plane) + 1
planeName = re.sub("S_MASK_", "", plane)
mask.addMaskPlane(planeName)
planeBitMask = afwImage.Mask.getPlaneBitMask(planeName)
listToSet.append((idx, planeBitMask))

for plane, bitmask in listToSet:
if len(hdr) < plane:
continue

if hdr[plane].data is None:
continue

nmask = len(hdr[plane].data)
for i in range(nmask):
frow = hdr[plane].data[i]
Objmask(frow, bitmask).setMask(mask)

return mask

Expand Down
5 changes: 2 additions & 3 deletions python/lsst/obs/sdss/convertpsField.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ def convertpsField(infile, filt, trim=True, rcscale=0.001, MAX_ORDER_B=5, LSST_O
print("INVALID FILTER", filt)
sys.exit(1)

buff = open(infile, "rb")
pstruct = fits.getdata(buff, ext=filtToHdu[filt])
with open(infile, "rb") as buff:
pstruct = fits.getdata(buff, ext=filtToHdu[filt])

spaParList = [[]]*len(pstruct)
kernelList = []
Expand Down Expand Up @@ -124,7 +124,6 @@ def convertpsField(infile, filt, trim=True, rcscale=0.001, MAX_ORDER_B=5, LSST_O
spaParamsTri = spaParamsTri[:nTerms]
spaParList[i] = spaParamsTri

buff.close()
spaFun = afwMath.PolynomialFunction2D(LSST_ORDER)
spatialKernel = afwMath.LinearCombinationKernel(kernelList, spaFun)
spatialKernel.setSpatialParameters(spaParList)
Expand Down
24 changes: 11 additions & 13 deletions python/lsst/obs/sdss/converttsField.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,19 +46,19 @@ def converttsField(infile, filt, exptime=53.907456):
- exptime: exposure time (sec)
- airmass: airmass
"""
ptr = fits.open(infile)
if ptr[0].header['NFIELDS'] != 1:
print("INVALID TSFIELD FILE")
sys.exit(1)
filts = ptr[0].header['FILTERS'].split()
idx = filts.index(filt)
with fits.open(infile) as ptr:
if ptr[0].header['NFIELDS'] != 1:
print("INVALID TSFIELD FILE")
sys.exit(1)
filts = ptr[0].header['FILTERS'].split()
idx = filts.index(filt)

mjdTaiStart = ptr[1].data.field('mjd')[0][idx] # MJD(TAI) when row 0 was read
airmass = ptr[1].data.field("airmass")[0][idx]
mjdTaiStart = ptr[1].data.field('mjd')[0][idx] # MJD(TAI) when row 0 was read
airmass = ptr[1].data.field("airmass")[0][idx]

gain = float(ptr[1].data.field('gain')[0][idx]) # comes out as numpy.float32
aa = ptr[1].data.field('aa')[0][idx] # f0 = 10**(-0.4*aa) counts/second
aaErr = ptr[1].data.field('aaErr')[0][idx]
gain = float(ptr[1].data.field('gain')[0][idx]) # comes out as numpy.float32
aa = ptr[1].data.field('aa')[0][idx] # f0 = 10**(-0.4*aa) counts/second
aaErr = ptr[1].data.field('aaErr')[0][idx]

# Conversions
dateAvg = dafBase.DateTime(mjdTaiStart + 0.5 * exptime / 3600 / 24)
Expand All @@ -67,8 +67,6 @@ def converttsField(infile, filt, exptime=53.907456):

photoCalib = afwImage.makePhotoCalibFromCalibZeroPoint(fluxMag0, dfluxMag0)

ptr.close()

return TsField(
photoCalib=photoCalib,
gain=gain,
Expand Down
27 changes: 14 additions & 13 deletions python/lsst/obs/sdss/makeCamera.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
from lsst.afw.cameraGeom import makeCameraFromCatalogs, CameraConfig, DetectorConfig, \
TransformMapConfig, SCIENCE, PIXELS, FIELD_ANGLE, FOCAL_PLANE, NullLinearityType
import lsst.afw.table as afwTable
import lsst.geom as geom
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a stray import of afwGeom on line 25.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I see. This is needed for the TransformConfig and transformRegistry.

from lsst.obs.sdss.convertOpECalib import SdssCameraState

#
Expand Down Expand Up @@ -59,22 +60,22 @@ def addAmp(ampCatalog, i, eparams):
#
record = ampCatalog.addNew()
xtot = width + nExtended + nOverclock
allPixels = afwGeom.BoxI(afwGeom.PointI(0, 0), afwGeom.ExtentI(xtot, height))
biasSec = afwGeom.BoxI(afwGeom.PointI(nExtended if i == 0 else width, 0),
afwGeom.ExtentI(nOverclock, height))
dataSec = afwGeom.BoxI(afwGeom.PointI(nExtended + nOverclock if i == 0 else 0, 0),
afwGeom.ExtentI(width, height))
emptyBox = afwGeom.BoxI()
bbox = afwGeom.BoxI(afwGeom.PointI(0, 0), afwGeom.ExtentI(width, height))
bbox.shift(afwGeom.Extent2I(width*i, 0))

shiftp = afwGeom.Extent2I(xtot*i, 0)
allPixels = geom.BoxI(geom.PointI(0, 0), geom.ExtentI(xtot, height))
biasSec = geom.BoxI(geom.PointI(nExtended if i == 0 else width, 0),
geom.ExtentI(nOverclock, height))
dataSec = geom.BoxI(geom.PointI(nExtended + nOverclock if i == 0 else 0, 0),
geom.ExtentI(width, height))
emptyBox = geom.BoxI()
bbox = geom.BoxI(geom.PointI(0, 0), geom.ExtentI(width, height))
bbox.shift(geom.Extent2I(width*i, 0))

shiftp = geom.Extent2I(xtot*i, 0)
allPixels.shift(shiftp)
biasSec.shift(shiftp)
dataSec.shift(shiftp)

record.setBBox(bbox)
record.setRawXYOffset(afwGeom.ExtentI(0, 0))
record.setRawXYOffset(geom.ExtentI(0, 0))
record.setName('left' if i == 0 else 'right')
record.setReadoutCorner(afwTable.LL if i == 0 else afwTable.LR)
record.setGain(eparams['gain'])
Expand Down Expand Up @@ -153,7 +154,7 @@ def makeCamera(name="SDSS", outputDir=None):
camConfig.name = name
camConfig.detectorList = {}
camConfig.plateScale = 16.5 # arcsec/mm
pScaleRad = afwGeom.arcsecToRad(camConfig.plateScale)
pScaleRad = geom.arcsecToRad(camConfig.plateScale)
radialDistortCoeffs = [0.0, 1.0/pScaleRad]
tConfig = afwGeom.TransformConfig()
tConfig.transform.name = 'inverted'
Expand All @@ -172,7 +173,7 @@ def makeCamera(name="SDSS", outputDir=None):
filters = "riuzg"
for j, c in enumerate(reversed(filters)):
ccdName = "%s%s" % (c, dewarName)
offsetPoint = afwGeom.Point2D(25.4*2.5*(2.5-i), 25.4*2.1*(2.0 - j))
offsetPoint = geom.Point2D(25.4*2.5*(2.5-i), 25.4*2.1*(2.0 - j))
ccdInfo = makeCcd(ccdName, ccdId, offsetPoint)
ampInfoCatDict[ccdName] = ccdInfo['ampInfo']
camConfig.detectorList[ccdId] = ccdInfo['ccdConfig']
Expand Down
6 changes: 3 additions & 3 deletions python/lsst/obs/sdss/sdssNullIsr.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
import lsst.afw.image as afwImage
import lsst.afw.geom as afwGeom
import lsst.geom as geom
from lsst.pipe.tasks.processCcd import ProcessCcdTask


Expand Down Expand Up @@ -130,8 +130,8 @@ def loadExposure(self, sensorRef):
bbox = mi.getBBox()
begin = bbox.getBegin()
extent = bbox.getDimensions()
extent -= afwGeom.Extent2I(0, self.config.overlapSize)
tbbox = afwGeom.BoxI(begin, extent)
extent -= geom.Extent2I(0, self.config.overlapSize)
tbbox = geom.BoxI(begin, extent)
mi = afwImage.MaskedImageF(mi, tbbox)

exposure = afwImage.ExposureF(mi, wcs)
Expand Down
3 changes: 2 additions & 1 deletion tests/test_sdss.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,9 @@
import lsst.daf.persistence as dafPersist
from lsst.daf.base import DateTime
import lsst.afw.image
from lsst.afw.geom import SkyWcs, SpherePoint, degrees
from lsst.afw.geom import SkyWcs
import lsst.afw.detection
from lsst.geom import SpherePoint, degrees


class SdssMapperTestCase(lsst.utils.tests.TestCase):
Expand Down
1 change: 1 addition & 0 deletions ups/obs_sdss.table
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ setupRequired(pipe_tasks)
setupRequired(utils)
setupRequired(log)
setupRequired(astropy)
setupRequired(geom)

envPrepend(PYTHONPATH, ${PRODUCT_DIR}/python)
envPrepend(PATH, ${PRODUCT_DIR}/bin)