From a46cd79022f116ba14a74df22fcc178f029778f8 Mon Sep 17 00:00:00 2001 From: Thomas Williams <32936518+thomaswilliamsastro@users.noreply.github.com> Date: Mon, 21 Sep 2026 11:54:39 +0200 Subject: [PATCH] Begin removal of analysisUtils This PR starts the long-overdue removal of the analysisUtils dependency. Keeping this up-to-date is difficult, complicates tests, and makes keeping the pipeline transferable annoying. This first step removes the analysisUtils dependencies for everything except the SD pipeline. The reason the SD pipeline has been left out is because there are *a lot* of functions here that would need to be translated, both across the current version and the legacy pipeline. There's a decision to be made here about whether now is the time to deprecate the legacy pipeline; I think it is, but happy to hear opinions. In any case, the shared functions have been moved into utilsSingleDish, so it'll be super easy to rip it out in the future. Major change summary: - Removed the calls to au.createCasaTool(). Since we now use a unified CASA version, this is no longer needed. - Add a new custom estimate_synthesised_beam function. This has modifications to be better for large datasets by optionally using the first field in the dataset, and will loop over each observation ID to find a minimum beam size from each individual, concatenated observation. It's about twice as fast, and is quieter. - Add a new custom pick_cell_size function. This is a little different from the au equivalent. First, it'll always return the image size, and automatically pads for an FFT (we were actually doing this twice, in effect, so now the image sizes should be smaller). Second, it does not take in information about the target or the SPW, since we only ever call this at a point where we've already split out the dataset to a single target and window. - Add get_spw_info function, which is essentially a direct translation of the analysisUtils function, with a bit of tidying up - Add get_scans_for_spw function, which is essentially a direct translation of the analysisUtils function, with a bit of tidying up - Add get_science_spws function, which is essentially a direct translation of the analysisUtils function, with a bit of tidying up - Add get_chan_widths function, which is essentially a direct translation of the analysisUtils function, with a bit of tidying up - Add clear_pointing_table function, which is essentially a direct translation of the analysisUtils function, with a bit of tidying up - A couple of functions, get_sourcename and read_source_coordinates, have been moved from to casaLegacySingleDishRoutines to utilsSingleDish - Updated tests to remove those no longer needed Detailed changes: - Remove au.createCasaTool calls from casaCubeRoutines - Remove au.createCasaTool calls from casaFeatherRoutines - Remove au.createCasaTool calls from casaImagingRoutines - Replaced au.pickCellSize with pick_cell_and_im_size in casaImagingRoutines - Replaced au.ValueMapping call with get_spw_info in casaImagingRoutines - Remove au.createCasaTool calls from casaLegacySingleDishRoutines - Move get_sourcename and read_source_coordinates from casaLegacySingleDishRoutines to utilsSingleDish - Remove au.createCasaTool calls from casaMaskingRoutines - Remove au.createCasaTool calls from casaMosaicRoutines - Remove unused analysisUtils import from casaRoutineTests - Update imports for get_sourcename and read_source_coordinates in casaSingleDishALMAWrapper - Remove au.createCasaTool calls from casaVisRoutines - Add estimate_synthesised_beam function to casaVisRoutines - Add pick_cell_and_im_size function to casaVisRoutines - Add get_representative_freq function to casaVisRoutines - Add get_obs_ids function to casaVisRoutines - Add get_chan_widths function to casaVisRoutines - Add get_spw_info function to casaVisRoutines - Add get_scans_for_spw function to casaVisRoutines - Add get_science_spws function to casaVisRoutines - Add get_observatory_name function to casaVisRoutines - Add clear_pointing_table function to casaVisRoutines - Replace estimateSynthesizedBeam in clean_call - Remove unused analysisUtils import from handlerImaging - Remove au.createCasaTool calls from handlerImagingChunked - Remove au.createCasaTool calls from handlerTestImaging - Remove au.createCasaTool calls from utilsFieldSelection - Replaced getScienceSpws in utilsFieldSelection - Remove au.createCasaTool calls from utilsSingleDish - Remove au.createCasaTool calls from utilsTestImagingPlots - Update tests - Update CHANGELOG.md --- CHANGELOG.md | 1 + phangsPipeline/casaCubeRoutines.py | 19 +- phangsPipeline/casaFeatherRoutines.py | 9 +- phangsPipeline/casaImagingRoutines.py | 81 +- .../casaLegacySingleDishRoutines.py | 108 +-- phangsPipeline/casaMaskingRoutines.py | 15 +- phangsPipeline/casaMosaicRoutines.py | 21 +- phangsPipeline/casaRoutineTests.py | 1 - phangsPipeline/casaSingleDishALMAWrapper.py | 7 +- phangsPipeline/casaVisRoutines.py | 823 +++++++++++++++--- phangsPipeline/clean_call.py | 5 +- phangsPipeline/handlerImaging.py | 3 - phangsPipeline/handlerImagingChunked.py | 8 +- phangsPipeline/handlerTestImaging.py | 5 +- phangsPipeline/utilsFieldSelection.py | 14 +- phangsPipeline/utilsSingleDish.py | 91 +- phangsPipeline/utilsTestImagingPlots.py | 3 +- tests/test_casa_imports.py | 23 - 18 files changed, 863 insertions(+), 374 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ea5d18c3..162ecbc0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Use spectral-cube common beam calculations for common beams (#385). - Mosaic weights are now generated after convolution and regridding, and will undo any primary beam corrections (#396). - Don't skip making large mosaics (but keep the log messages to inform users this is happening) (#405). +- Begin removal of analysisUtils (#402). ### Fixed diff --git a/phangsPipeline/casaCubeRoutines.py b/phangsPipeline/casaCubeRoutines.py index bd99f3ae..37b41eee 100644 --- a/phangsPipeline/casaCubeRoutines.py +++ b/phangsPipeline/casaCubeRoutines.py @@ -7,7 +7,6 @@ import logging import os -import analysisUtils as au import numpy as np import scipy.ndimage as nd from astropy.io import fits @@ -45,7 +44,7 @@ def check_getchunk_putchunk_memory_issue( Args: infile: input image file. - myia: optional, if not given then we will run `myia = au.createCasaTool(casaStuff.iatool)` then `myia.open(infile)`. + myia: optional, if not given then we will run `myia = casaStuff.iatool()` then `myia.open(infile)`. return_data: if True then we will return the data if `getchunk()` runs properly withut a memory issue. return_mask: if True then we will return the mask if `getchunk(getmask=True)` runs properly withut a memory issue. return_shape: if True then we will return the shape of the cube data. @@ -64,7 +63,7 @@ def check_getchunk_putchunk_memory_issue( cube_mask = None has_opened_file = False if myia is None: - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(infile) has_opened_file = True if not myia.isopen(): @@ -177,14 +176,14 @@ def get_mask(infile, huge_cube_workaround=True): # os.system('rm -rf ' + infile + '.temp.fits') # #else: - # myia = au.createCasaTool(casaStuff.iatool) + # myia = casaStuff.iatool() # myia.open(infile+'.temp_deg') # mask = myia.getchunk(getmask=True) # myia.close() # #os.system('rm -rf ' + infile + '.temp_deg') - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(infile) @@ -237,7 +236,7 @@ def copy_mask(infile, outfile, huge_cube_workaround=True): # use putregion to update pixel mask - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(outfile) @@ -259,7 +258,7 @@ def copy_mask(infile, outfile, huge_cube_workaround=True): else: # getchunk was unsuccessful, has memory issue # putregion channel by channel logger.debug('putregion channel by channel for known memory issue') - myrg = au.createCasaTool(casaStuff.rgtool) + myrg = casaStuff.rgtool() nx = out_shape[0] ny = out_shape[1] #logger.debug('copy_mask mask.shape {}'.format(mask.shape)) @@ -324,13 +323,13 @@ def multiply_cube_by_value(infile, value, brightness_unit, huge_cube_workaround= # overwrite=True) # os.system('rm -rf ' + infile + '.fits') #else: - # myia = au.createCasaTool(casaStuff.iatool) + # myia = casaStuff.iatool() # myia.open(infile) # vals = myia.getchunk() # vals *= value # myia.putchunk(vals) - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(infile) @@ -627,7 +626,7 @@ def trim_rind( target_file = infile # Figure out the extent of the image inside the cube - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(target_file) mask = myia.getchunk(getmask=True) elt = nd.generate_binary_structure(2,1) diff --git a/phangsPipeline/casaFeatherRoutines.py b/phangsPipeline/casaFeatherRoutines.py index 5b79feb9..d8042273 100644 --- a/phangsPipeline/casaFeatherRoutines.py +++ b/phangsPipeline/casaFeatherRoutines.py @@ -6,7 +6,6 @@ import logging import os -import analysisUtils as au import numpy as np from . import casaCubeRoutines as ccr @@ -115,12 +114,12 @@ def prep_sd_for_feather( overwrite=overwrite) # Get out and match axis orders from interferometric images - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(interf_file) interf_cs = myia.coordsys() myia.close() - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(current_outfile) sd_cs = myia.coordsys() myia.close() @@ -235,7 +234,7 @@ def feather_two_cubes( os.system('cp -rf '+interf_file+' '+current_interf_file) os.system('cp -rf '+sd_file+' '+current_sd_file) - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() # As noted in [https://casa.nrao.edu/docs/casaref/image.putchunk.html], # "If all the pixels didn't easily fit in memory, you would iterate through @@ -390,7 +389,7 @@ def feather_two_cubes( # Finally, if we don't have a 4th (Stokes I) axis, add this # here to avoid feather crashing - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(current_interf_file) interf_shape = myia.shape() diff --git a/phangsPipeline/casaImagingRoutines.py b/phangsPipeline/casaImagingRoutines.py index 00cafda5..06bf37f5 100644 --- a/phangsPipeline/casaImagingRoutines.py +++ b/phangsPipeline/casaImagingRoutines.py @@ -7,11 +7,11 @@ import logging import os -import analysisUtils as au import numpy as np from . import casaMaskingRoutines as cmr from . import casaStuff +from .casaVisRoutines import get_spw_info, pick_cell_and_im_size from .clean_call import CleanCall logger = logging.getLogger(__name__) @@ -50,65 +50,21 @@ def estimate_cell_and_imsize( else: pblevel = clean_call.get_param('pblimit') - # These are the CASA-preferred sizes for fast FFTs - - valid_sizes = [] - for ii in range(10): - for kk in range(3): - for jj in range(3): - valid_sizes.append(2 ** (ii + 1) * 5 ** jj * 3 ** kk) - valid_sizes = sorted(valid_sizes) - valid_sizes = np.array(valid_sizes) - - # Cell size implied by baseline distribution from analysis - # utilities. - - au_cellsize, au_imsize, _ = au.pickCellSize(infile, - imsize=True, - npix=oversamp, - intent='', - pblevel=pblevel, - ) - xextent = au_cellsize * au_imsize[0] * 1.2 - yextent = au_cellsize * au_imsize[1] * 1.2 - - # Make the cell size a nice round number - - if au_cellsize < 0.01: - cell_size = au_cellsize - elif 0.01 <= au_cellsize < 0.05: - cell_size = np.floor(au_cellsize / 0.005) * 0.005 - elif 0.05 <= au_cellsize < 0.1: - cell_size = np.floor(au_cellsize / 0.01) * 0.01 - elif 0.1 <= au_cellsize < 0.5: - cell_size = np.floor(au_cellsize / 0.05) * 0.05 - elif 0.5 <= au_cellsize < 1.0: - cell_size = np.floor(au_cellsize / 0.1) * 0.1 - elif 1.0 <= au_cellsize < 2.0: - cell_size = np.floor(au_cellsize / 0.25) * 0.25 - elif 2.0 <= au_cellsize < 5.0: - cell_size = np.floor(au_cellsize / 0.5) * 0.5 - else: - cell_size = np.floor(au_cellsize / 1.0) * 0.5 - - # Now make the image size a good number for the FFT - - need_cells_x = xextent / cell_size - need_cells_y = yextent / cell_size - - cells_x = np.min(valid_sizes[valid_sizes > need_cells_x]) - cells_y = np.min(valid_sizes[valid_sizes > need_cells_y]) + # Cell size implied by baseline distribution + cell_size, image_size = pick_cell_and_im_size( + infile, + npix=oversamp, + pblevel=pblevel, + ) # If requested, force the mosaic to be square. This avoids # pathologies in CASA versions 5.1 and 5.3. - if force_square: - if cells_y < cells_x: - cells_y = cells_x - if cells_x < cells_y: - cells_x = cells_y + if image_size[0] < image_size[1]: + image_size[1] = image_size[0] + if image_size[0] > image_size[1]: + image_size[0] = image_size[1] - image_size = [int(cells_x), int(cells_y)] cell_size_string = str(cell_size) + 'arcsec' return cell_size_string, image_size @@ -786,18 +742,17 @@ def clean_loop( working_call.set_param('calcpsf', False) # Note the number of channels, which is used in setting the number - # of iterations that we give to an individual clean call. - - vm = au.ValueMapping(working_call.get_param('vis')) - - # If we have a number of channels defined, then use that + # of iterations that we give to an individual clean call. If we + # have a number of channels defined, then use that nchan = working_call.get_param('nchan') # Otherwise, pull the total number of channels from the spwInfo if nchan is None: nchan = -1 if nchan <= 0: - nchan = vm.spwInfo[0]['numChannels'] + vis = working_call.get_param('vis') + spw_info = get_spw_info(vis) + nchan = spw_info[0]['numChannels'] # Create a text record of progress through successive clean calls. @@ -1016,7 +971,7 @@ def clean_loop( pb_image = (_imagename + '.pb') if _imagename else None if pb_image and os.path.isdir(pb_image): try: - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(pb_image) pb_data = myia.getchunk() myia.close() @@ -1115,7 +1070,7 @@ def calc_residual_statistics( logger.error('Error! The input file "' + mask_name + '" was not found!') return - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(mask_name) mask = myia.getchunk() diff --git a/phangsPipeline/casaLegacySingleDishRoutines.py b/phangsPipeline/casaLegacySingleDishRoutines.py index 53d2548b..60f37115 100644 --- a/phangsPipeline/casaLegacySingleDishRoutines.py +++ b/phangsPipeline/casaLegacySingleDishRoutines.py @@ -49,7 +49,13 @@ import numpy as np from . import casaStuff -from .utilsSingleDish import getTPSampling, get_first_arr_val +from .casaVisRoutines import get_science_spws +from .utilsSingleDish import ( + getTPSampling, + get_first_arr_val, + get_sourcename, + read_source_coordinates, +) es = au.stuffForScienceDataReduction() @@ -86,19 +92,11 @@ def checkpipeline(): logger.info("> Setting the variable 'pipeline' to False") return False -# Creating CASA tools -#def createCasaTool(mytool): -# -# if (type(casac.Quantity) != type): # casa 4.x -# myt = mytool() -# else: # casa 3.x -# myt = mytool.create() -# return(myt) # Retrieve name of the column def getDataColumnName(inputms): - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(inputms) colnames = mytb.colnames() if 'FLOAT_DATA' in colnames: @@ -112,7 +110,7 @@ def getDataColumnForSDBaseline(vis): """ Returns the names of the corrected data columns (corrected) in a measurement set. """ - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(vis) names = copy.copy(mytb.colnames()) mytb.close() @@ -178,8 +176,8 @@ def scaleAutocorr(vis, scale=1., antenna='', spw='', field='', scan=''): logger.warning("No table.dat. This does not appear to be an MS.") return - mymsmd = au.createCasaTool(casaStuff.msmdtool) - mytb = au.createCasaTool(casaStuff.tbtool) + mymsmd = casaStuff.msmdtool() + mytb = casaStuff.tbtool() conditions = ["ANTENNA1==ANTENNA2"] @@ -242,7 +240,7 @@ def scaleAutocorr(vis, scale=1., antenna='', spw='', field='', scan=''): # Create vector with antenna names def read_ants_names(filename): - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(filename + '/ANTENNA') vec_ants = mytb.getcol('NAME') mytb.close() @@ -266,7 +264,7 @@ def get_tsysmap(tsysmap,spws_scie,spws_tsys,freq_rep_scie,freq_rep_tsys): def read_spw(filename,source): # Tsys spws (index) - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(filename + '/SYSCAL') spwstsys = mytb.getcol('SPECTRAL_WINDOW_ID') @@ -285,7 +283,7 @@ def read_spw(filename,source): names = mytb.getcol('NAME') rest_freq_scie = [rest_freq_scie[i] for i in range(len(spws_scie)) if "FULL_RES" in names[spws_scie[i]]] spws_scie = [spw for spw in spws_scie if "FULL_RES" in names[spw]] - spws_scie = au.getScienceSpws(filename) + spws_scie = get_science_spws(filename) spws_scie = spws_scie.split(",") spws_scie = [int(i) for i in spws_scie] @@ -315,7 +313,7 @@ def read_spw(filename,source): # Get information of the source velocity def read_vel_source(filename,source): - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(filename + '/SOURCE') names = mytb.getcol('NAME') numli = mytb.getcol('NUM_LINES') @@ -456,7 +454,7 @@ def convert_vel2chan_line(filename_in,freq_rest,vel_line,spw_line,coords,date): freq2_topo = au.lsrkToTopo(freq2, date, ra, dec) if fsuffix == '.asap': - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(filename_in) nchan = mytb.getkeyword('nChan') if_eq = mytb.getcol('FREQ_ID',startrow=1,nrow=1) @@ -470,7 +468,7 @@ def convert_vel2chan_line(filename_in,freq_rest,vel_line,spw_line,coords,date): mytb.close() else: - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(filename_in+os.sep+'SPECTRAL_WINDOW') spwids = np.arange(mytb.nrows()) spw_line = -1 @@ -581,73 +579,6 @@ def extract_jyperk(filename, spw_line, pipeline): jyperk = kw['jyperk'] return jyperk -# Read source coordinates -def read_source_coordinates(filename,source): - - coord_source = au.getRADecForSource(filename,source) - RA_h = (coord_source.split(' ')[0]).split(':')[0] - RA_m = (coord_source.split(' ')[0]).split(':')[1] - RA_s = (coord_source.split(' ')[0]).split(':')[2] - DEC_d = (coord_source.split(' ')[1]).split(':')[0] - DEC_m = (coord_source.split(' ')[1]).split(':')[1] - DEC_s = (coord_source.split(' ')[1]).split(':')[2] - coord = "J2000 "+str(RA_h)+"h"+str(RA_m)+"m"+str(RA_s[0:6])+" "+str(DEC_d)+"d"+str(DEC_m)+"m"+str(DEC_s) - return coord - -# Get source name -def get_sourcename(filename, - source='all', - ): - - mytb = au.createCasaTool(casaStuff.msmdtool) - mytb.open(filename) - - # Get source field names - fieldnames = mytb.fieldnames() - - # If we're looking for all intents, then just take the (first) science target - if source == 'all': - field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')[0] - - # Else, key in on the source name. Do this in a regex-style - else: - - # For field selection, this can be done with forward slashes, - # so make sure to strip them - if source.startswith("/"): - source = source.strip("/") - - # Get science fieldnames as names - fieldnames_asname = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE', asnames=True) - - regex_matches = [re.fullmatch(source, x) for x in fieldnames_asname] - regex_idxs = [x is not None for x in regex_matches] - found_field = any(regex_idxs) - - # If we've found a match, update the source name with the first - # fieldname - if found_field: - source = fieldnames_asname[regex_idxs][0] - - if not found_field: - logger.warning(f"source {source} not in field names: {', '.join(fieldnames)}") - return None - - field = mytb.fieldsforname(source) - - # If we don't find anything, get the science targets - if len(field) == 0: - field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE') - - # Take the first result - field = field[0] - - source_name = fieldnames[field] - - mytb.close() - - return source_name - # Create string of spws to apply the Tsys def str_spw_apply_tsys(spws_info): #science spws @@ -1491,7 +1422,7 @@ def imaging(source, name_line, phcenter, vel_source, source_vel_kms, vwidth_kms, infile = ms_concat_filename # Read frequency - mymsmd = au.createCasaTool(casaStuff.msmdtool) + mymsmd = casaStuff.msmdtool() mymsmd.open(infile) freq = mymsmd.meanfreq(0) mymsmd.close() @@ -1585,7 +1516,7 @@ def imaging(source, name_line, phcenter, vel_source, source_vel_kms, vwidth_kms, #ia.open('ALMA_TP.'+source+'.'+name_line+suffix+'.image') #ia.setrestoringbeam(major = str(sfbeam)+'arcsec', minor = str(sfbeam)+'arcsec', pa = '0deg') #ia.done() - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(outimage) # 'ALMA_TP.'+source+'.'+name_line+suffix+'.image' myia.setrestoringbeam(major = str(sfbeam)+'arcsec', minor = str(sfbeam)+'arcsec', pa = '0deg') mask = myia.getchunk(getmask=True) @@ -1933,4 +1864,3 @@ def run_ALMA_TP_tools( os.chdir(ori_path) - diff --git a/phangsPipeline/casaMaskingRoutines.py b/phangsPipeline/casaMaskingRoutines.py index d59c996f..828ad1d1 100644 --- a/phangsPipeline/casaMaskingRoutines.py +++ b/phangsPipeline/casaMaskingRoutines.py @@ -6,7 +6,6 @@ import logging import os -import analysisUtils as au import numpy as np import scipy.ndimage as ndimage from astropy.io import fits @@ -162,7 +161,7 @@ def noise_for_cube( logger.error('maskfile specified but not found - ' + maskfile) return (None) - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(infile) has_memory_issue, cube_shape = ccr.check_getchunk_putchunk_memory_issue( @@ -191,7 +190,7 @@ def noise_for_cube( myia_mask = None if maskfile is not None: - myia_mask = au.createCasaTool(casaStuff.iatool) + myia_mask = casaStuff.iatool() myia_mask.open(maskfile) per_channel_noise = [] @@ -303,7 +302,7 @@ def read_cube(infile, huge_cube_workaround=True): # Remove intermediate fits file os.system('rm -rf ' + infile + '.fits') else: - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(infile) cube = myia.getchunk() myia.close() @@ -344,7 +343,7 @@ def write_mask(infile, outfile, mask, huge_cube_workaround=True): # Remove the intermediate fits file os.system('rm -rf ' + outfile + '.fits') else: - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(outfile) myia.putchunk(mask) myia.close() @@ -478,7 +477,7 @@ def apply_additional_mask( beam based masks by setting the PB file to new_mask_file and the pb_limit as new_thresh. """ - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(new_mask_file) new_mask = myia.getchunk() myia.close() @@ -515,8 +514,8 @@ def import_and_align_mask( overwrite=True) # Prepare analysis utility tool - myia = au.createCasaTool(casaStuff.iatool) - myim = au.createCasaTool(casaStuff.imtool) + myia = casaStuff.iatool() + myim = casaStuff.imtool() # Read mask data # myia.open(out_file+'.temp_copy') diff --git a/phangsPipeline/casaMosaicRoutines.py b/phangsPipeline/casaMosaicRoutines.py index 2d90f2ed..d103e3b6 100644 --- a/phangsPipeline/casaMosaicRoutines.py +++ b/phangsPipeline/casaMosaicRoutines.py @@ -7,7 +7,6 @@ import logging import os -import analysisUtils as au import numpy as np from . import casaMaskingRoutines as cma @@ -760,7 +759,7 @@ def generate_weight_file( # Copy the template and open it os.system("cp -r "+template+" "+outfile) - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(outfile) # Check for memory issues @@ -1004,7 +1003,7 @@ def mosaic_aligned_data( # We might need to add a 4th (Stokes) axis # to prevent crashing - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() for this_infile in infile_list: @@ -1142,8 +1141,8 @@ def mosaic_aligned_data( os.system("cp -r "+local_imlist[0]+" "+sum_file) os.system("cp -r "+local_imlist[0]+" "+weight_file) - myia_sum = au.createCasaTool(casaStuff.iatool) - myia_weight = au.createCasaTool(casaStuff.iatool) + myia_sum = casaStuff.iatool() + myia_weight = casaStuff.iatool() myia_sum.open(sum_file) myia_weight.open(weight_file) @@ -1171,7 +1170,7 @@ def mosaic_aligned_data( casaStuff.immath(imagename=temp_chan_images, mode='evalexpr', expr=lel_exp_sum, outfile=temp_sum) - myia_temp = au.createCasaTool(casaStuff.iatool) + myia_temp = casaStuff.iatool() myia_temp.open(temp_sum) sum_slice = myia_temp.getchunk() myia_temp.close() @@ -1237,7 +1236,7 @@ def mosaic_aligned_data( casaStuff.immath(imagename=temp_chan_images, mode='evalexpr', expr=lel_exp_sum, outfile=temp_sum) - myia_temp = au.createCasaTool(casaStuff.iatool) + myia_temp = casaStuff.iatool() myia_temp.open(temp_sum) sum_slice = myia_temp.getchunk() myia_temp.close() @@ -1267,7 +1266,7 @@ def mosaic_aligned_data( # Just to be safe, reset the masks on the two images. - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(sum_file) myia.set(pixelmask=1) myia.close() @@ -1298,8 +1297,8 @@ def mosaic_aligned_data( os.system("cp -r "+sum_file+" "+temp_file) os.system("cp -r "+weight_file+" "+local_maskfile) - myia_temp = au.createCasaTool(casaStuff.iatool) - myia_mask = au.createCasaTool(casaStuff.iatool) + myia_temp = casaStuff.iatool() + myia_mask = casaStuff.iatool() myia_temp.open(temp_file) myia_mask.open(local_maskfile) @@ -1368,7 +1367,7 @@ def mosaic_aligned_data( myia_mask.close() # Add in potential Stokes axis to the mask file - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(local_maskfile) in_shape = myia.shape() if len(in_shape) != 4: diff --git a/phangsPipeline/casaRoutineTests.py b/phangsPipeline/casaRoutineTests.py index 0f149c10..8cfd67e2 100644 --- a/phangsPipeline/casaRoutineTests.py +++ b/phangsPipeline/casaRoutineTests.py @@ -8,7 +8,6 @@ import logging -import analysisUtils as au import numpy as np from . import casaMaskingRoutines as cma diff --git a/phangsPipeline/casaSingleDishALMAWrapper.py b/phangsPipeline/casaSingleDishALMAWrapper.py index d5ea31e4..3a2f40b5 100644 --- a/phangsPipeline/casaSingleDishALMAWrapper.py +++ b/phangsPipeline/casaSingleDishALMAWrapper.py @@ -8,9 +8,8 @@ import numpy as np from astropy.table import Table -from . import casaLegacySingleDishRoutines as csdr from . import casaStuff -from .utilsSingleDish import getTPSampling +from .utilsSingleDish import getTPSampling, get_sourcename, read_source_coordinates logger = logging.getLogger(__name__) @@ -87,7 +86,7 @@ def SDImaging(filename, imsize = int(round(maxsize/cell)*1.5) if phcenter == '': - coord_phase = csdr.read_source_coordinates(filename, source) + coord_phase = read_source_coordinates(filename, source) else: coord_phase = phcenter @@ -374,7 +373,7 @@ def runALMAPipeline(path_galaxy, this_vis = ms_concat_filename # read the source name directly from the ms - source = csdr.get_sourcename(this_vis, source=in_source) + source = get_sourcename(this_vis, source=in_source) for this_product in product_dict: diff --git a/phangsPipeline/casaVisRoutines.py b/phangsPipeline/casaVisRoutines.py index 744e5e96..9242f8e5 100644 --- a/phangsPipeline/casaVisRoutines.py +++ b/phangsPipeline/casaVisRoutines.py @@ -7,7 +7,6 @@ import os import shutil -import analysisUtils as au import astropy.constants as const import astropy.units as u import numpy as np @@ -234,7 +233,7 @@ def split_science_targets( # Verify the column to use. If present, we use the corrected # column. If not, then we use the data column. - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(infile, nomodify = True) colnames = mytb.colnames() if 'CORRECTED_DATA' in colnames: @@ -571,14 +570,13 @@ def find_spws_for_line( spw_lowest_ghz = None spw_highest_ghz = None - logger.debug("... vm = au.ValueMapping(infile) ...") - vm = au.ValueMapping(infile) - logger.debug("... vm = au.ValueMapping(infile) done") + spw_info = get_spw_info(infile) + scans_for_spw = get_scans_for_spw(infile) - for this_spw in vm.spwInfo.keys(): + for this_spw in spw_info.keys(): - spw_high_ghz = np.max(vm.spwInfo[this_spw]['edgeChannels'])/1e9 - spw_low_ghz = np.min(vm.spwInfo[this_spw]['edgeChannels'])/1e9 + spw_high_ghz = np.max(spw_info[this_spw]['edgeChannels'])/1e9 + spw_low_ghz = np.min(spw_info[this_spw]['edgeChannels'])/1e9 logger.debug( "... spw: %s, freq: %.6f - %.6f GHz" % (this_spw, spw_low_ghz, spw_high_ghz)) @@ -590,12 +588,12 @@ def find_spws_for_line( continue if max_chanwidth_ghz is not None: - spw_chanwidth_ghz = abs(vm.spwInfo[this_spw]['chanWidth'])/1e9 + spw_chanwidth_ghz = abs(spw_info[this_spw]['chanWidth'])/1e9 if spw_chanwidth_ghz > max_chanwidth_ghz: continue if require_data: - if len(vm.scansForSpw[this_spw]) == 0: + if len(scans_for_spw[this_spw]) == 0: continue if require_full_line_coverage and not ( @@ -647,8 +645,7 @@ def find_spws_for_science( infile=None, require_data=False, exit_on_error=True, as_list=False): """ List all spectral windows that we judge likely to be used for - science. Mostly wraps analysisUtils rather than reinventing the - wheel. + science. """ # Check inputs @@ -667,25 +664,26 @@ def find_spws_for_science( 'Error! The input uv data measurement set "'+infile + '"does not exist!') - # Call the analysisUtil version. - - spw_string = au.getScienceSpws( - infile, intent='OBSERVE_TARGET*') + # Get science SPWs + spw_string = get_science_spws( + vis=infile, + intent='OBSERVE_TARGET*', + ) if spw_string is None or len(spw_string) == 0: - spw_string = au.getScienceSpws( - infile, intent='OBSERVE_TARGET#ON_SOURCE') + spw_string = get_science_spws( + vis=infile, + intent='OBSERVE_TARGET#ON_SOURCE', + ) spw_list = [] for this_spw_string in spw_string.split(','): spw_list.append(int(this_spw_string)) - # Shouldn't get here, I think, because of the analysisUtils logic - if require_data: - vm = au.ValueMapping(infile) + scans_for_spw = get_scans_for_spw(infile) for spw in spw_list: - if len(vm.scansForSpw[spw]) == 0: + if len(scans_for_spw[spw]) == 0: spw_list.remove(spw) # Return @@ -740,18 +738,18 @@ def spw_string_for_freq_ranges( if not isinstance(just_spw, list): just_spw = [just_spw] - vm = au.ValueMapping(infile) + spw_info = get_spw_info(infile) # Loop over spectral windows spw_flagging_string = '' first_string = True - for this_spw in vm.spwInfo: + for this_spw in spw_info: if len(just_spw) > 0: if this_spw not in just_spw: continue - freq_axis = vm.spwInfo[this_spw]['chanFreqs'] + freq_axis = spw_info[this_spw]['chanFreqs'] half_chan = abs(freq_axis[1]-freq_axis[0])*0.5 chan_axis = np.arange(len(freq_axis)) mask_axis = np.zeros_like(chan_axis, dtype='bool') @@ -843,7 +841,7 @@ def compute_common_chanwidth( this_infile, line, vsys_kms=vsys_kms, vwidth_kms=vwidth_kms, require_full_line_coverage=require_full_line_coverage) - chan_widths_hz = au.getChanWidths(this_infile, spw_list_string) + chan_widths_hz = get_chan_widths(this_infile, spw_list_string) # Convert to km/s and return for this_chan_width_hz in chan_widths_hz: @@ -951,10 +949,10 @@ def batch_extract_line( # This zaps the whole table: if os.path.exists(this_outfile+os.sep+'POINTING'): - au.clearPointingTable(this_outfile) + clear_pointing_table(this_outfile) else: copy_pointing = False - #logger.debug('Warning! Failed to run au.clearPointingTable(%r)'%(this_outfile)) + #logger.debug('Warning! Failed to run clear_pointing_table(%r)'%(this_outfile)) # Allow a small tolerance in the channel width if allow_freqtol_chanfrac: @@ -1045,7 +1043,7 @@ def suggest_extraction_scheme( for this_infile in infile_list: - vm = au.ValueMapping(this_infile) + spw_info = get_spw_info(this_infile) spw_list = find_spws_for_line( this_infile, restfreq_ghz=restfreq_ghz, @@ -1058,7 +1056,7 @@ def suggest_extraction_scheme( for this_spw in spw_list: - chan_width_ghz = np.abs(vm.spwInfo[this_spw]['chanWidth'])/1e9 + chan_width_ghz = np.abs(spw_info[this_spw]['chanWidth'])/1e9 # Using RADIO convention: chan_width_kms = chan_width_ghz/restfreq_ghz * sol_kms @@ -1072,7 +1070,7 @@ def suggest_extraction_scheme( # Figure out the binfactor this_binfactor = int(np.floor(target_chan_kms/chan_width_kms)) # clamp to nchan if the binfactor exceeds the number of channels in the spw - nchan_spw = vm.spwInfo[this_spw]['numChannels'] + nchan_spw = spw_info[this_spw]['numChannels'] if this_binfactor > nchan_spw: this_binfactor = nchan_spw @@ -1536,7 +1534,7 @@ def build_mstransform_call( # Determine the column to use if datacolumn is None: - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(infile, nomodify = True) colnames = mytb.colnames() if 'CORRECTED_DATA' in colnames: @@ -1596,7 +1594,7 @@ def build_mstransform_call( vlow_kms=vstart_kms, vhigh_kms=vstart_kms+vwidth_kms) # line_freq_ghz = (line_low_ghz+line_high_ghz)*0.5 - max_chan_ghz = np.max(np.abs(au.getChanWidths(infile, spw)))/1e9 + max_chan_ghz = np.max(np.abs(get_chan_widths(infile, spw)))/1e9 # Using RADIO velocity convention: current_chan_kms = max_chan_ghz/restfreq_ghz*sol_kms @@ -1740,7 +1738,7 @@ def reweight_data( # Determine column to use if datacolumn is None: - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(infile, nomodify = True) colnames = mytb.colnames() if 'CORRECTED_DATA' in colnames: @@ -1757,15 +1755,15 @@ def reweight_data( exclude_str = '' if (edge_chans is not None) or (edge_kms is not None): - vm = au.ValueMapping(infile) + spw_info = get_spw_info(infile) first = True - for this_spw in vm.spwInfo.keys(): + for this_spw in spw_info.keys(): if edge_kms is not None: - spw_high_ghz = np.max(vm.spwInfo[this_spw]['edgeChannels'])/1e9 - spw_low_ghz = np.min(vm.spwInfo[this_spw]['edgeChannels'])/1e9 - spw_chanwidth_ghz = abs(vm.spwInfo[this_spw]['chanWidth'])/1e9 + spw_high_ghz = np.max(spw_info[this_spw]['edgeChannels'])/1e9 + spw_low_ghz = np.min(spw_info[this_spw]['edgeChannels'])/1e9 + spw_chanwidth_ghz = abs(spw_info[this_spw]['chanWidth'])/1e9 mean_freq_ghz = 0.5*(spw_high_ghz+spw_low_ghz) @@ -1787,7 +1785,7 @@ def reweight_data( edge_chans = int(np.ceil(edge_kms / mean_chanwidth_kms)) - nchan = vm.spwInfo[this_spw]['numChannels'] + nchan = spw_info[this_spw]['numChannels'] if edge_chans*2 > nchan: logger.warning( @@ -1999,7 +1997,7 @@ def batch_extract_continuum( # all SPWs except the first one. if clear_pointing: - au.clearPointingTable(this_outfile) + clear_pointing_table(this_outfile) # Concatenate and combine the output data sets @@ -2107,7 +2105,7 @@ def extract_continuum( if not os.path.isdir(outfile + '.touch'): os.mkdir(outfile + '.touch') - mytb = au.createCasaTool(casaStuff.tbtool) + mytb = casaStuff.tbtool() mytb.open(infile + '.temp_copy', nomodify=True) colnames = mytb.colnames() if 'CORRECTED_DATA' in colnames: @@ -2316,9 +2314,9 @@ def noise_spectrum( return None # Note the number of channels in SPW 0 + spw_info = get_spw_info(vis) - vm = au.ValueMapping(vis) - nchan = vm.spwInfo[0]['numChannels'] + nchan = spw_info[0]['numChannels'] spec = np.zeros(nchan) for ii in range(nchan): if start_chan is not None: @@ -2402,89 +2400,12 @@ def estimate_mrs( tb_subset.close() - tb_subset = tb.query(" && ".join(subset_conditions)) - total_rows = tb_subset.nrows() - - uv_distance_m = [] - - # Loop over in chunks to reduce memory cost - for startrow in range(0, total_rows, chunk_size): - nrow = min(chunk_size, total_rows - startrow) - - uvw = np.asarray( - tb_subset.getcol("UVW", - startrow=startrow, - nrow=nrow, - rowincr=1, - ) - ) - antenna1 = np.asarray( - tb_subset.getcol( - "ANTENNA1", - startrow=startrow, - nrow=nrow, - rowincr=1, - ) - ) - antenna2 = np.asarray( - tb_subset.getcol( - "ANTENNA2", - startrow=startrow, - nrow=nrow, - rowincr=1, - ) - ) - ddid = np.asarray( - tb_subset.getcol( - "DATA_DESC_ID", - startrow=startrow, - nrow=nrow, - rowincr=1, - ) - ) - flag_row = np.asarray( - tb_subset.getcol( - "FLAG_ROW", - startrow=startrow, - nrow=nrow, - rowincr=1, - ), - dtype=bool, - ) - flags = np.asarray( - tb_subset.getcol( - "FLAG", - startrow=startrow, - nrow=nrow, - rowincr=1, - ), - dtype=bool, - ) - - # CASA normally returns UVW as (3, nrow). - if uvw.shape[0] != 3 and uvw.shape[-1] == 3: - uvw = uvw.T - - uvdm = np.hypot(uvw[0], uvw[1]) - - # FLAG normally has dimensions (ncorr, nchan, nrow). - # Keep a row if at least one correlation/channel is unflagged. - flag_axes = tuple(range(flags.ndim - 1)) - completely_flagged = np.all(flags, axis=flag_axes) - - valid = ( - ~flag_row - & ~completely_flagged - & (antenna1 != antenna2) - & np.isfinite(uvdm) - & (uvdm > 0) - & (ddid >= 0) - ) - - # Only take valid values - uv_distance_m.extend(uvdm[valid]) - - uv_distance_m = np.asarray(uv_distance_m) + # Get UV distances in meters + uv_distance_m = get_uv_distance_m( + tb=tb, + subset_conditions=subset_conditions, + chunk_size=chunk_size, + ) # Now we loop over, take the 5th percentile baseline baseline_percentiles = [] @@ -2516,6 +2437,218 @@ def estimate_mrs( return result +def estimate_synthesised_beam( + vis: str, + baseline_percentile: float = 80.0, + beam_factor: float = 0.574, + use_first_field: bool = True, + chunk_size: int = 100_000, +) -> float: + """ + Estimate the synthesised beam for a measurement set. + + This function is specifically designed to sidestep the biases that can arise + from concatenating measurement sets. It calculates some minimum baseline percentile + from each unique observation ID, and then uses the minimum of these to for the beam size. + By default, this function uses 0.574 * 80th baseline percentile, to match the analysisUtils default. + + Args: + vis (str): Path to measurement set + baseline_percentile (float): The percentile of the baseline distribution to use. + Defaults to 80. + beam_factor (float): Factor to convert from baseline percentile to synthesised beam. + Defaults to 0.574. + use_first_field (bool): Whether to only use the first field. + Defaults to True. + chunk_size (int): Number of rows to process in each chunk, to limit memory usage. + Defaults to 100,000. + + Returns: + float: The estimated synthesised beam in arcseconds + """ + + # Get the representative frequency + rep_freq = get_representative_freq(vis) + logger.debug(f"Representative frequency: {rep_freq}") + + # Get a list of the observation IDs + obs_ids = get_obs_ids(vis=vis) + logger.debug(f"Found {len(obs_ids)} observation IDs") + + tb = casaStuff.tbtool() + tb.open(vis) + + # Set up a dictionary to hold all the calculations + result = { + "representative_frequency": rep_freq.to(u.GHz).value, + "synthesised_beam_per_obs_id": {}, + } + + for obs_id in obs_ids: + + # Downselect on observation ID + subset_conditions = [ + f"OBSERVATION_ID == {obs_id}", + ] + tb_subset = tb.query(" && ".join(subset_conditions)) + + # If we're only using first field, get the first field ID for further downselecting + if use_first_field: + field_id = tb_subset.getcell("FIELD_ID", 0) + subset_conditions.append(f"FIELD_ID == {field_id}") + + tb_subset.close() + + # Get UV distances in meters + uv_distance_m = get_uv_distance_m( + tb=tb, + subset_conditions=subset_conditions, + chunk_size=chunk_size, + ) + baseline_for_beam = np.nanpercentile(uv_distance_m, baseline_percentile) * u.m + + # Calculate the synthesised beam in arcsec + synthesised_beam = ( + beam_factor + * const.c.to(u.m * u.Hz) + / (rep_freq.to(u.Hz) * baseline_for_beam.to(u.m)) + * u.rad + ) + synthesised_beam = synthesised_beam.to(u.arcsec).value + + result["synthesised_beam_per_obs_id"][obs_id] = synthesised_beam + + logger.debug(f"Calculated synthesised beam of {synthesised_beam} for obs ID {obs_id}") + + tb.close() + + # The synthesised beam is the minimum of the synthesised beams from each observation ID + synthesised_beam = np.nanmin([r[-1] for r in result["synthesised_beam_per_obs_id"].items()]) + + return synthesised_beam + + +def pick_cell_and_im_size( + vis: str, + npix: float = 5.0, + baseline_percentile: float = 80.0, + cellstring: bool = False, + roundcell: int = 2, + pblevel: float = 0.2, + beam_size_factor: float = 0.574, + beam_size_use_first_field: bool = True, + beam_size_chunk_size: int = 100_000, +) -> tuple[float | str, list[int]]: + """Pick a cell size and image size for a measurement set. + + Args: + vis (str): Path to measurement set. + npix (float): Number of pixels across the synthesised beam. + Defaults to 5.0. + baseline_percentile (float): The percentile of the baseline distribution to use for estimating the synthesised + beam. Defaults to 80.0. + cellstring (bool): Whether to return the cell size as a string with units. + Defaults to False. + roundcell (int): Number of decimal places to round the cell size to. + Defaults to 2. + pblevel (float): The primary beam level for calculating image size. + Defaults to 0.2. + beam_size_factor (float): Factor to convert from baseline percentile to synthesised beam. + Defaults to 0.574. + beam_size_use_first_field (bool): Whether to only use the first field for estimating the synthesised beam. + Defaults to True. + beam_size_chunk_size (int): Number of rows to process in each chunk when estimating the synthesised beam. + Defaults to 100,000. + + Returns: + float or str, list: Tuple of the cell size (potentially as a string) and an associated list of imsize. + """ + + # Get cell size from synthesised beam and oversample factor + synthesised_beam = estimate_synthesised_beam( + vis=vis, + baseline_percentile=baseline_percentile, + beam_factor=beam_size_factor, + use_first_field=beam_size_use_first_field, + chunk_size=beam_size_chunk_size, + ) + cell_size = synthesised_beam / npix + + # If selected, round to requested number of significant figures + if roundcell > 0: + cell_size = float(f"{cell_size:.{roundcell}g}") + + # Convert to a cell size with units + cell_size_unit = cell_size * u.arcsec + + if cellstring: + cell_size = f"{cell_size}arcsec" + + # Get dish diameter for image size calculation. Use max to get minimum FOV + tb = casaStuff.table() + tb.open(os.path.join(f"{vis}/ANTENNA")) + dish_dia = tb.getcol("DISH_DIAMETER").max() + dish_dia = dish_dia * u.m + tb.close() + + # Get representative frequency + rep_freq = get_representative_freq(vis) + + # Get the source name from the MS metadata + msmd = casaStuff.msmdtool() + msmd.open(vis) + + # Get a list of all science field IDs, turn into a source name by taking the first + # in the MS + intents = msmd.intents() + field_ids = msmd.fieldsforintent(intents[0]) + sourcename = msmd.namesforfields(field_ids)[0] + logger.debug(f"Using source name {sourcename} for image size calculation") + msmd.close() + + tb = casaStuff.table() + tb.open(f"{vis}/FIELD") + tb_subset = tb.query(f"NAME == '{sourcename}'") + phase_dirs = tb_subset.getcol("PHASE_DIR") + + tb_subset.close() + tb.close() + + # Extract RA and Dec arrays (in radians) + ra_fields = phase_dirs[0, 0, :] * u.rad + dec_fields = phase_dirs[1, 0, :] * u.rad + + # Get primary beam FWHM + fwhm = (1.14 * 1.22 * const.c.to(u.m * u.Hz) / rep_freq.to(u.Hz) / dish_dia) * u.rad + + # Apply the PB cutoff + radius_pb = (fwhm / 2.0) * np.sqrt(np.log(pblevel) / np.log(0.5)) + + ra_mins = ra_fields - radius_pb + ra_maxs = ra_fields + radius_pb + dec_mins = dec_fields - radius_pb + dec_maxs = dec_fields + radius_pb + + ra_extent = ra_maxs.max() - ra_mins.min() + dec_extent = dec_maxs.max() - dec_mins.min() + + x_extent = ra_extent / cell_size_unit + y_extent = dec_extent / cell_size_unit + + # Round this up to a good FFT number for the FFT, + # cast to int + x_extent = int(1.2 * x_extent) + y_extent = int(1.2 * y_extent) + + su = casaStuff.synthesisutils() + x_extent = int(su.getOptimumSize(x_extent)) + y_extent = int(su.getOptimumSize(y_extent)) + su.done() + + im_size = [x_extent, y_extent] + + return cell_size, im_size + def get_representative_freq( vis: str, ) -> u.Quantity: @@ -2592,3 +2725,439 @@ def get_obs_ids(vis): obs_ids = list(range(nrows)) return obs_ids + + +def get_chan_widths( + vis: str, + spw: str | int | list = "", + velocity: bool = False, +) -> np.ndarray: + """Get channel widths for a given measurement set and SPW selection. + + Args: + vis (str): Path to measurement set. + spw (str | int | list, optional): SPW selection. Can be a string, integer, or list of integers. + Defaults to "" (all SPWs). + velocity (bool, optional): If True, return channel widths in km/s. + Defaults to False. + + Returns: + np.ndarray: Array of channel widths for the selected SPWs. If velocity is True, + the widths are in velocity units; otherwise, they are in Hz. + """ + + if not os.path.exists(vis): + raise FileNotFoundError(f"Measurement set '{vis}' not found.") + + msmd = casaStuff.msmdtool() + msmd.open(vis) + + # If an empty string or list is passed, selected all SPWs + if spw in ["", []]: + spw = list(range(msmd.nspw())) + + # If we still have a string, turn to int + if isinstance(spw, str): + spw = [int(i) for i in spw.split(",")] + + # Make sure spw is a list + if not isinstance(spw, list): + spw = [spw] + + chan_widths = [] + mean_freqs = [] + for s in spw: + chan_width = np.median(msmd.chanwidths(s)) + mean_freq = msmd.meanfreq(s) + + chan_widths.append(chan_width) + mean_freqs.append(mean_freq) + + chan_widths = chan_widths * u.Hz + + # Convert to velocity units if requested + if velocity: + mean_freqs = mean_freqs * u.Hz + chan_widths *= const.c.to(u.m * u.Hz) / mean_freqs + chan_widths = chan_widths.to(u.km / u.s).value + else: + chan_widths = chan_widths.to(u.Hz).value + + msmd.close() + + return chan_widths + +def get_uv_distance_m( + tb: casaStuff.tbtool, + subset_conditions: list, + chunk_size: int = 100_000, +) -> np.ndarray: + """Get the UV distance in meters for a given table and subset conditions. + + Args: + tb (casaStuff.tbtool): The CASA table tool object. + subset_conditions (list): List of conditions to subset the table. + chunk_size (int): Number of rows to process in each chunk, to limit memory usage. + + Returns: + np.ndarray: Array of UV distances in meters. + """ + + tb_subset = tb.query(" && ".join(subset_conditions)) + total_rows = tb_subset.nrows() + + uv_distance_m = [] + + # Loop over in chunks to reduce memory cost + for startrow in range(0, total_rows, chunk_size): + nrow = min(chunk_size, total_rows - startrow) + + uvw = np.asarray( + tb_subset.getcol( + "UVW", + startrow=startrow, + nrow=nrow, + rowincr=1, + ) + ) + antenna1 = np.asarray( + tb_subset.getcol( + "ANTENNA1", + startrow=startrow, + nrow=nrow, + rowincr=1, + ) + ) + antenna2 = np.asarray( + tb_subset.getcol( + "ANTENNA2", + startrow=startrow, + nrow=nrow, + rowincr=1, + ) + ) + ddid = np.asarray( + tb_subset.getcol( + "DATA_DESC_ID", + startrow=startrow, + nrow=nrow, + rowincr=1, + ) + ) + flag_row = np.asarray( + tb_subset.getcol( + "FLAG_ROW", + startrow=startrow, + nrow=nrow, + rowincr=1, + ), + dtype=bool, + ) + + # Sometimes trying to pull out the FLAG column fails. + # In that case, just make a dummy array of all False flags + try: + flags = np.asarray( + tb_subset.getcol( + "FLAG", + startrow=startrow, + nrow=nrow, + rowincr=1, + ), + dtype=bool, + ) + except RuntimeError: + flags = np.zeros((1, 1, nrow), dtype=bool) + + # CASA normally returns UVW as (3, nrow). + if uvw.shape[0] != 3 and uvw.shape[-1] == 3: + uvw = uvw.T + + uvdm = np.hypot(uvw[0], uvw[1]) + + # FLAG normally has dimensions (ncorr, nchan, nrow). + # Keep a row if at least one correlation/channel is unflagged. + flag_axes = tuple(range(flags.ndim - 1)) + completely_flagged = np.all(flags, axis=flag_axes) + + valid = ( + ~flag_row + & ~completely_flagged + & (antenna1 != antenna2) + & np.isfinite(uvdm) + & (uvdm > 0) + & (ddid >= 0) + ) + + # Only take valid values + uv_distance_m.extend(uvdm[valid]) + + tb_subset.close() + + uv_distance_m = np.asarray(uv_distance_m) + + return uv_distance_m + +def get_spw_info( + vis: str, + ignore_wvr: bool = True, +) -> dict: + """Get SPW information from a measurement set. + + Args: + vis (str): Path to measurement set. + ignore_wvr (bool): Whether to ignore WVR SPWs. Defaults to True. + + Returns: + dict: Dictionary containing SPW information, including bandwidth, channel frequencies, + channel width, edge channels, sideband, mean frequency, and number of channels. + """ + + spw_info = {} + + mytb = casaStuff.tbtool() + mytb.open(f"{vis}/SPECTRAL_WINDOW") + + # Keep track of number of rows in the table + nrows = range(mytb.nrows()) + + for i in nrows: + + # Get out useful info + bandwidth = mytb.getcell("TOTAL_BANDWIDTH", i) + chan_freqs = mytb.getcell("CHAN_FREQ", i) + min_freq = min(chan_freqs) + max_freq = max(chan_freqs) + mean_freq = chan_freqs.mean() + num_channels = chan_freqs.shape[0] + chan_width = mytb.getcell("CHAN_WIDTH", i)[0] + net_sideband = mytb.getcell("NET_SIDEBAND", i) + + # Put this all into a dictionary + spw_info[i] = {} + spw_info[i]["bandwidth"] = bandwidth + spw_info[i]["chanFreqs"] = chan_freqs + spw_info[i]["chanWidth"] = chan_width + spw_info[i]["edgeChannels"] = [min_freq, max_freq] + if net_sideband == 2: + spw_info[i]["sideband"] = 1 + else: + spw_info[i]["sideband"] = -1 + spw_info[i]["meanFreq"] = mean_freq + spw_info[i]["numChannels"] = num_channels + + # If ignoring WVR, then remove and log this + if ignore_wvr and (num_channels == 4): + logger.debug(f"Ignoring spectral window {i} because it is WVR related") + spw_info.pop(i) + + mytb.close() + + return spw_info + +def get_scans_for_spw( + vis: str, +) -> dict: + """Get the unique scan numbers for each spectral window in a measurement set. + + Args: + vis (str): Path to measurement set. + + Returns: + dict: Dictionary mapping SPW IDs to unique scan numbers. + """ + + # Get out SPWs for each data description ID + tb = casaStuff.tbtool() + tb.open(f"{vis}/DATA_DESCRIPTION") + spw_for_data_desc_id = tb.getcol("SPECTRAL_WINDOW_ID") + tb.close() + + tb = casaStuff.tbtool() + tb.open(vis) + + data_desc_id = tb.getcol('DATA_DESC_ID') + scans = tb.getcol('SCAN_NUMBER') + + scans_for_spw = {} + for i in spw_for_data_desc_id: + spw = spw_for_data_desc_id[i] + indices = np.where(data_desc_id == i) + scans_for_spw[spw] = np.unique(scans[indices]) + + tb.close() + + return scans_for_spw + +def get_science_spws( + vis: str, + intent: str = "OBSERVE_TARGET#ON_SOURCE", + return_string: bool = True, + return_list_of_strings: bool = False, + return_freq_ranges: bool = False, + tdm: bool = True, + fdm: bool = True, + sqld: bool = False, + chavg: bool = False, +) -> str | list[int | str] | dict: + """ Get science SPWs from a measurement set based. + + Will select SPWs based on specified intent. For ALMA data, + it also can ignore channel-averaged and SQLD SPWs. + + Args: + vis (str): Path to measurement set. + intent (str): The intent to filter by. + Defaults to "OBSERVE_TARGET#ON_SOURCE". + return_string (bool): If True, return a comma-separated string of SPWs. + Defaults to True. + return_list_of_strings (bool): If True, return a list of SPWs as strings. + Defaults to False. + return_freq_ranges (bool): If True, return a dictionary with frequency ranges. + Defaults to False. + tdm (bool): If True, include TDM SPWs. + Defaults to True. + fdm (bool): If True, include FDM SPWs. + Defaults to True. + sqld (bool): If True, include SQLD SPWs. + Defaults to False. + chavg (bool): If True, include channel-averaged SPWs. + Defaults to False. + + Returns: + str: Comma-separated string of science SPWs if return_string is True. + list[int]: List of science SPWs as integers if return_list_of_strings is False. + list[str]: List of science SPWs as strings if return_list_of_strings is True. + dict: Dictionary of science SPWs with frequency ranges if return_freq_ranges is True. + """ + + if return_string and return_list_of_strings: + raise ValueError("You can only specify one of: return_string, return_list_of_strings") + + msmd = casaStuff.msmdtool() + msmd.open(vis) + + all_intents = msmd.intents() + + if intent not in all_intents and intent != "": + for i in all_intents: + if i.find(intent) >= 0: + intent = i + logger.debug(f"Translated intent to {i}") + break + + # Minimum match OBSERVE_TARGET to OBSERVE_TARGET#UNSPECIFIED + value = [i.find(intent.replace("*", "")) for i in all_intents] + + # If any intent gives a match, the mean value of the location list will be > -1 + if np.mean(value) == -1 and intent != "": + logger.warning(f"{intent} not found in this dataset. Available intents: {all_intents}") + + if return_string: + science_spws = "" + else: + science_spws = [] + + else: + # If we don't have an intent, match on wildcards + if intent == "": + intent = "*" + + # Get SPWs for intent, cast back to int + spws = msmd.spwsforintent(intent) + spws = [int(i) for i in spws] + + # Get observatory name, if ALMA data then perform further filtering + observatory_name = get_observatory_name(vis) + + if observatory_name.find("ALMA") >= 0 or observatory_name.find("OSF") >= 0: + + logger.debug( + f"Calling almaspws with:\n chavg={chavg}\n tdm={tdm}\n fdm={fdm}\n sqld={sqld}" + ) + alma_spws = msmd.almaspws( + chavg=chavg, + tdm=tdm, + fdm=fdm, + sqld=sqld, + ) + alma_spws = [int(i) for i in alma_spws] + + if chavg and not sqld: + sqld_spws = msmd.almaspws(sqld=True) + sqld_spws = [int(i) for i in sqld_spws] + + logger.debug(f"removing SQLD ({sqld_spws}) from list") + alma_spws = list(set(alma_spws) - set(sqld_spws)) + if len(spws) == 0 or len(alma_spws) == 0: + science_spws = [] + else: + science_spws = [int(i) for i in np.intersect1d(spws, alma_spws)] + else: + science_spws = spws + + science_spws_dict = {} + for spw in science_spws: + science_spws_dict[spw] = sorted([msmd.chanfreqs(spw)[0], msmd.chanfreqs(spw)[-1]]) + + msmd.close() + + if return_freq_ranges: + science_spws = science_spws_dict + if return_string: + science_spws = ",".join(str(i) for i in science_spws) + elif return_list_of_strings: + science_spws = list([str(i) for i in science_spws]) + else: + science_spws = list(science_spws) + + return science_spws + +def get_observatory_name( + vis: str, +) -> str: + """Get observatory name for given measurement set. + + Args: + vis (str): The path to the measurement set. + + Returns: + str: The observatory name, if found; otherwise, an empty string. + """ + + tb = casaStuff.tbtool() + tb.open(f"{vis}/OBSERVATION") + + try: + observatory_name = tb.getcell("TELESCOPE_NAME") + except RuntimeError: + observatory_name = "" + + tb.close() + + return observatory_name + +def clear_pointing_table( + vis: str, +): + """Removes all rows from the POINTING table of a measurement set. + + Args: + vis (str): Path to the measurement set. + + Returns: + bool: True if the POINTING table was cleared successfully. + """ + + tb = casaStuff.tbtool() + tb.open( + f"{vis}/POINTING", + nomodify=False, + ) + row_numbers = tb.rownumbers() + tb.removerows(row_numbers) + tb.close() + + logger.info(f"Cleared the POINTING table for {vis}") + + return True diff --git a/phangsPipeline/clean_call.py b/phangsPipeline/clean_call.py index 12870db0..23b66c06 100644 --- a/phangsPipeline/clean_call.py +++ b/phangsPipeline/clean_call.py @@ -6,11 +6,10 @@ import os import re -import analysisUtils as au import numpy as np from .casaStuff import imhead -from .casaVisRoutines import estimate_mrs +from .casaVisRoutines import estimate_mrs, estimate_synthesised_beam logger = logging.getLogger(__name__) @@ -276,7 +275,7 @@ def estimate_synthesised_beam( logger.info("Estimating beam from visibilities") vis = self.get_param("vis") - beam = au.estimateSynthesizedBeam(vis) + beam = estimate_synthesised_beam(vis) return beam diff --git a/phangsPipeline/handlerImaging.py b/phangsPipeline/handlerImaging.py index e0084b69..7ba4ebfc 100644 --- a/phangsPipeline/handlerImaging.py +++ b/phangsPipeline/handlerImaging.py @@ -57,9 +57,6 @@ if casa_enabled: - # Analysis utilities - import analysisUtils as au - from .clean_call import CleanCall, CleanCallFunctionDecorator from . import casaImagingRoutines as imr diff --git a/phangsPipeline/handlerImagingChunked.py b/phangsPipeline/handlerImagingChunked.py index 2ab62ee9..183c9dee 100644 --- a/phangsPipeline/handlerImagingChunked.py +++ b/phangsPipeline/handlerImagingChunked.py @@ -29,8 +29,6 @@ if casa_enabled: - import analysisUtils as au - from . import casaImagingRoutines as imr from . import casaMaskingRoutines as msr from . import casaStuff @@ -539,7 +537,7 @@ def return_valid_chunks(self, chunk_num=None): # chunks_iter = self.return_valid_chunks(chunk_num=chunk_num) # # TODO: move to helper function. - # mytb = au.createCasaTool(casaStuff.tbtool) + # mytb = casaStuff.tbtool() # mytb.open(self.full_vis_file, nomodify = True) # colnames = mytb.colnames() # if 'CORRECTED_DATA' in colnames: @@ -590,8 +588,8 @@ def task_split_cube_to_chunks(self, imagename, overwrite=False, chunk_num=None, chunks_iter = self.return_valid_chunks(chunk_num=chunk_num) - myia = au.createCasaTool(casaStuff.iatool) - myrg = au.createCasaTool(casaStuff.rgtool) + myia = casaStuff.iatool() + myrg = casaStuff.rgtool() myia.open(imagename) diff --git a/phangsPipeline/handlerTestImaging.py b/phangsPipeline/handlerTestImaging.py index 4d6ef7b1..0d7c7b63 100644 --- a/phangsPipeline/handlerTestImaging.py +++ b/phangsPipeline/handlerTestImaging.py @@ -55,9 +55,6 @@ if casa_enabled: - # Analysis utilities - import analysisUtils as au - from .clean_call import CleanCall from . import utilsLines as lines @@ -257,7 +254,7 @@ def _get_beam_properties( logger.warning(f"PSF file not found: {psf_file}") return None - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(psf_file) beam_info = {} diff --git a/phangsPipeline/utilsFieldSelection.py b/phangsPipeline/utilsFieldSelection.py index e4b3ee1a..7ac6081e 100644 --- a/phangsPipeline/utilsFieldSelection.py +++ b/phangsPipeline/utilsFieldSelection.py @@ -5,12 +5,12 @@ import re import shutil -import analysisUtils as aU import numpy as np from . import casaStuff +from .casaVisRoutines import get_science_spws -tb = aU.createCasaTool(casaStuff.tbtool) +tb = casaStuff.tbtool() split = casaStuff.split # @@ -108,8 +108,14 @@ def extract_field_selections( # # compute pb if pb_fwhm is None: - spw_dict = aU.getScienceSpws(vis, intent='OBSERVE_TARGET#ON_SOURCE', returnFreqRanges=True) - #print('spw_dict.values()', spw_dict.values()) + spw_dict = get_science_spws( + vis, + intent='OBSERVE_TARGET#ON_SOURCE', + return_freq_ranges=True, + ) + if not isinstance(spw_dict, dict): + raise TypeError('Expecting spw_dict to be a dictionary') + min_freq = np.min(np.array(list(spw_dict.values())).ravel()) pb_fwhm = 1.14*1.22*(3e8/min_freq)*3600*180/(ant_diam*3.1415926) match_radius = pb_factor * pb_fwhm / 2.0 # arcsec diff --git a/phangsPipeline/utilsSingleDish.py b/phangsPipeline/utilsSingleDish.py index b4be1425..5aea6943 100644 --- a/phangsPipeline/utilsSingleDish.py +++ b/phangsPipeline/utilsSingleDish.py @@ -1,5 +1,7 @@ import copy +import logging import os +import re import analysisUtils as au import astropy.units as u @@ -8,19 +10,17 @@ import pylab as pb from analysisUtils import mjdSecondsListToDateTime, mjdsecToUT +from . import casaStuff from .casaStuff import ( casa_version, - tbtool, - msmdtool, - metool, - qatool, - iatool, imhead, importfits, imregrid, imtrans, ) +logger = logging.getLogger(__name__) + casaVersion = "{0}.{1}.{2}".format(*casa_version) def get_first_arr_val(arr): @@ -69,7 +69,7 @@ def prepare_sd_image( os.system('mv -f ' + sd_image_file + '_reorder ' + sd_image_file) # Regrid this to the input measurement set to avoid any weirdness with overlap. - mytb = au.createCasaTool(tbtool) + mytb = casaStuff.tbtool() mytb.open(clean_call.get_param('vis') + '/SPECTRAL_WINDOW') freq = mytb.getcol('CHAN_FREQ') n_chan = mytb.getcol('NUM_CHAN')[0] @@ -108,7 +108,7 @@ def prepare_sd_image( n_chan = cube_info['shape'][-1] n_pol = cube_info['shape'][-2] - myia = au.createCasaTool(iatool) + myia = casaStuff.iatool() myia.open(sd_image_file) restoring_beam = myia.restoringbeam() myia.setrestoringbeam(remove=True) @@ -137,7 +137,7 @@ def get_dish_diameter( freq = hdr["crval4"] * u.Unit(hdr["cunit4"]) # Get restoring beam in radians - myia = au.createCasaTool(iatool) + myia = casaStuff.iatool() myia.open(infile=sdimage) # Beam size is the same for all channels, so just take the first @@ -308,7 +308,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False, if field == '': field = None - mytb = au.createCasaTool(tbtool) + mytb = casaStuff.tbtool() mytb.open(vis) allObsIDs = mytb.getcol('OBSERVATION_ID') obsIDs = np.unique(allObsIDs) @@ -329,7 +329,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False, print("This function requires casa >= 4.1.0.") return - mymsmd = au.createCasaTool(msmdtool) + mymsmd = casaStuff.msmdtool() mymsmd.open(vis) timeranges = {} timecenters = [] @@ -653,7 +653,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False, if (coordSys.find('AZEL') >= 0 and convert): coordSys = 'J2000' print("Converting %d coordinates from AZELGEO to J2000..." % (len(x))) - my_metool = au.createCasaTool(metool) + my_metool = casaStuff.metool() for i in range(len(x)): if ((i+1) % 10000 == 0): print("%d/%d" % (i+1,len(x))) ra, dec = au.computeRADecFromAzElMJD([xrad[i],yrad[i]], mjd=times[i]/86400., @@ -1175,7 +1175,7 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False, sample = i*labelIncrement pb.text(x[sample],y[sample],str(sample),size=8) - myqa = au.createCasaTool(qatool) + myqa = casaStuff.qatool() if (coordSys.find('AZEL') >= 0): azimString = myqa.formxxx('%frad'%(rightAscension), format='deg', prec=5) elevString = myqa.formxxx('%frad'%(declination), format='deg', prec=5) @@ -1252,3 +1252,70 @@ def getTPSampling(vis, obsid=0, showplot=False, plotfile='', debug=False, else: return xSampling, ySampling, largestDimension # end of getTPSampling + +# Get source name +def get_sourcename(filename, + source='all', + ): + + mytb = casaStuff.msmdtool() + mytb.open(filename) + + # Get source field names + fieldnames = mytb.fieldnames() + + # If we're looking for all intents, then just take the (first) science target + if source == 'all': + field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE')[0] + + # Else, key in on the source name. Do this in a regex-style + else: + + # For field selection, this can be done with forward slashes, + # so make sure to strip them + if source.startswith("/"): + source = source.strip("/") + + # Get science fieldnames as names + fieldnames_asname = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE', asnames=True) + + regex_matches = [re.fullmatch(source, x) for x in fieldnames_asname] + regex_idxs = [x is not None for x in regex_matches] + found_field = any(regex_idxs) + + # If we've found a match, update the source name with the first + # fieldname + if found_field: + source = fieldnames_asname[regex_idxs][0] + + if not found_field: + logger.warning(f"source {source} not in field names: {', '.join(fieldnames)}") + return None + + field = mytb.fieldsforname(source) + + # If we don't find anything, get the science targets + if len(field) == 0: + field = mytb.fieldsforintent('OBSERVE_TARGET#ON_SOURCE') + + # Take the first result + field = field[0] + + source_name = fieldnames[field] + + mytb.close() + + return source_name + +# Read source coordinates +def read_source_coordinates(filename,source): + + coord_source = au.getRADecForSource(filename,source) + RA_h = (coord_source.split(' ')[0]).split(':')[0] + RA_m = (coord_source.split(' ')[0]).split(':')[1] + RA_s = (coord_source.split(' ')[0]).split(':')[2] + DEC_d = (coord_source.split(' ')[1]).split(':')[0] + DEC_m = (coord_source.split(' ')[1]).split(':')[1] + DEC_s = (coord_source.split(' ')[1]).split(':')[2] + coord = "J2000 "+str(RA_h)+"h"+str(RA_m)+"m"+str(RA_s[0:6])+" "+str(DEC_d)+"d"+str(DEC_m)+"m"+str(DEC_s) + return coord diff --git a/phangsPipeline/utilsTestImagingPlots.py b/phangsPipeline/utilsTestImagingPlots.py index ba28b4f0..36b694cc 100644 --- a/phangsPipeline/utilsTestImagingPlots.py +++ b/phangsPipeline/utilsTestImagingPlots.py @@ -61,7 +61,6 @@ casa_enabled = is_casa_installed() if casa_enabled: - import analysisUtils as au from . import casaStuff from . import utilsFilenames @@ -107,7 +106,7 @@ def _read_image_data(image_file): # Try CASA image if casa_enabled and os.path.isdir(image_file): try: - myia = au.createCasaTool(casaStuff.iatool) + myia = casaStuff.iatool() myia.open(image_file) data = myia.getchunk() csys = myia.coordsys() diff --git a/tests/test_casa_imports.py b/tests/test_casa_imports.py index 3e2c3361..8dd1cafd 100644 --- a/tests/test_casa_imports.py +++ b/tests/test_casa_imports.py @@ -390,27 +390,20 @@ class TestAnalysisUtilsImports: au_funcs = [ "Atmcal", "MAD", - "ValueMapping", "angularSeparationRadians", - "clearPointingTable", "computeAzElFromRADecMJD", "computeRADecFromAzElMJD", - "createCasaTool", "dateStringToMJD", "getAntennaPadXYZ", - "getBaselineStats", - "getChanWidths", "getObservationStart", "getObservationStartDate", "getOffSourceTimes", "getRADecForField", "getRADecForSource", - "getScienceSpws", "lsrkToTopo", "mjdSecondsListToDateTime", "mjdSecondsToMJDandUT", "parseTimerangeArgument", - "pickCellSize", "radec2rad", "stuffForScienceDataReduction", "timeOnSource", @@ -473,22 +466,6 @@ class TestAnalysisUtilsArguments: "usemstool", ], ), - ( - "getScienceSpws", - [ - "intent", - "returnFreqRanges", - ], - ), - ( - "pickCellSize", - [ - "imsize", - "intent", - "npix", - "pblevel", - ], - ), ] @pytest.mark.casa