diff --git a/.github/workflows/ci-master.yml b/.github/workflows/ci-master.yml new file mode 100644 index 0000000..be2eef1 --- /dev/null +++ b/.github/workflows/ci-master.yml @@ -0,0 +1,54 @@ +name: CI (master) + +on: + push: + branches: [master] + pull_request: + branches: [master] + +jobs: + test: + runs-on: ubuntu-latest + + strategy: + matrix: + python-version: ["3.10", "3.11"] + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Cache pip + uses: actions/cache@v4 + with: + path: ~/.cache/pip + key: ${{ runner.os }}-pip-${{ hashFiles('**/pyproject.toml') }} + restore-keys: | + ${{ runner.os }}-pip- + + - name: Upgrade pip and install build tools + run: python -m pip install --upgrade pip build + + - name: Install package and test deps + run: | + python -m pip install --upgrade pip + pip install -e ".[test]" + # install formatting tools used in CI + python -m pip install black isort + + - name: Black format check + run: python -m black --check . + + - name: Lint + continue-on-error: true + run: flake8 . + + - name: Test + run: | + # Run tests in the tests/ folder explicitly + python -m pytest tests -q diff --git a/pyproject.toml b/pyproject.toml index 86dc90d..874917f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,4 +44,11 @@ spd_scubesml = "splusdata.scubes.scripts:scubesml" where = ["."] # python -m build -# python -m twine upload dist/* \ No newline at end of file +# python -m twine upload dist/* + +[project.optional-dependencies] +test = [ + "pytest>=7.0,<8", + "flake8>=6.0", + "coverage>=7.0" +] \ No newline at end of file diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..6d164d7 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,12 @@ +# requirements.txt gerado a partir de pyproject.toml [project].dependencies +# Use: python -m pip install -r requirements.txt +adss>=1.38 +requests +astropy +astroquery +numpy +pandas +scipy +pillow +pyyaml +tqdm diff --git a/splusdata/connect.py b/splusdata/connect.py index 5ad5c2b..9d750f0 100644 --- a/splusdata/connect.py +++ b/splusdata/connect.py @@ -11,11 +11,12 @@ from astropy.io.votable import from_table, writeto + class connect: - """Class that logs in into splus.cloud and perform all types of download operation. - """ + """Class that logs in into splus.cloud and perform all types of download operation.""" + def __init__(self, username=None, password=None): - """Class may be initialized with splus.cloud user and password. If left in blanck, user and password will be asked on the runtime. + """Class may be initialized with splus.cloud user and password. If left in blanck, user and password will be asked on the runtime. Once logged in you may use all main tools from this class. Examples: @@ -27,64 +28,95 @@ def __init__(self, username=None, password=None): Args: username (str, optional): splus.cloud username password (str, optional): splus.cloud pasword - """ + """ if not username or not password: username = input("splus.cloud username: ") password = getpass("splus.cloud password: ") - - data = {'username':username, 'password':password} - res = requests.post("https://splus.cloud/api/auth/login", data = data) + + data = {"username": username, "password": password} + res = requests.post("https://splus.cloud/api/auth/login", data=data) usr = json.loads(res.content) - self.token = usr['token'] - self.headers = {'Authorization': 'Token ' + self.token} - - res = requests.post("https://splus.cloud/api/auth/collab", headers = self.headers) + self.token = usr["token"] + self.headers = {"Authorization": "Token " + self.token} + + res = requests.post("https://splus.cloud/api/auth/collab", headers=self.headers) collab = json.loads(res.content) - - if collab['collab'] == 'yes': + + if collab["collab"] == "yes": self.collab = True - print('You have access to internal data') + print("You have access to internal data") else: self.collab = False pass - - self.lastres = '' - - - def twelve_band_img(self, ra, dec, radius, noise=0.15, saturation=0.15, R="R,I,F861,Z", G="G,F515,F660", B="U,F378,F395,F410,F430", option=1): - """Function to get twelve band composed images. + + self.lastres = "" + + def twelve_band_img( + self, + ra, + dec, + radius, + noise=0.15, + saturation=0.15, + R="R,I,F861,Z", + G="G,F515,F660", + B="U,F378,F395,F410,F430", + option=1, + ): + """Function to get twelve band composed images. Args: - ra (float): RA coordinate in degrees. + ra (float): RA coordinate in degrees. dec (float): DEC coordinate in degrees. radius (int): Image size in pixels. Final size will be (radius X radius) noise (float, optional): Image noise value. Defaults to 0.15. - saturation (float, optional): Image saturation value. Defaults to 0.15. + saturation (float, optional): Image saturation value. Defaults to 0.15. R (str, optional): Combinations of bands to compose red. Defaults to "R,I,F861,Z". G (str, optional): Combinations of bands to compose green. Defaults to "G,F515,F660". B (str, optional): Combinations of bands to compose blue. Defaults to "U,F378,F395,F410,F430". option (str, optional): in case a coordinate overlap over two fields, use option = 2 if you want the second field. Returns: PIL.Image: Image requested - """ - - - res = requests.get("https://splus.cloud/api/get_image/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + R + "-" + G + "-" + B + "/" + str(noise) + "/" + str(saturation) + "/" + str(option), headers=self.headers) + """ + + res = requests.get( + "https://splus.cloud/api/get_image/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + R + + "-" + + G + + "-" + + B + + "/" + + str(noise) + + "/" + + str(saturation) + + "/" + + str(option), + headers=self.headers, + ) source = json.loads(res.content) - source['filename'] - res = requests.get("https://splus.cloud" + source['filename'], headers=self.headers) + source["filename"] + res = requests.get( + "https://splus.cloud" + source["filename"], headers=self.headers + ) image = Image.open(io.BytesIO(res.content)) - + self.lastcontent = image - self.lastres = '12img' + self.lastres = "12img" return image def get_img(self, ra, dec, radius, R="I", G="R", B="G", stretch=3, Q=8, option=1): - """Function to get three band composed images made by lupton. + """Function to get three band composed images made by lupton. Args: - ra (float): RA coordinate in degrees. + ra (float): RA coordinate in degrees. dec (float): DEC coordinate in degrees. radius (int): Image size in pixels. Final size will be (radius X radius) R (str, optional): Band to compose red. Defaults to "I". @@ -95,22 +127,44 @@ def get_img(self, ra, dec, radius, R="I", G="R", B="G", stretch=3, Q=8, option=1 option (str, optional): in case a coordinate overlap over two fields, use option = 2 if you want the second field. Returns: PIL.Image: Image requested - """ - res = requests.get("https://splus.cloud/api/get_lupton_image/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + str(R) + "/" + str(G) + "/" + str(B) + "/" + str(stretch) + "/" + str(Q) + "/" + str(option), headers=self.headers) + """ + res = requests.get( + "https://splus.cloud/api/get_lupton_image/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + str(R) + + "/" + + str(G) + + "/" + + str(B) + + "/" + + str(stretch) + + "/" + + str(Q) + + "/" + + str(option), + headers=self.headers, + ) source = json.loads(res.content) - source['filename'] - res = requests.get("https://splus.cloud" + source['filename'], headers=self.headers) + source["filename"] + res = requests.get( + "https://splus.cloud" + source["filename"], headers=self.headers + ) image = Image.open(io.BytesIO(res.content)) - + self.lastcontent = image - self.lastres = 'get_cut' + self.lastres = "get_cut" return image - def get_band_img(self, ra, dec, radius, band='R', mode='linear', option=1): + def get_band_img(self, ra, dec, radius, band="R", mode="linear", option=1): """Get image composed with one band. Args: - ra (float): RA coordinate in degrees. + ra (float): RA coordinate in degrees. dec (float): DEC coordinate in degrees. radius (int): Image size in pixels. Final size will be (radius X radius) band (str, optional): Band to compose image. Defaults to 'R'. @@ -118,133 +172,203 @@ def get_band_img(self, ra, dec, radius, band='R', mode='linear', option=1): option (str, optional): in case a coordinate overlap over two fields, use option = 2 if you want the second field. Returns: PIL.Image: Image requested - """ - res = requests.get("https://splus.cloud/api/get_band_image/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + str(band) + "/" + str(mode) + "/" + str(option), headers=self.headers) + """ + res = requests.get( + "https://splus.cloud/api/get_band_image/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + str(band) + + "/" + + str(mode) + + "/" + + str(option), + headers=self.headers, + ) source = json.loads(res.content) - source['filename'] - res = requests.get("https://splus.cloud" + source['filename'], headers=self.headers) + source["filename"] + res = requests.get( + "https://splus.cloud" + source["filename"], headers=self.headers + ) image = Image.open(io.BytesIO(res.content)) - + self.lastcontent = image - self.lastres = 'get_band_image' + self.lastres = "get_band_image" return image - def get_cut(self, ra, dec, radius, band, filepath=None, option=1): """Get fits cut. Args: - ra (float): RA coordinate in degrees. + ra (float): RA coordinate in degrees. dec (float): DEC coordinate in degrees. radius (int): Image size in pixels. Final size will be (radius X radius) band (str): Band requested filepath (str, optional): file path to save result. Defaults to None. option (str, optional): in case a coordinate overlap over two fields, use option = 2 if you want the second field. Returns: - astropy.io.fits: Fits file with image. - """ - if band.upper() == 'ALL': + astropy.io.fits: Fits file with image. + """ + if band.upper() == "ALL": if filepath == None: return 'You must save the file while getting "all" bands' - - elif filepath != None: + + elif filepath != None: filepath = filepath + ".tar.gz" - res = requests.get("https://splus.cloud/api/get_direct_cut/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + "ALL", headers=self.headers) - - with open(filepath, 'wb') as f: + res = requests.get( + "https://splus.cloud/api/get_direct_cut/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + "ALL", + headers=self.headers, + ) + + with open(filepath, "wb") as f: f.write(res.content) f.close() - - return 'File saved to ' + filepath - - res = requests.get("https://splus.cloud/api/get_direct_cut/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + str(band) + "/" + str(option), headers=self.headers) + + return "File saved to " + filepath + + res = requests.get( + "https://splus.cloud/api/get_direct_cut/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + str(band) + + "/" + + str(option), + headers=self.headers, + ) hdu = fits.open(io.BytesIO(res.content)) if filepath != None: hdu.writeto(filepath + ".fz") - + self.lastcontent = hdu - self.lastres = 'cut' + self.lastres = "cut" return hdu def get_cut_weight(self, ra, dec, radius, band, filepath=None, option=1): - """Get weight image fits cut. + """Get weight image fits cut. Args: - ra (float): RA coordinate in degrees. + ra (float): RA coordinate in degrees. dec (float): DEC coordinate in degrees. radius (int): Image size in pixels. Final size will be (radius X radius) band (str): Band requested filepath (str, optional): file path to save result. Defaults to None. option (str, optional): in case a coordinate overlap over two fields, use option = 2 if you want the second field. Returns: - astropy.io.fits: Fits file with image. - """ - if band.upper() == 'ALL': + astropy.io.fits: Fits file with image. + """ + if band.upper() == "ALL": if filepath == None: return 'You must save the file while getting "all" bands' - - elif filepath != None: + + elif filepath != None: filepath = filepath + ".tar.gz" - res = requests.get("https://splus.cloud/api/get_direct_cut_weight/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + "ALL" + "/" + str(option), headers=self.headers) - - with open(filepath, 'wb') as f: + res = requests.get( + "https://splus.cloud/api/get_direct_cut_weight/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + "ALL" + + "/" + + str(option), + headers=self.headers, + ) + + with open(filepath, "wb") as f: f.write(res.content) f.close() - - return 'File saved to ' + filepath - - res = requests.get("https://splus.cloud/api/get_direct_cut_weight/" + str(ra) + "/" + str(dec) + "/" + str(radius) + "/" + str(band) + "/" + str(option), headers=self.headers) + + return "File saved to " + filepath + + res = requests.get( + "https://splus.cloud/api/get_direct_cut_weight/" + + str(ra) + + "/" + + str(dec) + + "/" + + str(radius) + + "/" + + str(band) + + "/" + + str(option), + headers=self.headers, + ) hdu = fits.open(io.BytesIO(res.content)) if filepath != None: hdu.writeto(filepath + ".fz") - + self.lastcontent = hdu - self.lastres = 'cut' + self.lastres = "cut" return hdu - + def get_field(self, field, band): """Get whole 11k field fits. Args: - field (str): field name. + field (str): field name. band (str): Band. Returns: - astropy.io.fits: Fits file with image. - """ - res = requests.get("https://splus.cloud/api/get_direct_field/" + str(field) + "/" + str(band) , headers=self.headers) + astropy.io.fits: Fits file with image. + """ + res = requests.get( + "https://splus.cloud/api/get_direct_field/" + str(field) + "/" + str(band), + headers=self.headers, + ) hdu = fits.open(io.BytesIO(res.content)) - - self.lastres = 'field' + + self.lastres = "field" self.lastcontent = hdu - return hdu + return hdu def get_field_weight(self, field, band): """Get whole 11k weight field fits. Args: - field (str): field name. + field (str): field name. band (str): Band. Returns: - astropy.io.fits: Fits file with image. - """ - res = requests.get("https://splus.cloud/api/get_direct_field_weight/" + str(field) + "/" + str(band) , headers=self.headers) + astropy.io.fits: Fits file with image. + """ + res = requests.get( + "https://splus.cloud/api/get_direct_field_weight/" + + str(field) + + "/" + + str(band), + headers=self.headers, + ) hdu = fits.open(io.BytesIO(res.content)) - - self.lastres = 'field' + + self.lastres = "field" self.lastcontent = hdu - return hdu + return hdu def get_tap_tables(self): - """Get info about available tables. - """ + """Get info about available tables.""" print("Tables and columns available at https://splus.cloud/query/") - + def query(self, query, table_upload=None, publicdata=None): from xml.dom import minidom from astropy.table import Table + """Perform async queries on splus cloud TAP service. Args: @@ -254,30 +378,28 @@ def query(self, query, table_upload=None, publicdata=None): Returns: astropy.table.Table: result table. - """ + """ if self.collab: baselink = "https://splus.cloud/tap/tap/async/" else: baselink = "https://splus.cloud/public-TAP/tap/async/" - if publicdata and self.collab: baselink = "https://splus.cloud/public-TAP/tap/async/" data = { - "request": 'doQuery', - "version": '1.0', - "lang": 'ADQL', - "phase": 'run', + "request": "doQuery", + "version": "1.0", + "lang": "ADQL", + "phase": "run", "query": query, - "format": 'fits' + "format": "fits", } - - + if str(type(table_upload)) != "": - if 'astropy.table' in str(type(table_upload)): + if "astropy.table" in str(type(table_upload)): if len(table_upload) > 6000: - print('Cutting to the first 6000 objects!') + print("Cutting to the first 6000 objects!") table_upload = table_upload[0:6000] table_upload = from_table(table_upload) @@ -293,17 +415,17 @@ def query(self, query, table_upload=None, publicdata=None): IObytes.seek(0) - elif 'astropy.io.votable' in str(type(table_upload)): + elif "astropy.io.votable" in str(type(table_upload)): if table_upload.get_first_table().nrows > 6000: - return 'votable bigger than 6000' + return "votable bigger than 6000" else: IObytes = io.BytesIO() writeto(table_upload, IObytes) IObytes.seek(0) - elif 'DataFrame' in str(type(table_upload)): + elif "DataFrame" in str(type(table_upload)): if len(table_upload) > 6000: - print('Cutting to the first 6000 objects!') + print("Cutting to the first 6000 objects!") table_upload = table_upload[0:6000] table_upload = Table.from_pandas(table_upload) table_upload = from_table(table_upload) @@ -316,78 +438,88 @@ def query(self, query, table_upload=None, publicdata=None): IObytes = io.BytesIO() writeto(table_upload, IObytes) IObytes.seek(0) - else: - return 'Table type not supported' + return "Table type not supported" - data['upload'] = 'upload,param:uplTable' - res = requests.post(baselink , data = data, headers=self.headers, files={'uplTable': IObytes.read()}) + data["upload"] = "upload,param:uplTable" + res = requests.post( + baselink, + data=data, + headers=self.headers, + files={"uplTable": IObytes.read()}, + ) if not table_upload: - res = requests.post(baselink , data = data, headers=self.headers) + res = requests.post(baselink, data=data, headers=self.headers) xmldoc = minidom.parse(io.BytesIO(res.content)) try: - item = xmldoc.getElementsByTagName('phase')[0] + item = xmldoc.getElementsByTagName("phase")[0] process = item.firstChild.data - item = xmldoc.getElementsByTagName('jobId')[0] + item = xmldoc.getElementsByTagName("jobId")[0] jobID = item.firstChild.data - while process == 'EXECUTING': + while process == "EXECUTING": res = requests.get(baselink + jobID, headers=self.headers) xmldoc = minidom.parse(io.BytesIO(res.content)) - item = xmldoc.getElementsByTagName('phase')[0] + item = xmldoc.getElementsByTagName("phase")[0] process = item.firstChild.data time.sleep(5) - if process == 'COMPLETED': - item = xmldoc.getElementsByTagName('result')[0] - link = item.attributes['xlink:href'].value + if process == "COMPLETED": + item = xmldoc.getElementsByTagName("result")[0] + link = item.attributes["xlink:href"].value - link = link.replace("http://192.168.10.23:8080", "https://splus.cloud").replace("http://10.180.0.209:8080", "https://splus.cloud").replace("http://10.180.0.207:8080", "https://splus.cloud").replace("http://10.180.0.219:8080", "https://splus.cloud") + link = ( + link.replace("http://192.168.10.23:8080", "https://splus.cloud") + .replace("http://10.180.0.209:8080", "https://splus.cloud") + .replace("http://10.180.0.207:8080", "https://splus.cloud") + .replace("http://10.180.0.219:8080", "https://splus.cloud") + ) res = requests.get(link, headers=self.headers) - - self.lastres = 'query' + + self.lastres = "query" self.lastcontent = Table.read(io.BytesIO(res.content)) - print('finished') - + print("finished") + return self.lastcontent - if process == 'ERROR': - item = xmldoc.getElementsByTagName('message')[0] + if process == "ERROR": + item = xmldoc.getElementsByTagName("message")[0] message = item.firstChild.data print("Error: ", message) except: - item = xmldoc.getElementsByTagName('INFO') - print(item[0].attributes['value'].value, ": ", item[0].firstChild.data) - + item = xmldoc.getElementsByTagName("INFO") + print(item[0].attributes["value"].value, ": ", item[0].firstChild.data) def checkcoords(self, ra, dec): """Check if coords are in footprint Args: - ra (float): RA coordinate in degrees. + ra (float): RA coordinate in degrees. dec (float): DEC coordinate in degrees. Returns: dict: result. - """ - res = requests.get("https://splus.cloud/api/whichdr/" + str(ra) + "/" + str(dec) , headers=self.headers) + """ + res = requests.get( + "https://splus.cloud/api/whichdr/" + str(ra) + "/" + str(dec), + headers=self.headers, + ) res = res.content.decode("utf-8") - + res = ast.literal_eval(res) return res - def get_last_result(self): """If you missed to save your last query or image into a variable, you may get it here. Returns: - _type_: Last result, may be image or query. - """ + _type_: Last result, may be image or query. + """ return self.lastcontent diff --git a/splusdata/core.py b/splusdata/core.py index bb4b8ac..5887381 100644 --- a/splusdata/core.py +++ b/splusdata/core.py @@ -10,6 +10,7 @@ from splusdata.features.io import print_level from splusdata.features.find_pointings import find_pointing + class SplusdataError(Exception): """Custom exception type for S-PLUS data errors raised by this helper module. @@ -39,6 +40,7 @@ def open_image(image_bytes): If Pillow cannot identify or open the image. """ from PIL import Image + im = Image.open(io.BytesIO(image_bytes)) return im @@ -64,7 +66,7 @@ def save_image(image_bytes, filename): """ im = open_image(image_bytes) im.save(filename) - + # field frame class Core: @@ -84,7 +86,14 @@ class Core: * All methods pass through to a single `adss.ADSSClient` instance. """ - def __init__(self, username=None, password=None, SERVER_IP=f"https://splus.cloud", auto_renew=False, verbose=0): + def __init__( + self, + username=None, + password=None, + SERVER_IP=f"https://splus.cloud", + auto_renew=False, + verbose=0, + ): """Initialize a Core client. Parameters @@ -114,9 +123,9 @@ def __init__(self, username=None, password=None, SERVER_IP=f"https://splus.cloud """ if not username: username = input("splus.cloud username: ") - if not password: + if not password: password = getpass.getpass("splus.cloud password: ") - + self.client = adss.ADSSClient( SERVER_IP, username=username, @@ -124,7 +133,7 @@ def __init__(self, username=None, password=None, SERVER_IP=f"https://splus.cloud ) self.collections = [] self.verbose = verbose - + def _load_collections(self): """Fetch and cache image collections from the server. @@ -149,7 +158,7 @@ def check_available_images_releases(self): Collection names, e.g., ["dr4", "dr5", "dr6", ...]. """ collections = self.client.get_collections() - names = [col['name'] for col in collections] + names = [col["name"] for col in collections] return names def get_collection_id_by_pattern(self, pattern): @@ -172,19 +181,23 @@ def get_collection_id_by_pattern(self, pattern): """ self._load_collections() for col in self.collections: - if pattern in col['name']: + if pattern in col["name"]: return col raise SplusdataError("Collection not found") - + def get_file_metadata(self, field, band, pattern="", data_release="dr4"): collection = self.get_collection_id_by_pattern(data_release) collection_id = collection["id"] - candidates = self.client.list_files(collection_id, filter_str=field, filter_name=band) + candidates = self.client.list_files( + collection_id, filter_str=field, filter_name=band + ) if not candidates and ("-" in field or "_" in field): alt = field.replace("-", "_") if "-" in field else field.replace("_", "-") - candidates = self.client.list_files(collection_id, filter_str=alt, filter_name=band) + candidates = self.client.list_files( + collection_id, filter_str=alt, filter_name=band + ) field = alt if not candidates: @@ -224,14 +237,14 @@ def ok(c): return fz[0] if fz else pool[0] def field_frame( - self, - field, - band, - weight=False, - outfile=None, - data_release="dr4", - timeout = 60, - verbose = False + self, + field, + band, + weight=False, + outfile=None, + data_release="dr4", + timeout=60, + verbose=False, ): """Download and open a full field FITS image. @@ -265,21 +278,31 @@ def field_frame( pattern = "" if verbose: - print(field,band,pattern,data_release) + print(field, band, pattern, data_release) final_candidate = self.get_file_metadata(field, band, pattern, data_release) - + if verbose: print(final_candidate) - + image_bytes = self.client.download_file( - final_candidate['id'], - output_path=outfile, - timeout=timeout + final_candidate["id"], output_path=outfile, timeout=timeout ) - + return fits.open(io.BytesIO(image_bytes)) - - def stamp(self, ra, dec, size, band, weight=False, field_name=None, size_unit="pixels", outfile=None, data_release="dr4", timeout = 60): + + def stamp( + self, + ra, + dec, + size, + band, + weight=False, + field_name=None, + size_unit="pixels", + outfile=None, + data_release="dr4", + timeout=60, + ): """Create and open a FITS stamp (cutout) by coordinates or by object name. Parameters @@ -315,8 +338,8 @@ def stamp(self, ra, dec, size, band, weight=False, field_name=None, size_unit="p If the collection cannot be resolved. """ collection = self.get_collection_id_by_pattern(data_release) - collection_id = collection['id'] - + collection_id = collection["id"] + if weight: weight = "weight" if not field_name: @@ -329,7 +352,7 @@ def stamp(self, ra, dec, size, band, weight=False, field_name=None, size_unit="p size_unit=size_unit, pattern=weight if weight else "", output_path=outfile, - timeout=timeout + timeout=timeout, ) else: stamp_bytes = self.client.stamp_images.create_stamp_by_object( @@ -342,12 +365,26 @@ def stamp(self, ra, dec, size, band, weight=False, field_name=None, size_unit="p size_unit=size_unit, pattern=weight if weight else "", output_path=outfile, - timeout=timeout + timeout=timeout, ) - + return fits.open(io.BytesIO(stamp_bytes)) - def lupton_rgb(self, ra, dec, size, R="I", G="R", B="G", Q=8, stretch=3, field_name=None, size_unit="pixels", outfile=None, data_release="dr4"): + def lupton_rgb( + self, + ra, + dec, + size, + R="I", + G="R", + B="G", + Q=8, + stretch=3, + field_name=None, + size_unit="pixels", + outfile=None, + data_release="dr4", + ): """Create a Lupton RGB composite and return a PIL image. Parameters @@ -377,7 +414,7 @@ def lupton_rgb(self, ra, dec, size, R="I", G="R", B="G", Q=8, stretch=3, field_n Composite RGB image. """ collection = self.get_collection_id_by_pattern(data_release) - collection_id = collection['id'] + collection_id = collection["id"] if not field_name: stamp_bytes = self.client.create_rgb_image_by_coordinates( @@ -391,7 +428,7 @@ def lupton_rgb(self, ra, dec, size, R="I", G="R", B="G", Q=8, stretch=3, field_n Q=Q, size_unit=size_unit, stretch=stretch, - output_path=outfile + output_path=outfile, ) else: stamp_bytes = self.client.lupton_images.create_rgb_by_object( @@ -406,12 +443,27 @@ def lupton_rgb(self, ra, dec, size, R="I", G="R", B="G", Q=8, stretch=3, field_n Q=Q, size_unit=size_unit, stretch=stretch, - output_path=outfile + output_path=outfile, ) return Image.open(io.BytesIO(stamp_bytes)) - def trilogy_image(self, ra, dec, size, R=["R", "I", "F861", "Z"], G=["G", "F515", "F660"], B=["U", "F378", "F395", "F410", "F430"], noiselum=0.15, satpercent=0.15, colorsatfac=2, size_unit="pixels", field_name=None, outfile=None, data_release="dr4"): + def trilogy_image( + self, + ra, + dec, + size, + R=["R", "I", "F861", "Z"], + G=["G", "F515", "F660"], + B=["U", "F378", "F395", "F410", "F430"], + noiselum=0.15, + satpercent=0.15, + colorsatfac=2, + size_unit="pixels", + field_name=None, + outfile=None, + data_release="dr4", + ): """Create a Trilogy RGB composite (multi-filter blend) and return a PIL image. Parameters @@ -443,7 +495,7 @@ def trilogy_image(self, ra, dec, size, R=["R", "I", "F861", "Z"], G=["G", "F515" Composite RGB image (Trilogy method). """ collection = self.get_collection_id_by_pattern(data_release) - collection_id = collection['id'] + collection_id = collection["id"] if not field_name: stamp_bytes = self.client.trilogy_images.create_trilogy_rgb_by_coordinates( @@ -458,7 +510,7 @@ def trilogy_image(self, ra, dec, size, R=["R", "I", "F861", "Z"], G=["G", "F515" noiselum=noiselum, satpercent=satpercent, colorsatfac=colorsatfac, - output_path=outfile + output_path=outfile, ) else: stamp_bytes = self.client.trilogy_images.create_trilogy_rgb_by_object( @@ -474,12 +526,21 @@ def trilogy_image(self, ra, dec, size, R=["R", "I", "F861", "Z"], G=["G", "F515" size_unit=size_unit, satpercent=satpercent, colorsatfac=colorsatfac, - output_path=outfile + output_path=outfile, ) return Image.open(io.BytesIO(stamp_bytes)) - - def query(self, query, table_upload=None, table_name=None, verbose = False, timeout = 3200, execution_mode = "async", lang = "astroql"): + + def query( + self, + query, + table_upload=None, + table_name=None, + verbose=False, + timeout=3200, + execution_mode="async", + lang="astroql", + ): """Execute a server-side query; optionally upload a small table first. Parameters @@ -511,36 +572,40 @@ def query(self, query, table_upload=None, table_name=None, verbose = False, time if table_upload is not None and table_name is not None: import pandas as pd from astropy.table import Table - + table_upload_bytes = None if isinstance(table_upload, pd.DataFrame): table_upload_bytes = table_upload.to_csv(index=False).encode() elif isinstance(table_upload, Table): - table_upload_bytes = table_upload.to_pandas().to_csv(index=False).encode() + table_upload_bytes = ( + table_upload.to_pandas().to_csv(index=False).encode() + ) else: - raise ValueError("table_upload must be a pandas DataFrame or an astropy Table") + raise ValueError( + "table_upload must be a pandas DataFrame or an astropy Table" + ) if execution_mode == "async": response = self.client.query_and_wait( query_text=query, table_name=table_name, - file=table_upload_bytes, + file=table_upload_bytes, verbose=verbose, - timeout = timeout, - mode = lang, + timeout=timeout, + mode=lang, ) else: response = self.client.query( query_text=query, table_name=table_name, - file=table_upload_bytes, - timeout = 10, - mode = lang, + file=table_upload_bytes, + timeout=10, + mode=lang, ) - + return response.data - def get_zp_file(self, field, band, data_release = "dr6"): + def get_zp_file(self, field, band, data_release="dr6"): """Download and parse the per-field zero-point model (DR6). Parameters @@ -565,22 +630,25 @@ def get_zp_file(self, field, band, data_release = "dr6"): If the downloaded bytes are not valid JSON. """ import json + collection = self.get_collection_id_by_pattern(data_release) - collection_id = collection['id'] - + collection_id = collection["id"] + files = self.client.list_files( - collection_id, - filter_str=f"{field}_{band}_zp", + collection_id, + filter_str=f"{field}_{band}_zp", ) if len(files) == 0: - raise SplusdataError(f"No zp model found for field {field} in band {band} in {data_release}") + raise SplusdataError( + f"No zp model found for field {field} in band {band} in {data_release}" + ) file = files[0] - + print_level(f"Downloading zp_model {file['filename']}", 1, self.verbose) - json_bytes = self.client.download_file(file["id"], timeout = 20) + json_bytes = self.client.download_file(file["id"], timeout=20) json_data = json.loads(json_bytes) return json_data - + def get_zp(self, field, band, ra, dec): """Evaluate the local zero point at a sky position using the field model. @@ -606,11 +674,23 @@ def get_zp(self, field, band, ra, dec): Any error propagated from `zp_at_coord` evaluation. """ model = self.get_zp_file(field, band) - + from splusdata.features.zeropoints.zp_map import zp_at_coord + return zp_at_coord(model, ra, dec) - - def calibrated_stamp(self, ra, dec, size, band, weight=False, field_name=None, size_unit="pixels", outfile=None, data_release="dr6"): + + def calibrated_stamp( + self, + ra, + dec, + size, + band, + weight=False, + field_name=None, + size_unit="pixels", + outfile=None, + data_release="dr6", + ): """Create a stamp and return a photometrically calibrated PrimaryHDU. This computes a cutout via `stamp(...)`, then loads the appropriate DR6+ @@ -650,23 +730,39 @@ def calibrated_stamp(self, ra, dec, size, band, weight=False, field_name=None, s Exception Propagates any calibration errors from `calibrate_hdu_with_zpmodel`. """ - stamp = self.stamp(ra, dec, size, band, weight=weight, field_name=field_name, size_unit=size_unit, data_release=data_release) - - if not weight: - from splusdata.features.zeropoints.zp_image import calibrate_hdu_with_zpmodel - zp_model = self.get_zp_file(stamp[1].header["FIELD"], stamp[1].header["FILTER"], data_release=data_release) - + stamp = self.stamp( + ra, + dec, + size, + band, + weight=weight, + field_name=field_name, + size_unit=size_unit, + data_release=data_release, + ) + + if not weight: + from splusdata.features.zeropoints.zp_image import ( + calibrate_hdu_with_zpmodel, + ) + + zp_model = self.get_zp_file( + stamp[1].header["FIELD"], + stamp[1].header["FILTER"], + data_release=data_release, + ) + calibrated_hdu, factor_map = calibrate_hdu_with_zpmodel( stamp[1], zp_model, in_place=False, return_factor=True ) stamp[1] = calibrated_hdu stamp.append(fits.ImageHDU(factor_map, name="ZP_FACTOR")) - + if outfile: stamp.writeto(outfile, overwrite=True) return stamp - + def get_zps_field(self, ras, decs, field, band, data_release="dr6"): """ Compute zero-point (ZP) values for a set of coordinates in a specific field and band. @@ -690,10 +786,11 @@ def get_zps_field(self, ras, decs, field, band, data_release="dr6"): Array of zero-point values (in magnitudes) for each input (RA, Dec) pair. """ zp_model = self.get_zp_file(field, band, data_release=data_release) - + from splusdata.features.zeropoints.zp_image import compute_zp_for_coords_array + return compute_zp_for_coords_array(ras, decs, zp_model) - + def check_coords(self, ra, dec, radius=1 * u.degree): """Check which DR contains a pointing within `radius` of (ra, dec). @@ -715,22 +812,21 @@ def check_coords(self, ra, dec, radius=1 * u.degree): Exception Propagates any errors from `find_pointing`. """ - + return find_pointing(ra, dec, radius=radius) - + def check_coords_query(self, ra, dec): res = self.query( f"SELECT top 10 field from idr6.idr6 where cone(ra, dec, {ra},{dec}, 0.01)", - execution_mode = "sync" + execution_mode="sync", ) - + # res is a df with a column 'field' - fields = res['field'].tolist() + fields = res["field"].tolist() # get unique fields fields = list(set(fields)) return fields - def download_collection( self, collection, @@ -738,7 +834,7 @@ def download_collection( path_key="full_path", root_marker="splus", max_workers=2, - **kwargs + **kwargs, ): import os from concurrent.futures import ThreadPoolExecutor, as_completed @@ -780,9 +876,7 @@ def _download_one(f): print(f"Downloading {rel_path}") self.client.download_file( - file_id=f["id"], - output_path=final_path, - timeout=180 + file_id=f["id"], output_path=final_path, timeout=180 ) return f"Downloaded {rel_path}" @@ -797,10 +891,7 @@ def _download_one(f): skip = 0 while True: files = self.client.list_files( - collection_id=collection["id"], - skip=skip, - limit=200, - **kwargs + collection_id=collection["id"], skip=skip, limit=200, **kwargs ) if len(files) == 0: @@ -812,4 +903,4 @@ def _download_one(f): for future in as_completed(futures): print(future.result()) - skip += 200 \ No newline at end of file + skip += 200 diff --git a/splusdata/features/extinction.py b/splusdata/features/extinction.py index c148500..848a2e2 100644 --- a/splusdata/features/extinction.py +++ b/splusdata/features/extinction.py @@ -1,7 +1,7 @@ # Import packages and configure the folder to save the dust maps from astropy.coordinates import SkyCoord import pandas as pd -import numpy as np +import numpy as np import os @@ -12,9 +12,10 @@ except: pass + class SplusExtinction: """ - A class to apply extinction correction to astronomical magnitudes based on the + A class to apply extinction correction to astronomical magnitudes based on the Cardelli, Clayton, & Mathis (1989) law using CSFD dust maps. Parameters @@ -22,10 +23,10 @@ class SplusExtinction: data_dir : str, optional Directory to save the downloaded dust maps, by default "../data/outputs/dustmaps". bands : list of str, optional - List of band names for which extinction correction will be applied, + List of band names for which extinction correction will be applied, by default ['u', 'J0378', 'J0395', 'J0410', 'J0430', 'g', 'J0515', 'r', 'J0660', 'i', 'J0861', 'z']. wavelengths : list of int, optional - List of wavelengths (in Angstroms) corresponding to the bands, by default + List of wavelengths (in Angstroms) corresponding to the bands, by default [3536, 3770, 3940, 4094, 4292, 4751, 5133, 6258, 6614, 7690, 8611, 8831]. Attributes @@ -45,16 +46,16 @@ class SplusExtinction: ------- prepare_map() Prepares the extinction map for querying by downloading and configuring the CSFD map. - + get_extinction(df) Calculates the extinction for given coordinates and adds extinction values to the DataFrame. - + apply(df, mag_cols) Applies extinction correction to the magnitude columns in the DataFrame. ```python import SplusExtinction - + df = pd.DataFrame({ 'ra': [0, 1, 2], 'dec': [0, 1, 2], @@ -63,44 +64,70 @@ class SplusExtinction: ... #12 bands of magnitudes }) - + ext = SplusExtinction("dustmaps_dir/") ext.prepare_map() - - df = ext.get_extinction(df) + + df = ext.get_extinction(df) df = ext.apply(df, [band + "_auto" for band in ext.bands]) ``` """ - + def __init__( - self, - dustmap_dir = "../data/outputs/dustmaps", - bands = ['u', 'J0378', 'J0395', 'J0410', 'J0430', 'g', 'J0515', 'r', 'J0660', 'i', 'J0861', 'z'], - wavelengths = [3536, 3770, 3940, 4094, 4292, 4751, 5133, 6258, 6614, 7690, 8611, 8831] + self, + dustmap_dir="../data/outputs/dustmaps", + bands=[ + "u", + "J0378", + "J0395", + "J0410", + "J0430", + "g", + "J0515", + "r", + "J0660", + "i", + "J0861", + "z", + ], + wavelengths=[ + 3536, + 3770, + 3940, + 4094, + 4292, + 4751, + 5133, + 6258, + 6614, + 7690, + 8611, + 8831, + ], ): - self.data_dir = dustmap_dir + self.data_dir = dustmap_dir self.bands = bands self.wavelengths = wavelengths self.extinction_columns = [f"ext_{band}" for band in self.bands] - + self.extinction_map = None - + if not extinction: raise ImportError("The 'extinction' package is required for this class.") if not dustmaps: raise ImportError("The 'dustmaps' package is required for this class.") - + def prepare_map(self): """Prepares the extinction map by downloading the CSFD map and configuring the file paths.""" - config['data_dir'] = self.data_dir # Folder to save the dust maps - dustmaps.csfd.fetch() # Downloads the dust map if it is not already downloaded + config["data_dir"] = self.data_dir # Folder to save the dust maps + dustmaps.csfd.fetch() # Downloads the dust map if it is not already downloaded self.extinction_map = dustmaps.csfd.CSFDQuery( - map_fname=os.path.join(self.data_dir, 'csfd/csfd_ebv.fits'), - mask_fname=os.path.join(self.data_dir, 'csfd/mask.fits') + map_fname=os.path.join(self.data_dir, "csfd/csfd_ebv.fits"), + mask_fname=os.path.join(self.data_dir, "csfd/mask.fits"), ) - - def get_extinction(self, df, ra_col='ra', dec_col='dec'): + + def get_extinction(self, df, ra_col="ra", dec_col="dec"): """ Calculates extinction values based on RA and DEC coordinates using the CCM89 Law. @@ -115,9 +142,9 @@ def get_extinction(self, df, ra_col='ra', dec_col='dec'): DataFrame with added columns for extinction values in each band. """ # Obtaining E(B-V) and Av in a given RA, DEC position - input_file_coords = SkyCoord(df[ra_col], df[dec_col], frame='icrs', unit='deg') + input_file_coords = SkyCoord(df[ra_col], df[dec_col], frame="icrs", unit="deg") ebv = self.extinction_map(input_file_coords) - av = 3.1*ebv + av = 3.1 * ebv # Calculating the extinction on the S-PLUS bands using the Cardelli, Clayton & Mathis law. lambdas = np.array(self.wavelengths).astype(float) @@ -125,13 +152,13 @@ def get_extinction(self, df, ra_col='ra', dec_col='dec'): extinctions = [] for i in range(len(av)): extinctions.append(extinction.ccm89(lambdas, av[i], 3.1)) - + extinction_df = pd.DataFrame(extinctions, columns=self.extinction_columns) - + df = pd.concat([df, extinction_df], axis=1) - + return df - + def apply(self, df, mag_cols): """ Applies extinction correction to the magnitude columns in the DataFrame. @@ -147,17 +174,19 @@ def apply(self, df, mag_cols): ------- pandas.DataFrame DataFrame with new columns for extinction-corrected magnitudes (_extcorr cols). - + Raises ------ ValueError If the number of columns in `mag_cols` and extinction columns do not match. """ if len(mag_cols) != len(self.extinction_columns): - raise ValueError("The number of columns in mag_cols and ext_cols must be the same.") - + raise ValueError( + "The number of columns in mag_cols and ext_cols must be the same." + ) + # Correct the magnitudes for extinction for mag_col, ext_col in zip(mag_cols, self.extinction_columns): df[mag_col + "_extcorr"] = df[mag_col] - df[ext_col] - - return df \ No newline at end of file + + return df diff --git a/splusdata/features/filterbw.py b/splusdata/features/filterbw.py index 06db613..746f6f9 100644 --- a/splusdata/features/filterbw.py +++ b/splusdata/features/filterbw.py @@ -7,232 +7,244 @@ from scipy import interpolate from scipy import fft + ##Butterworth def filter_bw(hdu): - #frequencia de corte da filtragem a ser utilizada. Para os dados do SPLUS 0.5 é um valor ótimo. - cutoff_frequency_x=0.5 - cutoff_frequency_y=0.5 + # frequencia de corte da filtragem a ser utilizada. Para os dados do SPLUS 0.5 é um valor ótimo. + cutoff_frequency_x = 0.5 + cutoff_frequency_y = 0.5 - #ordem da filtragem, para os dados do SPLUS a melhor ordem é 1. Se usar 2, começam a surgir anéis. - order=1.0 + # ordem da filtragem, para os dados do SPLUS a melhor ordem é 1. Se usar 2, começam a surgir anéis. + order = 1.0 - #formato do filtro, não mudar. Para SPLUS elipse2= yes - elipse='no' - elipse2='yes' - elipsevssquare='no' + # formato do filtro, não mudar. Para SPLUS elipse2= yes + elipse = "no" + elipse2 = "yes" + elipsevssquare = "no" - #tamanho dos pixels - CDELT1=0.55 - CDELT2=0.55 + # tamanho dos pixels + CDELT1 = 0.55 + CDELT2 = 0.55 try: - data=hdu[0].data + data = hdu[0].data if data == None: - data=hdu[1].data + data = hdu[1].data except: - data=hdu[1].data - - #y-eixo vertical da imagem - #-eixo horizontal da imagem - - y,x =data.shape - - - #imagem zerada com tamanho maior do que o original, ela recebe a informação da imagem original - imagemzero=numpy.zeros((y+6,x+6)) - - for i in range(0,x): - for j in range(0,y): - imagemzero[3+j,3+i]=data[j,i] - - #replicação da informação da imagem nas bordas(borda esquerda do cubo) - - for i in range(0,3): - for j in range(0,y): - imagemzero[3+j,i]=data[j,0] - - #replicação da informação da imagem nas bordas(borda direita do cubo) - - for i in range(0,3): - for j in range(0,y): - imagemzero[3+j,3+x+i]=data[j,x-1] + data = hdu[1].data - #replicação da informação da imagem nas bordas(borda de baixo do cubo) + # y-eixo vertical da imagem + # -eixo horizontal da imagem + + y, x = data.shape - for i in range(0,x): - for j in range(0,3): - imagemzero[j,3+i]=data[0,i] - - #replicação da informação da imagem nas bordas(borda de cima do cubo) - - for i in range(0,x): - for j in range(0,3): - imagemzero[3+y+j,3+i]= data[y-1,i] - - #sobrou os 4 cantos da imagem com zeros, precisamos preenche-los. - - #canto inferior esquerdo - for i in range(0,3): - for j in range(0,3): - imagemzero[j,i]=imagemzero[j,3] - - #canto inferior direito - for i in range(0,3): - for j in range(0,3): - imagemzero[j,3+i+x]=imagemzero[j,3+x-1] - - #canto superior esquerdo - for i in range (0,3): - for j in range(0,3): - imagemzero[3+y+j,i]=imagemzero[3+y+j,3] - - #canto superior direito - for i in range(0,3): - for j in range(0,3): - imagemzero[3+y+j,3+x+i]=imagemzero[3+y+j,3+x-1] - - y1,x1=imagemzero.shape - - #cria uma imagem maior com bordas zeradas mas com a imagemzero no meio (segunda ampliação do cubo) - imagemfinalexpandida=numpy.zeros((y1+30,x1+30)) - imagemfinalexpandida[15:y1+15,15:x1+15]=imagemzero - - interpol_values=numpy.zeros((2,2)) - interpol_values1=numpy.zeros((2,2)) - interpol_values2=numpy.zeros((2,2)) - interpol_values3=numpy.zeros((2,2)) - - temporary_column=numpy.zeros((y1+2,2)) - temporary_column[y1+1,0]=y1+30 - temporary_column[1:y1+1,0]=range(15,y1+15) - - #as bordas desta imagem devem ter valores caindo para zero. - - interpol_values[1,0]=15 - for i in range(15,x1+30): - interpol_values[1,1]=imagemfinalexpandida[15,i] - interpolate_function=interpolate.interp1d(interpol_values[0:2,0],interpol_values[0:2,1], kind='linear') - imagemfinalexpandida[0:16,i]=interpolate_function(numpy.arange(16)) - - interpol_values1[0,0]=y1+15-1 - interpol_values1[1,0]=y1+30 - interpol_values1[1,1]=0.0 - for i in range(15,x1+30): - interpol_values1[0,1]=imagemfinalexpandida[y1+15-1,i] - interpolate_function1=interpolate.interp1d(interpol_values1[0:2,0],interpol_values1[0:2,1], kind='linear') - imagemfinalexpandida[y1+15-1:y1+30,i]=interpolate_function1(y1+15-1+numpy.arange(16)) - - interpol_values2[0,0]=0.0 - interpol_values2[1,0]=15 - interpol_values2[0,1]=0.0 - for j in range(0,y1+30): - interpol_values2[1,1]=imagemfinalexpandida[j,15] - interpolate_function2=interpolate.interp1d(interpol_values2[0:2,0],interpol_values2[0:2,1], kind='linear') - imagemfinalexpandida[j,0:16]=interpolate_function2(numpy.arange(16)) - - interpol_values3[0,0]=x1+15-1 - interpol_values3[1,0]=x1+30 - interpol_values3[1,1]=0.0 - for j in range(0, y1+30): - interpol_values3[0,1]=imagemfinalexpandida[j,x1+15-1] - interpolate_function3=interpolate.interp1d(interpol_values3[0:2,0],interpol_values3[0:2,1], kind='linear') - imagemfinalexpandida[j,x1+15-1:x1+30]=interpolate_function3(x1+15-1+numpy.arange(16)) - - yf,xf= imagemfinalexpandida.shape - - #criar o filtro com o dobro do tamanho do imagemfinalexpandida (que vai ser o tamanho da transformada de fourier) - filtro=numpy.zeros((2*yf,2*xf)) - - aa=cutoff_frequency_x*xf - bb=cutoff_frequency_y*yf - cc=cutoff_frequency_x*xf - dd=cutoff_frequency_y*yf - - exporder=2.0*order - - imfilter1=filtro - xm=xf - ym=yf - - for i in range(0,2*xf): - for j in range(0,2*yf): - conta= math.sqrt(((i-xm)/aa)**2+((j-ym)/bb)**2) - imfilter1[j,i]=1.0/(1.0+conta**exporder) - - imfilter2=filtro - for i in range(0,2*xf): - for j in range(0,2*yf): - conta1=abs(i-xm)/cc - conta2=abs(j-ym)/dd - imfilter2[j,i]=(1.0/(1.0+conta1**exporder))*(1.0/(1.0+conta2**exporder)) - - #define o filtro final (baseando-se no seu formato) - if elipse == 'yes': - filtro=imfilter1 - elif elipse2 == 'yes': - filtro=imfilter1*imfilter1 - elif elipsevssquare == 'yes': - filtro=imfilter1*imfilter2 - - #imagem que vai receber a filtragem - imagemsaida=imagemfinalexpandida - npadd=numpy.zeros((2*yf,2*xf)) - - - #faz a filtragem - - imagein=imagemfinalexpandida[0:yf,0:xf] - imageflag=npadd - #padding - procedimento de multiplicar certos pixels da imagem por -1, para que a transformada de Fourier fique com - #frequência = 0 no centro. - #faz o padding para cada imagem - for i in range(0,2*xf): - for j in range(0,2*yf): - if (i<=xf-1) and (j > yf-1): - if((i+j)%2) != 0: - flag=1 + # imagem zerada com tamanho maior do que o original, ela recebe a informação da imagem original + imagemzero = numpy.zeros((y + 6, x + 6)) + + for i in range(0, x): + for j in range(0, y): + imagemzero[3 + j, 3 + i] = data[j, i] + + # replicação da informação da imagem nas bordas(borda esquerda do cubo) + + for i in range(0, 3): + for j in range(0, y): + imagemzero[3 + j, i] = data[j, 0] + + # replicação da informação da imagem nas bordas(borda direita do cubo) + + for i in range(0, 3): + for j in range(0, y): + imagemzero[3 + j, 3 + x + i] = data[j, x - 1] + + # replicação da informação da imagem nas bordas(borda de baixo do cubo) + + for i in range(0, x): + for j in range(0, 3): + imagemzero[j, 3 + i] = data[0, i] + + # replicação da informação da imagem nas bordas(borda de cima do cubo) + + for i in range(0, x): + for j in range(0, 3): + imagemzero[3 + y + j, 3 + i] = data[y - 1, i] + + # sobrou os 4 cantos da imagem com zeros, precisamos preenche-los. + + # canto inferior esquerdo + for i in range(0, 3): + for j in range(0, 3): + imagemzero[j, i] = imagemzero[j, 3] + + # canto inferior direito + for i in range(0, 3): + for j in range(0, 3): + imagemzero[j, 3 + i + x] = imagemzero[j, 3 + x - 1] + + # canto superior esquerdo + for i in range(0, 3): + for j in range(0, 3): + imagemzero[3 + y + j, i] = imagemzero[3 + y + j, 3] + + # canto superior direito + for i in range(0, 3): + for j in range(0, 3): + imagemzero[3 + y + j, 3 + x + i] = imagemzero[3 + y + j, 3 + x - 1] + + y1, x1 = imagemzero.shape + + # cria uma imagem maior com bordas zeradas mas com a imagemzero no meio (segunda ampliação do cubo) + imagemfinalexpandida = numpy.zeros((y1 + 30, x1 + 30)) + imagemfinalexpandida[15 : y1 + 15, 15 : x1 + 15] = imagemzero + + interpol_values = numpy.zeros((2, 2)) + interpol_values1 = numpy.zeros((2, 2)) + interpol_values2 = numpy.zeros((2, 2)) + interpol_values3 = numpy.zeros((2, 2)) + + temporary_column = numpy.zeros((y1 + 2, 2)) + temporary_column[y1 + 1, 0] = y1 + 30 + temporary_column[1 : y1 + 1, 0] = range(15, y1 + 15) + + # as bordas desta imagem devem ter valores caindo para zero. + + interpol_values[1, 0] = 15 + for i in range(15, x1 + 30): + interpol_values[1, 1] = imagemfinalexpandida[15, i] + interpolate_function = interpolate.interp1d( + interpol_values[0:2, 0], interpol_values[0:2, 1], kind="linear" + ) + imagemfinalexpandida[0:16, i] = interpolate_function(numpy.arange(16)) + + interpol_values1[0, 0] = y1 + 15 - 1 + interpol_values1[1, 0] = y1 + 30 + interpol_values1[1, 1] = 0.0 + for i in range(15, x1 + 30): + interpol_values1[0, 1] = imagemfinalexpandida[y1 + 15 - 1, i] + interpolate_function1 = interpolate.interp1d( + interpol_values1[0:2, 0], interpol_values1[0:2, 1], kind="linear" + ) + imagemfinalexpandida[y1 + 15 - 1 : y1 + 30, i] = interpolate_function1( + y1 + 15 - 1 + numpy.arange(16) + ) + + interpol_values2[0, 0] = 0.0 + interpol_values2[1, 0] = 15 + interpol_values2[0, 1] = 0.0 + for j in range(0, y1 + 30): + interpol_values2[1, 1] = imagemfinalexpandida[j, 15] + interpolate_function2 = interpolate.interp1d( + interpol_values2[0:2, 0], interpol_values2[0:2, 1], kind="linear" + ) + imagemfinalexpandida[j, 0:16] = interpolate_function2(numpy.arange(16)) + + interpol_values3[0, 0] = x1 + 15 - 1 + interpol_values3[1, 0] = x1 + 30 + interpol_values3[1, 1] = 0.0 + for j in range(0, y1 + 30): + interpol_values3[0, 1] = imagemfinalexpandida[j, x1 + 15 - 1] + interpolate_function3 = interpolate.interp1d( + interpol_values3[0:2, 0], interpol_values3[0:2, 1], kind="linear" + ) + imagemfinalexpandida[j, x1 + 15 - 1 : x1 + 30] = interpolate_function3( + x1 + 15 - 1 + numpy.arange(16) + ) + + yf, xf = imagemfinalexpandida.shape + + # criar o filtro com o dobro do tamanho do imagemfinalexpandida (que vai ser o tamanho da transformada de fourier) + filtro = numpy.zeros((2 * yf, 2 * xf)) + + aa = cutoff_frequency_x * xf + bb = cutoff_frequency_y * yf + cc = cutoff_frequency_x * xf + dd = cutoff_frequency_y * yf + + exporder = 2.0 * order + + imfilter1 = filtro + xm = xf + ym = yf + + for i in range(0, 2 * xf): + for j in range(0, 2 * yf): + conta = math.sqrt(((i - xm) / aa) ** 2 + ((j - ym) / bb) ** 2) + imfilter1[j, i] = 1.0 / (1.0 + conta**exporder) + + imfilter2 = filtro + for i in range(0, 2 * xf): + for j in range(0, 2 * yf): + conta1 = abs(i - xm) / cc + conta2 = abs(j - ym) / dd + imfilter2[j, i] = (1.0 / (1.0 + conta1**exporder)) * ( + 1.0 / (1.0 + conta2**exporder) + ) + + # define o filtro final (baseando-se no seu formato) + if elipse == "yes": + filtro = imfilter1 + elif elipse2 == "yes": + filtro = imfilter1 * imfilter1 + elif elipsevssquare == "yes": + filtro = imfilter1 * imfilter2 + + # imagem que vai receber a filtragem + imagemsaida = imagemfinalexpandida + npadd = numpy.zeros((2 * yf, 2 * xf)) + + # faz a filtragem + + imagein = imagemfinalexpandida[0:yf, 0:xf] + imageflag = npadd + # padding - procedimento de multiplicar certos pixels da imagem por -1, para que a transformada de Fourier fique com + # frequência = 0 no centro. + # faz o padding para cada imagem + for i in range(0, 2 * xf): + for j in range(0, 2 * yf): + if (i <= xf - 1) and (j > yf - 1): + if ((i + j) % 2) != 0: + flag = 1 else: - flag=-1 - imageflag[j,i]=flag*imagein[j-yf,i] + flag = -1 + imageflag[j, i] = flag * imagein[j - yf, i] else: - imageflag[j,i]=0.0 - - #faz a transformada de fourier - twodfft=fft.fft2(imageflag) - #multiplica pelo filtro - filtragem=twodfft*filtro - #faz a transformada de fourier inversa - twodifft=fft.ifft2(filtragem) - invfft=numpy.real(twodifft) - paddout=npadd - imageout=imagein - - #remove o padding - for i in range(0,2*xf): - for j in range(0,2*yf): - if((i+j)%2) != 0: - flag=1 + imageflag[j, i] = 0.0 + + # faz a transformada de fourier + twodfft = fft.fft2(imageflag) + # multiplica pelo filtro + filtragem = twodfft * filtro + # faz a transformada de fourier inversa + twodifft = fft.ifft2(filtragem) + invfft = numpy.real(twodifft) + paddout = npadd + imageout = imagein + + # remove o padding + for i in range(0, 2 * xf): + for j in range(0, 2 * yf): + if ((i + j) % 2) != 0: + flag = 1 else: - flag=-1 - paddout[j,i]=flag*invfft[j,i] - if (i<=xf-1) and (j >yf-1): - imageout[j-yf,i]=paddout[j,i] - imagemsaida[0:yf,0:xf]=imageout - - - #cria a imagem final, cortando a imagem filtrada (tirando as bordas) - imagemsaida_final=numpy.zeros((y,x)) - imagemsaida_final=imagemsaida[18:yf-18,18:xf-18] - - #header - header=hdu[0].header + flag = -1 + paddout[j, i] = flag * invfft[j, i] + if (i <= xf - 1) and (j > yf - 1): + imageout[j - yf, i] = paddout[j, i] + imagemsaida[0:yf, 0:xf] = imageout + + # cria a imagem final, cortando a imagem filtrada (tirando as bordas) + imagemsaida_final = numpy.zeros((y, x)) + imagemsaida_final = imagemsaida[18 : yf - 18, 18 : xf - 18] + + # header + header = hdu[0].header new_hdu = fits.PrimaryHDU() - new_hdu.header=header + new_hdu.header = header image_hdu = fits.ImageHDU(imagemsaida_final) - new_hdu.header.set('SIMPLE','T') - new_hdu.header.set('CDELT1',0.55) - new_hdu.header.set('CDELT2',0.55) - new_hdu.verify('fix') + new_hdu.header.set("SIMPLE", "T") + new_hdu.header.set("CDELT1", 0.55) + new_hdu.header.set("CDELT2", 0.55) + new_hdu.verify("fix") hdul = fits.HDUList([new_hdu, image_hdu]) - - return hdul \ No newline at end of file + + return hdul diff --git a/splusdata/features/find_pointings.py b/splusdata/features/find_pointings.py index e1c1173..d317009 100644 --- a/splusdata/features/find_pointings.py +++ b/splusdata/features/find_pointings.py @@ -4,31 +4,36 @@ import astropy.units as u from splusdata.vars import DR_POINTINGS + # uses your DR_POINTINGS exactly as given + @lru_cache(maxsize=None) def _load_dr(dr: str): """Load and normalize a DR table once; results are cached.""" info = DR_POINTINGS[dr] df = pd.read_csv(info["link"]) # keep only needed cols with consistent names - df = df.rename(columns={ - info["ra_col"]: "ra", - info["dec_col"]: "dec", - info["field_col"]: "field", - })[["ra", "dec", "field"]].copy() + df = df.rename( + columns={ + info["ra_col"]: "ra", + info["dec_col"]: "dec", + info["field_col"]: "field", + } + )[["ra", "dec", "field"]].copy() # make sure ra/dec are numeric; drop bad rows df["ra"] = pd.to_numeric(df["ra"], errors="coerce") df["dec"] = pd.to_numeric(df["dec"], errors="coerce") df = df.dropna(subset=["ra", "dec"]) return df.reset_index(drop=True) -def find_pointing(ra, dec, radius=1*u.degree): + +def find_pointing(ra, dec, radius=1 * u.degree): """ Return the first DR (by dr_order) whose field center lies within `radius`. Output: {'dr': 'dr4', 'field': '...', 'distance': } """ - target = SkyCoord(ra=ra*u.deg, dec=dec*u.deg, frame="icrs") + target = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame="icrs") radius = radius if isinstance(radius, u.Quantity) else radius * u.arcmin for dr in DR_POINTINGS.keys(): @@ -40,12 +45,14 @@ def find_pointing(ra, dec, radius=1*u.degree): return {"dr": dr, "field": df["field"].iloc[idx], "distance": sep[idx]} return None + # optional helpers: -def warm_cache(dr_order=("dr4","dr5","dr6")): +def warm_cache(dr_order=("dr4", "dr5", "dr6")): """Preload all DR tables into cache (optional).""" for dr in dr_order: _ = _load_dr(dr) + def clear_cache(): """Clear cached tables if you need to refresh.""" - _load_dr.cache_clear() \ No newline at end of file + _load_dr.cache_clear() diff --git a/splusdata/features/hipscat.py b/splusdata/features/hipscat.py index 1851b60..c2ed6c4 100644 --- a/splusdata/features/hipscat.py +++ b/splusdata/features/hipscat.py @@ -2,6 +2,7 @@ import requests import re + def _match_patterns(patterns, strings): matches = {} for string in strings: @@ -13,56 +14,59 @@ def _match_patterns(patterns, strings): matches[string] = matched_patterns return matches -def _get_hips_n_margin_links(pattern, dirs_schema, starting_path = "/"): + +def _get_hips_n_margin_links(pattern, dirs_schema, starting_path="/"): patterns = pattern.split("/") if patterns[-1] == "": patterns = patterns[:-1] res = [] - + for key in dirs_schema.keys(): if key == "hipscats": match_res = _match_patterns(patterns, dirs_schema[key]) match_res = list(match_res.keys()) - + if len(match_res) == 0: continue - + for match in match_res: - local_res = [str(starting_path)+ str(match), None] - + local_res = [str(starting_path) + str(match), None] + if "margins" in dirs_schema: for margin in dirs_schema["margins"]: if margin.startswith(match): margin = str(starting_path) + str(match) local_res[1] = margin break - + res.append(local_res) - + elif key != "margins": - res += _get_hips_n_margin_links(pattern, dirs_schema[key], starting_path = str(starting_path) + str(key)) + res += _get_hips_n_margin_links( + pattern, dirs_schema[key], starting_path=str(starting_path) + str(key) + ) return res -def get_hipscats(pattern=None, headers=None, SERVER_IP = f"https://splus.cloud"): + +def get_hipscats(pattern=None, headers=None, SERVER_IP=f"https://splus.cloud"): link = "/api/get_hipscat_available" - + if headers is None: res = requests.get(SERVER_IP + link) HIPS_IP = f"{SERVER_IP}/HIPS/catalogs" else: res = requests.get(SERVER_IP + link, headers=headers) HIPS_IP = f"{SERVER_IP}/internalhips/catalogs" - - + if res.status_code != 200: raise ValueError(f"Error: {res.status_code}") - + data = res.json() if pattern is not None: res = _get_hips_n_margin_links(pattern, data) data = res - + for i in range(len(data)): data[i] = [HIPS_IP + data[i][0], HIPS_IP + data[i][1]] - - return data \ No newline at end of file + + return data diff --git a/splusdata/features/io.py b/splusdata/features/io.py index 362e304..7acd915 100644 --- a/splusdata/features/io.py +++ b/splusdata/features/io.py @@ -3,8 +3,9 @@ from datetime import datetime from astropy.coordinates import SkyCoord + def print_level(msg, level=0, verbose=0): - ''' + """ Print a message with a specified verbosity level. Parameters @@ -18,18 +19,20 @@ def print_level(msg, level=0, verbose=0): verbose : int, optional The overall verbosity level. Messages with verbosity levels less than or equal to this value will be printed. Defaults to 0. - ''' + """ import __main__ as main + try: __script_name__ = basename(main.__file__) except AttributeError: - __script_name__ = '' + __script_name__ = "" if verbose >= level: - print(f'[{datetime.now().isoformat()}] - {__script_name__}: {msg}') + print(f"[{datetime.now().isoformat()}] - {__script_name__}: {msg}") + def check_units(ra, dec): - ''' + """ Check and add units to the input coordinates if units are missing. Parameters @@ -44,25 +47,26 @@ def check_units(ra, dec): ------- tuple Tuple containing the checked and possibly modified right ascension and declination. - ''' + """ + # Pattern to match: any letter (a-z, A-Z) or "°" - # Default unit is degrees, so if no unit, then is assumed as degs. + # Default unit is degrees, so if no unit, then is assumed as degs. def check_deg(input_string): - pattern = r'[a-zA-Z°]' + pattern = r"[a-zA-Z°]" # re.search returns a match object if the pattern is found, None otherwise if re.search(pattern, input_string): return input_string else: - return input_string + '°' + return input_string + "°" def check(input_string): - pattern = r'[hdms]' + pattern = r"[hdms]" x = re.search(pattern, input_string) if x is None: return check_deg(input_string) else: return input_string - + if not isinstance(ra, str): ra = str(ra) @@ -71,8 +75,9 @@ def check(input_string): return check(ra), check(dec) -def convert_coord_to_degrees(ra, dec, frame='icrs'): - ''' + +def convert_coord_to_degrees(ra, dec, frame="icrs"): + """ Convert the input celestial coordinates to degrees. Parameters @@ -90,11 +95,11 @@ def convert_coord_to_degrees(ra, dec, frame='icrs'): ------- tuple Tuple containing the right ascension and declination converted to degrees. - ''' + """ ra, dec = check_units(ra, dec) # Create a SkyCoord object c = SkyCoord(ra, dec, frame=frame) # Return RA and Dec in degrees - return c.ra.deg, c.dec.deg \ No newline at end of file + return c.ra.deg, c.dec.deg diff --git a/splusdata/features/zeropoints/zp_image.py b/splusdata/features/zeropoints/zp_image.py index caa5853..7ca2b82 100644 --- a/splusdata/features/zeropoints/zp_image.py +++ b/splusdata/features/zeropoints/zp_image.py @@ -6,9 +6,13 @@ import pandas as pd from astropy.table import Table -from splusdata.features.zeropoints.zp_map import _reconstruct_centers_from_model as _reconstruct_centers +from splusdata.features.zeropoints.zp_map import ( + _reconstruct_centers_from_model as _reconstruct_centers, +) + # ---------- helpers: robust interpolator from model (dict) ---------- + def _build_zp_interpolator_from_model(model): """ Returns callable f(ra, dec) -> zp (mag). @@ -18,7 +22,9 @@ def _build_zp_interpolator_from_model(model): global_median = float(model.get("global_median", 0.0)) grid = np.asarray(model.get("grid", []), dtype=float) if grid.ndim != 2: - return lambda ra, dec: np.full_like(np.asarray(ra, dtype=float), global_median, dtype=float) + return lambda ra, dec: np.full_like( + np.asarray(ra, dtype=float), global_median, dtype=float + ) ra_centers = np.asarray(model.get("ra_centers", []), dtype=float) dec_centers = np.asarray(model.get("dec_centers", []), dtype=float) @@ -32,15 +38,17 @@ def _build_zp_interpolator_from_model(model): # Try (ra, dec) with grid as-is try: interp_rd = RegularGridInterpolator( - (ra_centers, dec_centers), grid, - bounds_error=False, fill_value=np.nan + (ra_centers, dec_centers), grid, bounds_error=False, fill_value=np.nan ) + def f_rd(ra, dec): - pts = np.column_stack([np.asarray(ra, float).ravel(), - np.asarray(dec, float).ravel()]) + pts = np.column_stack( + [np.asarray(ra, float).ravel(), np.asarray(dec, float).ravel()] + ) vals = interp_rd(pts).reshape(np.shape(ra)) vals = np.where(np.isnan(vals), 0.0, vals) # local deviation fallback = 0 return global_median + vals + _ = f_rd(np.mean(ra_centers), np.mean(dec_centers)) # sanity check return f_rd except Exception: @@ -48,20 +56,26 @@ def f_rd(ra, dec): # Fallback: (dec, ra) with transposed grid interp_dr = RegularGridInterpolator( - (dec_centers, ra_centers), grid.T, - bounds_error=False, fill_value=np.nan + (dec_centers, ra_centers), grid.T, bounds_error=False, fill_value=np.nan ) + def f_dr(ra, dec): - pts = np.column_stack([np.asarray(dec, float).ravel(), - np.asarray(ra, float).ravel()]) + pts = np.column_stack( + [np.asarray(dec, float).ravel(), np.asarray(ra, float).ravel()] + ) vals = interp_dr(pts).reshape(np.shape(ra)) vals = np.where(np.isnan(vals), 0.0, vals) return global_median + vals + return f_dr + # ---------- main: in-memory calibration from HDU & model (dict) ---------- -def calibrate_hdu_with_zpmodel(hdu, model_dict, *, in_place=False, return_factor=False, safe_global_fallback=True): + +def calibrate_hdu_with_zpmodel( + hdu, model_dict, *, in_place=False, return_factor=False, safe_global_fallback=True +): """ Calibrate an already-open FITS HDU (PrimaryHDU or ImageHDU) in memory using a ZP model dict. @@ -121,18 +135,20 @@ def calibrate_hdu_with_zpmodel(hdu, model_dict, *, in_place=False, return_factor if in_place: hdu.data = calibrated # annotate header - hdu.header['ZPCALIB'] = True - hdu.header['ZPCDATE'] = datetime.datetime.now().isoformat() - hdu.header['ZPCMED'] = global_median + hdu.header["ZPCALIB"] = True + hdu.header["ZPCDATE"] = datetime.datetime.now().isoformat() + hdu.header["ZPCMED"] = global_median return (hdu, factor) if return_factor else hdu else: from astropy.io.fits import ImageHDU + new_hdu = ImageHDU(data=calibrated, header=header.copy()) - new_hdu.header['ZPCALIB'] = True - new_hdu.header['ZPCDATE'] = datetime.datetime.now().isoformat() - new_hdu.header['ZPCMED'] = global_median + new_hdu.header["ZPCALIB"] = True + new_hdu.header["ZPCDATE"] = datetime.datetime.now().isoformat() + new_hdu.header["ZPCMED"] = global_median return (new_hdu, factor) if return_factor else new_hdu - + + def compute_zp_for_coords_array( ra: np.ndarray, dec: np.ndarray, @@ -168,7 +184,9 @@ def compute_zp_for_coords_array( dec = np.asarray(dec, dtype=float) if ra.shape != dec.shape: - raise ValueError(f"ra and dec must have the same shape. Got {ra.shape} and {dec.shape}.") + raise ValueError( + f"ra and dec must have the same shape. Got {ra.shape} and {dec.shape}." + ) global_median = float(model_dict.get("global_median", 0.0)) zp_fn = _build_zp_interpolator_from_model(model_dict) @@ -182,10 +200,12 @@ def compute_zp_for_coords_array( zp[bad] = global_median else: if not np.all(np.isfinite(zp)): - raise ValueError("ZP interpolation returned NaN/inf and safe fallback is disabled.") + raise ValueError( + "ZP interpolation returned NaN/inf and safe fallback is disabled." + ) except Exception: if not safe_global_fallback: raise zp = np.full_like(ra, global_median, dtype=out_dtype) - return zp \ No newline at end of file + return zp diff --git a/splusdata/features/zeropoints/zp_map.py b/splusdata/features/zeropoints/zp_map.py index a6a25ad..95b636c 100644 --- a/splusdata/features/zeropoints/zp_map.py +++ b/splusdata/features/zeropoints/zp_map.py @@ -5,6 +5,7 @@ import warnings from scipy.interpolate import RegularGridInterpolator + def _reconstruct_centers_from_model(model, axis="ra", grid_len=None): """ Rebuild centers to match the (possibly padded) grid length using the model's @@ -64,6 +65,7 @@ def _reconstruct_centers_from_model(model, axis="ra", grid_len=None): return centers + def zp_at_coord(model, ra, dec, margin=0.1): """ Get zero-point correction(s) for coordinate(s). @@ -103,15 +105,19 @@ def zp_at_coord(model, ra, dec, margin=0.1): # Rebuild centers if needed (mismatch / empty) need_rebuild = ( - ra_centers.size == 0 or - dec_centers.size == 0 or - grid.ndim != 2 or - grid.shape[0] != ra_centers.size or - grid.shape[1] != dec_centers.size + ra_centers.size == 0 + or dec_centers.size == 0 + or grid.ndim != 2 + or grid.shape[0] != ra_centers.size + or grid.shape[1] != dec_centers.size ) if need_rebuild: - ra_centers = _reconstruct_centers_from_model(model, "ra", grid_len=grid.shape[0]) - dec_centers = _reconstruct_centers_from_model(model, "dec", grid_len=grid.shape[1]) + ra_centers = _reconstruct_centers_from_model( + model, "ra", grid_len=grid.shape[0] + ) + dec_centers = _reconstruct_centers_from_model( + model, "dec", grid_len=grid.shape[1] + ) # Normalize inputs + broadcast ra_is_scalar = np.isscalar(ra) @@ -126,8 +132,10 @@ def zp_at_coord(model, ra, dec, margin=0.1): dec_min, dec_max = float(np.min(dec_centers)), float(np.max(dec_centers)) in_bounds = ( - (ra_b >= (ra_min - margin)) & (ra_b <= (ra_max + margin)) & - (dec_b >= (dec_min - margin)) & (dec_b <= (dec_max + margin)) + (ra_b >= (ra_min - margin)) + & (ra_b <= (ra_max + margin)) + & (dec_b >= (dec_min - margin)) + & (dec_b <= (dec_max + margin)) ) if not np.all(in_bounds): @@ -188,7 +196,11 @@ def zp_at_coord(model, ra, dec, margin=0.1): global_median = float(model.get("global_median", 0.0)) # If no grid, always fallback - if not ("grid" in model and ("ra_centers" in model or True) and ("dec_centers" in model or True)): + if not ( + "grid" in model + and ("ra_centers" in model or True) + and ("dec_centers" in model or True) + ): # keep original behavior: return global median only if np.isscalar(ra) and np.isscalar(dec) and return_scalar_if_scalar: return global_median @@ -204,15 +216,19 @@ def zp_at_coord(model, ra, dec, margin=0.1): # Rebuild centers if needed need_rebuild = ( - ra_centers.size == 0 or - dec_centers.size == 0 or - grid.ndim != 2 or - grid.shape[0] != ra_centers.size or - grid.shape[1] != dec_centers.size + ra_centers.size == 0 + or dec_centers.size == 0 + or grid.ndim != 2 + or grid.shape[0] != ra_centers.size + or grid.shape[1] != dec_centers.size ) if need_rebuild: - ra_centers = _reconstruct_centers_from_model(model, "ra", grid_len=grid.shape[0]) - dec_centers = _reconstruct_centers_from_model(model, "dec", grid_len=grid.shape[1]) + ra_centers = _reconstruct_centers_from_model( + model, "ra", grid_len=grid.shape[0] + ) + dec_centers = _reconstruct_centers_from_model( + model, "dec", grid_len=grid.shape[1] + ) # Interpolator (build once) interpolator = RegularGridInterpolator( @@ -235,8 +251,10 @@ def zp_at_coord(model, ra, dec, margin=0.1): dec_min, dec_max = float(np.min(dec_centers)), float(np.max(dec_centers)) in_bounds = ( - (ra_b >= (ra_min - margin)) & (ra_b <= (ra_max + margin)) & - (dec_b >= (dec_min - margin)) & (dec_b <= (dec_max + margin)) + (ra_b >= (ra_min - margin)) + & (ra_b <= (ra_max + margin)) + & (dec_b >= (dec_min - margin)) + & (dec_b <= (dec_max + margin)) ) if not np.all(in_bounds): @@ -268,4 +286,4 @@ def zp_at_coord(model, ra, dec, margin=0.1): if return_scalar_if_scalar and ra_is_scalar and dec_is_scalar: return float(zp) - return zp \ No newline at end of file + return zp diff --git a/splusdata/features/zeropointsdr4.py b/splusdata/features/zeropointsdr4.py index cef858c..f9d6add 100644 --- a/splusdata/features/zeropointsdr4.py +++ b/splusdata/features/zeropointsdr4.py @@ -1,16 +1,19 @@ import pandas as pd try: - source_cat = "https://splus.cloud/files/documentation/iDR4/tabelas/iDR4_zero-points.csv" + source_cat = ( + "https://splus.cloud/files/documentation/iDR4/tabelas/iDR4_zero-points.csv" + ) import requests - import io - + import io + r = requests.get(source_cat, timeout=5) zps = pd.read_csv(io.StringIO(r.text)) except Exception as e: - #print(f"Error loading zero points data: {e}") + # print(f"Error loading zero points data: {e}") pass + def get_zeropoint(conn, ra, dec, band): """ Retrieve the zero point for a given sky position and band. @@ -24,7 +27,7 @@ def get_zeropoint(conn, ra, dec, band): dec : float Declination of the target position (in degrees). band : str - Photometric band for which to retrieve the zero point (e.g., 'g', 'r', 'i'). + Photometric band for which to retrieve the zero point (e.g., 'g', 'r', 'i'). If set to 'all', returns zero points for all bands. Returns @@ -37,7 +40,7 @@ def get_zeropoint(conn, ra, dec, band): ------ KeyError If the provided field or band is not found in the zero point table. - + Examples -------- >>> conn = some_connection_object() @@ -51,16 +54,15 @@ def get_zeropoint(conn, ra, dec, band): Field ZP_u ZP_g ZP_r ... 1234_56789 25.13 24.98 25.02 ... """ - - field = list(conn.checkcoords(ra, dec).values())[0] - + + field = list(conn.checkcoords(ra, dec).values())[0] + zps["Field"] = zps["Field"].str.replace("-", "_") field = field.replace("-", "_") - + if band.lower() == "all": zp = zps[zps["Field"] == field] else: zp = zps[zps["Field"] == field]["ZP_" + band].values[0] - + return zp - \ No newline at end of file diff --git a/splusdata/models/star_gal_quasar.py b/splusdata/models/star_gal_quasar.py index c86f701..8fbf170 100644 --- a/splusdata/models/star_gal_quasar.py +++ b/splusdata/models/star_gal_quasar.py @@ -9,25 +9,30 @@ __author__ = "Lilianne Nakazono" __version__ = "0.1.0" + class ClassifyObj: - _morph = ['FWHM_n', 'A', 'B', 'KRON_RADIUS'] - _feat = ['u_iso', - 'J0378_iso', - 'J0395_iso', - 'J0410_iso', - 'J0430_iso', - 'g_iso', - 'J0515_iso', - 'r_iso', - 'J0660_iso', - 'i_iso', - 'J0861_iso', - 'z_iso'] - _feat_wise = ['w1mpro', 'w2mpro'] - _error_wise = ['w1snr', 'w2snr', 'w1sigmpro', 'w2sigmpro'] + _morph = ["FWHM_n", "A", "B", "KRON_RADIUS"] + _feat = [ + "u_iso", + "J0378_iso", + "J0395_iso", + "J0410_iso", + "J0430_iso", + "g_iso", + "J0515_iso", + "r_iso", + "J0660_iso", + "i_iso", + "J0861_iso", + "z_iso", + ] + _feat_wise = ["w1mpro", "w2mpro"] + _error_wise = ["w1snr", "w2snr", "w1sigmpro", "w2sigmpro"] _pos = ["RA", "DEC"] - def __init__(self, data, model = 'RF16', return_prob = True, match_irsa = False, verbose=True): + def __init__( + self, data, model="RF16", return_prob=True, match_irsa=False, verbose=True + ): self.data = data self.model = model self.return_prob = return_prob @@ -37,7 +42,7 @@ def __init__(self, data, model = 'RF16', return_prob = True, match_irsa = False, @staticmethod def get_wise(df, list_names): - ''' + """ Deletes objects without WISE counterpart (> 2 arcsecond) from data Keywords arguments: @@ -45,7 +50,7 @@ def get_wise(df, list_names): list_names -- names of the columns corresponding to W1 and W2 magnitudes/error/snr from WISE returns a copy of a subset from data with objects that have WISE counterpart - ''' + """ copy = df.copy(deep=True) copy = copy.query("d2d <= 2") @@ -54,19 +59,28 @@ def get_wise(df, list_names): return copy def irsa_query(self): - ''' + """ Query ALLWISE catalogue from IRSA database Keywords arguments: data -- dataframe containing information of RA and DEC for each object return dataframe containing WISE sources - ''' + """ ramin = np.min(self.data["RA"]) ramax = np.max(self.data["RA"]) decmin = np.min(self.data["DEC"]) decmax = np.max(self.data["DEC"]) - col = ['ra', 'dec', 'w1mpro', 'w2mpro', 'w1snr', 'w2snr', 'w1sigmpro', 'w2sigmpro'] + col = [ + "ra", + "dec", + "w1mpro", + "w2mpro", + "w1snr", + "w2snr", + "w1sigmpro", + "w2sigmpro", + ] if ramin < 2 and ramax > 358: ramax = np.max(self.data.query("RA<2")) ramin = np.min(self.data.query("RA>358")) @@ -91,17 +105,25 @@ def irsa_query(self): # Renaming any original column in data that has the same name from query for feat in col: if feat in self.data.columns: - self.data.rename(columns={feat: feat+"_"+str(0)}, inplace=True) + self.data.rename(columns={feat: feat + "_" + str(0)}, inplace=True) df_wise = df_wise.rename( - columns={"col_0": col[0], "col_1": col[1], "col_2": col[2], - "col_3": col[3], "col_4": col[4], "col_5": col[5], - "col_6": col[6], "col_7": col[7]}) + columns={ + "col_0": col[0], + "col_1": col[1], + "col_2": col[2], + "col_3": col[3], + "col_4": col[4], + "col_5": col[5], + "col_6": col[6], + "col_7": col[7], + } + ) return df_wise def crossmatch(self, df_wise): - ''' + """ Perform cross match between two catalogues based on RA and DEC data_ra -- RA from S-PLUS catalogue @@ -110,40 +132,51 @@ def crossmatch(self, df_wise): df_wise_dec -- DEC from ALLWISE catalogue returns dataframe containing all original columns from data cross-matched with the ALLWISE query catalogue - ''' + """ - coo_splus = SkyCoord(self.data["RA"]* u.deg, self.data["DEC"] * u.deg) - coo_wise = SkyCoord(df_wise["ra"]* u.deg, df_wise["dec"] * u.deg) + coo_splus = SkyCoord(self.data["RA"] * u.deg, self.data["DEC"] * u.deg) + coo_wise = SkyCoord(df_wise["ra"] * u.deg, df_wise["dec"] * u.deg) idx, d2d, d3d = coo_splus.match_to_catalog_sky(coo_wise) - self.data['d2d'] = pd.Series(d2d.arcsec) + self.data["d2d"] = pd.Series(d2d.arcsec) self.data = pd.concat([self.data, df_wise.iloc[idx].reset_index()], axis=1) return self.data def check_match_irsa(self): - ''' + """ Check statement of match_irsa and performs cross-match if True data -- dataframe containing S-PLUS information match_irsa -- If true, query ALLWISE catalogue and performs cross-match. If false, checks if all necessary columns are in data - ''' - + """ if self.match_irsa == False: try: for element in self._feat_wise + self._error_wise: i = self.data[element] except ValueError: - print("Please ensure that ", element, - "column is in data. Please provide a ALLWISE cross-matched data or set match_irsa == True. Use model == 'opt' if the provided data do not have any source with WISE counteerpart.") + print( + "Please ensure that ", + element, + "column is in data. Please provide a ALLWISE cross-matched data or set match_irsa == True. Use model == 'opt' if the provided data do not have any source with WISE counteerpart.", + ) else: try: - for element in self._pos: #check if RA and DEC are in data to perfom the cross-match + for ( + element + ) in ( + self._pos + ): # check if RA and DEC are in data to perfom the cross-match i = self.data[element] except ValueError: - print("Please ensure that ", element, - "column is in data, otherwise the cross-match with ALLWISE cannot be done.") + print( + "Please ensure that ", + element, + "column is in data, otherwise the cross-match with ALLWISE cannot be done.", + ) if self.verbose: - print("Querying ALLWISE catalogue via TapPlus(url='https://irsa.ipac.caltech.edu/TAP')...") + print( + "Querying ALLWISE catalogue via TapPlus(url='https://irsa.ipac.caltech.edu/TAP')..." + ) df_wise = self.irsa_query() if self.verbose: print("Starting cross-match...") @@ -151,20 +184,19 @@ def check_match_irsa(self): return self.data def classificator(self): - - ''' + """ Create classifications for sources with or without counterpart - + Keywords arguments: data -- dataframe containing information of the 12 S-PLUS ISO magnitudes already extincted corrected prob -- if true, estimates the probabilities for each class model -- options are "both", "RF16" or "RF18". If "opt", returns classification from model trained with 12 S-PLUS bands + 4 morphological features. If "RF18", returns classification from model trained with 12 S-PLUS bands + 2 WISE bands+ 4 morphological features, only for sources with WISE counterpart. If "both", return classification from the same model as "RF18" (flagged as 0) if - the source has WISE counterpart, otherwise returns classification from model "RF16" (flagged as 1). + the source has WISE counterpart, otherwise returns classification from model "RF16" (flagged as 1). match_irsa -- determines if cross-match with ALLWISE catalogue will be performed. If model == "RF16", match_irsa == False. - returns a dataframe with classes - ''' + returns a dataframe with classes + """ if self.model == "RF16": if self.match_irsa: @@ -174,41 +206,50 @@ def classificator(self): try: i = self.data.iloc[0,] except: - raise(ValueError)("Data input is empty.") + raise (ValueError)("Data input is empty.") try: - if self.model == 'RF16': + if self.model == "RF16": print("Loading model") model = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S4M") model = model.content model = pickle.loads(model) - elif self.model == 'RF18': + elif self.model == "RF18": print("Loading model") - model_wise = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S2W4M") + model_wise = requests.get( + "https://splus.cloud/files/models/iDR3n4_RF_12S2W4M" + ) model_wise = model_wise.content model_wise = pickle.loads(model_wise) - else: + else: print("Loading both models") model = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S4M") model = model.content model = pickle.loads(model) - model_wise = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S2W4M") + model_wise = requests.get( + "https://splus.cloud/files/models/iDR3n4_RF_12S2W4M" + ) model_wise = model_wise.content model_wise = pickle.loads(model_wise) except: - raise(ValueError)("Error downloading model from splus.cloud") - + raise (ValueError)("Error downloading model from splus.cloud") try: - for element in self._morph+self._feat: + for element in self._morph + self._feat: i = self.data[element] except ValueError: - print("Please ensure that ", element, "column is in data. Use splusdata.connect.query() to retrieve it." ) + print( + "Please ensure that ", + element, + "column is in data. Use splusdata.connect.query() to retrieve it.", + ) self.results = pd.DataFrame() if self.verbose: - print("Note that this function assumes that the input data were previously extinction corrected. Starting classification... ") + print( + "Note that this function assumes that the input data were previously extinction corrected. Starting classification... " + ) if self.model == "RF16": y_pred = pd.DataFrame(model.predict(self.data[self._morph + self._feat])) @@ -221,7 +262,9 @@ def classificator(self): if self.return_prob: if self.verbose: print("Calculating probabilities...") - prob_df = pd.DataFrame(model.predict_proba(self.data[self._morph + self._feat])) + prob_df = pd.DataFrame( + model.predict_proba(self.data[self._morph + self._feat]) + ) prob_df.index = self.data.index results = pd.concat([results, prob_df], axis=1) @@ -229,8 +272,12 @@ def classificator(self): elif self.model == "RF18": self.check_match_irsa() - data_wise = self.get_wise(self.data, self._feat_wise+self._error_wise) - ypred_wise = pd.DataFrame(model_wise.predict(data_wise[self._morph + self._feat + self._feat_wise])) # classify sources with WISE + data_wise = self.get_wise(self.data, self._feat_wise + self._error_wise) + ypred_wise = pd.DataFrame( + model_wise.predict( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) # classify sources with WISE ypred_wise.index = data_wise.index results = ypred_wise results.columns = ["CLASS"] @@ -239,14 +286,18 @@ def classificator(self): if self.verbose: print("Calculating probabilities...") - prob_wise_df = pd.DataFrame(model_wise.predict_proba(data_wise[self._morph + self._feat + self._feat_wise])) + prob_wise_df = pd.DataFrame( + model_wise.predict_proba( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) prob_wise_df.index = data_wise.index results = pd.concat([results, prob_wise_df], axis=1) results.columns = ["CLASS", "PROB_QSO", "PROB_STAR", "PROB_GAL"] - elif self.model=="both": + elif self.model == "both": self.check_match_irsa() - data_wise = self.get_wise(self.data, self._feat_wise+self._error_wise) + data_wise = self.get_wise(self.data, self._feat_wise + self._error_wise) data_nowise = self.data.drop(data_wise.index) y_pred = pd.DataFrame(model.predict(data_nowise[self._morph + self._feat])) @@ -254,7 +305,11 @@ def classificator(self): y_pred.index = data_nowise.index y_pred["model_flag"] = "1" # S-PLUS - ypred_wise = pd.DataFrame(model_wise.predict(data_wise[self._morph + self._feat + self._feat_wise])) # classify sources with WISE + ypred_wise = pd.DataFrame( + model_wise.predict( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) # classify sources with WISE ypred_wise.index = data_wise.index ypred_wise["model_flag"] = "0" # S-PLUS + WISE results = pd.concat([y_pred, ypred_wise], axis=0) @@ -263,18 +318,30 @@ def classificator(self): if self.return_prob: if self.verbose: print("Calculating probabilities...") - prob_df = pd.DataFrame(model.predict_proba(data_nowise[self._morph + self._feat])) - prob_wise_df = pd.DataFrame(model_wise.predict_proba(data_wise[self._morph + self._feat + self._feat_wise])) + prob_df = pd.DataFrame( + model.predict_proba(data_nowise[self._morph + self._feat]) + ) + prob_wise_df = pd.DataFrame( + model_wise.predict_proba( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) prob_df.index = data_nowise.index prob_wise_df.index = data_wise.index prob_results = pd.concat([prob_df, prob_wise_df], axis=0) results = pd.concat([results, prob_results], axis=1) - results.columns = ["CLASS", "model_flag", "PROB_QSO", "PROB_STAR", "PROB_GAL"] + results.columns = [ + "CLASS", + "model_flag", + "PROB_QSO", + "PROB_STAR", + "PROB_GAL", + ] else: - raise(ValueError)("Parameter 'model' should be 'RF16', 'RF18' or 'both'.") + raise (ValueError)("Parameter 'model' should be 'RF16', 'RF18' or 'both'.") if self.verbose: print("Finished process.") - return results.sort_index(axis=0) \ No newline at end of file + return results.sort_index(axis=0) diff --git a/splusdata/readconf.py b/splusdata/readconf.py index ba33043..51cb601 100644 --- a/splusdata/readconf.py +++ b/splusdata/readconf.py @@ -7,20 +7,21 @@ from astropy.table import Table from yaml import load, dump + try: from yaml import CLoader as Loader, CDumper as Dumper except ImportError: from yaml import Loader, Dumper - + operations = [ - 'stamp', - 'lupton_rgb', - 'trilogy_image' - 'field_frame', - 'checkcoords', - 'query' + "stamp", + "lupton_rgb", + "trilogy_image" "field_frame", + "checkcoords", + "query", ] + def filter_info(info, function): # Get the signature of the function func_signature = inspect.signature(function) @@ -30,99 +31,113 @@ def filter_info(info, function): # Filter the dictionary filtered_dict = {k: info[k] for k in expected_keys if k in info} - + return filtered_dict + def create_stamp_name(info, output_folder): - filename = os.path.join(output_folder, f"{info['ra']}_{info['dec']}_{info['size']}_{info['band']}.fits") + filename = os.path.join( + output_folder, f"{info['ra']}_{info['dec']}_{info['size']}_{info['band']}.fits" + ) return filename + def create_image_name(info, output_folder): - filename = os.path.join(output_folder, f"{info['ra']}_{info['dec']}_{info['size']}.png") + filename = os.path.join( + output_folder, f"{info['ra']}_{info['dec']}_{info['size']}.png" + ) return filename + def create_field_name(info, output_folder): filename = os.path.join(output_folder, f"{info['field']}_{info['band']}.fits") return filename - -def handle_operation(operation : str, info, conn): - output_folder = info['output_folder'] + + +def handle_operation(operation: str, info, conn): + output_folder = info["output_folder"] if not os.path.exists(output_folder): os.makedirs(output_folder) - - if operation == 'stamp': + + if operation == "stamp": filename = create_stamp_name(info, output_folder) info = filter_info(info, conn.stamp) - conn.stamp( filename=filename, **info ) - - elif operation == 'lupton_rgb': + conn.stamp(filename=filename, **info) + + elif operation == "lupton_rgb": filename = create_image_name(info, output_folder) info = filter_info(info, conn.lupton_rgb) - conn.lupton_rgb( filename=filename, **info ) - - elif operation == 'trilogy_image': + conn.lupton_rgb(filename=filename, **info) + + elif operation == "trilogy_image": filename = create_image_name(info, output_folder) info = filter_info(info, conn.trilogy_image) - conn.trilogy_image( filename=filename, **info ) + conn.trilogy_image(filename=filename, **info) - elif operation == 'field_frame': + elif operation == "field_frame": filename = create_field_name(info, output_folder) info = filter_info(info, conn.field_frame) - conn.field_frame( filename=filename, **info ) + conn.field_frame(filename=filename, **info) - elif operation == 'checkcoords': + elif operation == "checkcoords": info = filter_info(info, conn.checkcoords) try: - data = conn.checkcoords( **info ) + data = conn.checkcoords(**info) data["ra"] = info["ra"] data["dec"] = info["dec"] except: - data = {"ra": info["ra"], "dec": info["dec"], "field": None, "distance": None, "public": None, "internal": None} + data = { + "ra": info["ra"], + "dec": info["dec"], + "field": None, + "distance": None, + "public": None, + "internal": None, + } return data - elif operation == 'query': + elif operation == "query": if "query" not in info: raise ValueError("Query not specified") - + filename = os.path.join(output_folder, "query.fits") if "upload_table" in info: if not os.path.exists(info["upload_table"]): raise ValueError("File {} does not exist".format(info["upload_table"])) - + try: file = pd.read_csv(info["upload_table"]) info = filter_info(info, conn.query) - data = conn.query(table_upload=file, **info ) + data = conn.query(table_upload=file, **info) data.write(filename, overwrite=True) except Exception as e: print(e) else: - + info = filter_info(info, conn.query) try: - data = conn.query( **info ) + data = conn.query(**info) data.write(filename, overwrite=True) except Exception as e: print(e) - -def get_coordinates_col(df : pd.DataFrame): +def get_coordinates_col(df: pd.DataFrame): ra = None dec = None - + for col in df.columns: if col.lower() == "ra": ra = col elif col.lower() == "dec": dec = col - + if ra is None or dec is None: raise ValueError("File does not have columns 'ra' and 'dec'") - + return ra, dec - + def handle_operation_type(operation_data, conn): checkcoords_list = [] @@ -132,11 +147,11 @@ def handle_operation_type(operation_data, conn): operations = operation_data["type"] else: raise ValueError("Operation type must be a string or a list of strings") - + for op in operations: if op not in operations: raise ValueError("Operation type not supported: {}".format(op)) - + if "file_path" in operation_data: try: df = pd.read_csv(operation_data["file_path"]) @@ -155,7 +170,7 @@ def handle_operation_type(operation_data, conn): checkcoords_list.append(handle_operation(op, info, conn)) else: handle_operation(op, info, conn) - + else: for op in operations: if op == "checkcoords": @@ -165,30 +180,35 @@ def handle_operation_type(operation_data, conn): if len(checkcoords_list) > 0: checkcoords_df = pd.DataFrame(checkcoords_list) - checkcoords_df.to_csv(os.path.join(operation_data["output_folder"], "checkcoords.csv"), index=False) - + checkcoords_df.to_csv( + os.path.join(operation_data["output_folder"], "checkcoords.csv"), + index=False, + ) def main(): - parser = argparse.ArgumentParser(description='splusdata - Download SPLUS catalogs, FITS and more') - parser.add_argument('config_file', metavar='config_file', type=str, help='Configuration file') - + parser = argparse.ArgumentParser( + description="splusdata - Download SPLUS catalogs, FITS and more" + ) + parser.add_argument( + "config_file", metavar="config_file", type=str, help="Configuration file" + ) + args = parser.parse_args() - + configfile = os.path.join(os.getcwd(), args.config_file) - stream = open(configfile, 'rb') + stream = open(configfile, "rb") data = load(stream, Loader=Loader) try: conn = splusdata.Core(data["user"], data["password"]) except: conn = splusdata.Core() - + handle_operation_type(data["operation"], conn) return 0 - -if '__name__' == '__main__': - main() \ No newline at end of file +if "__name__" == "__main__": + main() diff --git a/splusdata/scripts/args.py b/splusdata/scripts/args.py index 09ca64d..6ee292e 100644 --- a/splusdata/scripts/args.py +++ b/splusdata/scripts/args.py @@ -1,13 +1,14 @@ from argparse import ArgumentParser, RawDescriptionHelpFormatter + class readFileArgumentParser(ArgumentParser): - ''' + """ A class that extends :class:`argparse.ArgumentParser` to read arguments from file. Methods ------- convert_arg_line_to_args: yeld's arguments reading a file. - ''' + """ def __init__(self, *args, **kwargs): super(readFileArgumentParser, self).__init__(*args, **kwargs) @@ -16,16 +17,17 @@ def convert_arg_line_to_args(self, line): for arg in line.split(): if not arg.strip(): continue - if arg[0] == '#': + if arg[0] == "#": break yield arg + def create_parser(args_dict, program_description=None): - ''' - Create the parser for the arguments keys on `args_dict` to `s-cubes` + """ + Create the parser for the arguments keys on `args_dict` to `s-cubes` entry-points console scripts. - It uses :class:`readFileArgumentParser` to include the option to read the + It uses :class:`readFileArgumentParser` to include the option to read the arguments from a file using the prefix @ before the filename:: entrypoint @file.args @@ -33,8 +35,8 @@ def create_parser(args_dict, program_description=None): Parameters ---------- args_dict : dict - Dictionary to with the long options as keys and a list containing the short - option char at the first position and the `kwargs` for the argument + Dictionary to with the long options as keys and a list containing the short + option char at the first position and the `kwargs` for the argument configuration. program_description : str, optional @@ -48,30 +50,30 @@ def create_parser(args_dict, program_description=None): See Also -------- `argparse.ArgumentParser.add_argument()` - ''' + """ _formatter = lambda prog: RawDescriptionHelpFormatter(prog, max_help_position=30) - + parser = readFileArgumentParser( - fromfile_prefix_chars='@', - description=program_description, - formatter_class=_formatter + fromfile_prefix_chars="@", + description=program_description, + formatter_class=_formatter, ) for k, v in args_dict.items(): long_option = k short_option, kwargs = v option_string = [] positional = False - if short_option != '': - if short_option == 'pos': + if short_option != "": + if short_option == "pos": positional = True else: - option_string.append(f'-{short_option}') + option_string.append(f"-{short_option}") if positional: option_string.append(long_option) else: - option_string.append(f'--{long_option}') - if (long_option != 'verbose') and (kwargs.get('default') is not None): - _tmp = kwargs['help'] - kwargs['help'] = f'{_tmp} Default value is %(default)s' + option_string.append(f"--{long_option}") + if (long_option != "verbose") and (kwargs.get("default") is not None): + _tmp = kwargs["help"] + kwargs["help"] = f"{_tmp} Default value is %(default)s" parser.add_argument(*option_string, **kwargs) - return parser \ No newline at end of file + return parser diff --git a/splusdata/scripts/utils.py b/splusdata/scripts/utils.py index 637f337..d423d0a 100644 --- a/splusdata/scripts/utils.py +++ b/splusdata/scripts/utils.py @@ -8,26 +8,38 @@ from splusdata.scripts.args import create_parser from splusdata import Core - FIELD_FRAME_ARGS = { # optional arguments - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS cloud password.')], - 'weight': ['w', dict(action='store_true', default=False, help='Download the weight map associated to the stamp.')], - 'data_release': ['D', dict(default='DR4', help='S-PLUS data release version.')], - + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS cloud password.")], + "weight": [ + "w", + dict( + action="store_true", + default=False, + help="Download the weight map associated to the stamp.", + ), + ], + "data_release": ["D", dict(default="DR4", help="S-PLUS data release version.")], # positional arguments - 'field': ['pos', dict(metavar='FIELD', help='Field identifier, e.g., "SPLUS-n01s10"')], - 'band': ['pos', dict(metavar='BAND', help='S-PLUS band, e.g., "R", "I", "F660", "U".')], + "field": [ + "pos", + dict(metavar="FIELD", help='Field identifier, e.g., "SPLUS-n01s10"'), + ], + "band": [ + "pos", + dict(metavar="BAND", help='S-PLUS band, e.g., "R", "I", "F660", "U".'), + ], } -FIELD_FRAME_PROG_DESC = f''' +FIELD_FRAME_PROG_DESC = f""" Download and open a full field FITS image. -''' +""" + def field_frame(): - ''' + """ Script function for downloading a S-PLUS field FITS image. Raises @@ -38,40 +50,79 @@ def field_frame(): Returns ------- None - ''' - parser = create_parser(args_dict=FIELD_FRAME_ARGS, program_description=FIELD_FRAME_PROG_DESC) + """ + parser = create_parser( + args_dict=FIELD_FRAME_ARGS, program_description=FIELD_FRAME_PROG_DESC + ) args = parser.parse_args(args=sys.argv[1:]) conn = Core(args.username, args.password) - print_level(f'Downloading field {args.field} - band {args.band}', 1, args.verbose) - outfile = f'{args.field}_{args.band}.fits.fz' - conn.field_frame(args.field, args.band, weight=args.weight, data_release=args.data_release, outfile=outfile) + print_level(f"Downloading field {args.field} - band {args.band}", 1, args.verbose) + outfile = f"{args.field}_{args.band}.fits.fz" + conn.field_frame( + args.field, + args.band, + weight=args.weight, + data_release=args.data_release, + outfile=outfile, + ) + STAMP_ARGS = { # optional arguments - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS cloud password.')], - 'weight': ['w', dict(action='store_true', default=False, help='Download the weight map associated to the stamp.')], - 'field_name': ['F', dict(metavar='FIELD', help='Field identifier, e.g., "SPLUS-n01s10"')], - 'size_unit': ['u', dict(metavar='SIZE_UNIT', choices=['arcsec', 'arcmin', 'pixels'], default='pix', help='Unit of the size parameter (arcsec, arcmin, pixels).')], - 'data_release': ['D', dict(default='DR4', help='S-PLUS data release version.')], - 'outfile': ['o', dict(default='stamp.fits.fz', help='Output filename for the stamp FITS image.')], - + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS cloud password.")], + "weight": [ + "w", + dict( + action="store_true", + default=False, + help="Download the weight map associated to the stamp.", + ), + ], + "field_name": [ + "F", + dict(metavar="FIELD", help='Field identifier, e.g., "SPLUS-n01s10"'), + ], + "size_unit": [ + "u", + dict( + metavar="SIZE_UNIT", + choices=["arcsec", "arcmin", "pixels"], + default="pix", + help="Unit of the size parameter (arcsec, arcmin, pixels).", + ), + ], + "data_release": ["D", dict(default="DR4", help="S-PLUS data release version.")], + "outfile": [ + "o", + dict(default="stamp.fits.fz", help="Output filename for the stamp FITS image."), + ], # positional arguments - 'ra': ['pos', dict(metavar='RA', help="Object's right ascension")], - 'dec': ['pos', dict(metavar='DEC', help="Object's declination")], - 'size': ['pos', dict(metavar='SIZE', help='Size of the stamp in size_unit (defaults to pixels) (see --size_unit).')], - 'band': ['pos', dict(metavar='BAND', help='S-PLUS band, e.g., "R", "I", "F660", "U".')], + "ra": ["pos", dict(metavar="RA", help="Object's right ascension")], + "dec": ["pos", dict(metavar="DEC", help="Object's declination")], + "size": [ + "pos", + dict( + metavar="SIZE", + help="Size of the stamp in size_unit (defaults to pixels) (see --size_unit).", + ), + ], + "band": [ + "pos", + dict(metavar="BAND", help='S-PLUS band, e.g., "R", "I", "F660", "U".'), + ], } -STAMP_PROG_DESC = f''' +STAMP_PROG_DESC = f""" Download a FITS stamp (cutout) by coordinates or by object name. -''' +""" + def stamp(): - ''' + """ Script function for downloading a S-PLUS stamp FITS image. Raises @@ -82,41 +133,86 @@ def stamp(): Returns ------- None - ''' + """ parser = create_parser(args_dict=STAMP_ARGS, program_description=STAMP_PROG_DESC) args = parser.parse_args(args=sys.argv[1:]) conn = Core(args.username, args.password) - print_level(f'Downloading stamp at RA: {args.ra}, DEC: {args.dec} - band {args.band}', 1, args.verbose) - conn.stamp(args.ra, args.dec, args.size, args.band, - weight=args.weight, field_name=args.field_name, - size_unit=args.size_unit, data_release=args.data_release, outfile=args.outfile, + print_level( + f"Downloading stamp at RA: {args.ra}, DEC: {args.dec} - band {args.band}", + 1, + args.verbose, + ) + conn.stamp( + args.ra, + args.dec, + args.size, + args.band, + weight=args.weight, + field_name=args.field_name, + size_unit=args.size_unit, + data_release=args.data_release, + outfile=args.outfile, ) -LUPTON_RGB_ARGS = { - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS cloud password.')], - 'rgb': ['b', dict(metavar='BAND1,BAND2,BAND3', nargs=3, default='I,R,G', help='Space-separated bands for RGB channels (default: I R G).')], - 'field_name': ['F', dict(metavar='FIELD', help='Field identifier, e.g., "SPLUS-n01s10"')], - 'size_unit': ['u', dict(metavar='SIZE_UNIT', choices=['arcsec', 'arcmin', 'pixels'], default='pix', help='Unit of the size parameter (arcsec, arcmin, pixels).')], - 'data_release': ['D', dict(default='DR4', help='S-PLUS data release version.')], - 'outfile': ['o', dict(default='lupton_RGB.png', help='Output filename for the Lupton RGB image.')], - 'Q': ['Q', dict(type=float, default=8, help='Lupton Q parameter (default: 8).')], - 'stretch': ['s', dict(type=float, default=3, help='Lupton stretch parameter (default: 3).')], +LUPTON_RGB_ARGS = { + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS cloud password.")], + "rgb": [ + "b", + dict( + metavar="BAND1,BAND2,BAND3", + nargs=3, + default="I,R,G", + help="Space-separated bands for RGB channels (default: I R G).", + ), + ], + "field_name": [ + "F", + dict(metavar="FIELD", help='Field identifier, e.g., "SPLUS-n01s10"'), + ], + "size_unit": [ + "u", + dict( + metavar="SIZE_UNIT", + choices=["arcsec", "arcmin", "pixels"], + default="pix", + help="Unit of the size parameter (arcsec, arcmin, pixels).", + ), + ], + "data_release": ["D", dict(default="DR4", help="S-PLUS data release version.")], + "outfile": [ + "o", + dict( + default="lupton_RGB.png", help="Output filename for the Lupton RGB image." + ), + ], + "Q": ["Q", dict(type=float, default=8, help="Lupton Q parameter (default: 8).")], + "stretch": [ + "s", + dict(type=float, default=3, help="Lupton stretch parameter (default: 3)."), + ], # positional arguments - 'ra': ['pos', dict(metavar='RA', help="Object's right ascension")], - 'dec': ['pos', dict(metavar='DEC', help="Object's declination")], - 'size': ['pos', dict(metavar='SIZE', help='Size of the stamp in size_unit (defaults to pixels) (see --size_unit).')], + "ra": ["pos", dict(metavar="RA", help="Object's right ascension")], + "dec": ["pos", dict(metavar="DEC", help="Object's declination")], + "size": [ + "pos", + dict( + metavar="SIZE", + help="Size of the stamp in size_unit (defaults to pixels) (see --size_unit).", + ), + ], } -LUPTON_RGB_PROG_DESC = f''' +LUPTON_RGB_PROG_DESC = f""" Create a Lupton RGB composite. -''' +""" + def lupton_rgb_argparse(args): - ''' + """ A particular parser of the command-line arguments for `lupton_rgb` script. Parameters @@ -128,13 +224,14 @@ def lupton_rgb_argparse(args): ------- :class:`argparse.Namespace` Command-line arguments parsed. - ''' + """ # split RGB bands args.R, args.G, args.B = args.rgb return args + def lupton_rgb(): - ''' + """ Script function for downloading a S-PLUS Lupton RGB image. Raises @@ -145,46 +242,122 @@ def lupton_rgb(): Returns ------- None - ''' - parser = create_parser(args_dict=LUPTON_RGB_ARGS, program_description=LUPTON_RGB_PROG_DESC) + """ + parser = create_parser( + args_dict=LUPTON_RGB_ARGS, program_description=LUPTON_RGB_PROG_DESC + ) args = lupton_rgb_argparse(args=parser.parse_args(args=sys.argv[1:])) conn = Core(args.username, args.password) - print_level(f'Downloading Lupton RGB image at RA: {args.ra}, DEC: {args.dec}', 1, args.verbose) - conn.lupton_rgb(args.ra, args.dec, args.size, R=args.R, G=args.G, B=args.B, - field_name=args.field_name, size_unit=args.size_unit, - data_release=args.data_release, outfile=args.outfile, - Q=args.Q, stretch=args.stretch + print_level( + f"Downloading Lupton RGB image at RA: {args.ra}, DEC: {args.dec}", + 1, + args.verbose, + ) + conn.lupton_rgb( + args.ra, + args.dec, + args.size, + R=args.R, + G=args.G, + B=args.B, + field_name=args.field_name, + size_unit=args.size_unit, + data_release=args.data_release, + outfile=args.outfile, + Q=args.Q, + stretch=args.stretch, ) -TRILOGY_IMAGE_ARGS = { - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS cloud password.')], - 'red': ['R', dict(metavar='RED_BAND', default='I', nargs='+', help='Band(s) for the red channel (default: I R G).')], - 'green': ['G', dict(metavar='GREEN_BAND', default='R', nargs='+', help='Band(s) for the green channel (default: R G B).')], - 'blue': ['B', dict(metavar='BLUE_BAND', default='G', nargs='+', help='Band(s) for the blue channel (default: G B I).')], - 'field_name': ['F', dict(metavar='FIELD', help='Field identifier, e.g., "SPLUS-n01s10"')], - 'size_unit': ['u', dict(metavar='SIZE_UNIT', choices=['arcsec', 'arcmin', 'pixels'], default='pix', help='Unit of the size parameter (arcsec, arcmin, pixels).')], - 'data_release': ['D', dict(default='DR4', help='S-PLUS data release version.')], - 'outfile': ['o', dict(default='trilogy_image.png', help='Output filename for the Lupton RGB image.')], - 'noiselum': ['n', dict(type=float, default=0.15, help='Controls noise luminance suppression')], - 'satpercent': ['s', dict(type=float, default=0.15, help='Percentile value for saturation clipping (default: 0.15).')], - 'colorsatfac': ['c', dict(type=float, default=2, help='Factor for color saturation (default: 2).')], +TRILOGY_IMAGE_ARGS = { + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS cloud password.")], + "red": [ + "R", + dict( + metavar="RED_BAND", + default="I", + nargs="+", + help="Band(s) for the red channel (default: I R G).", + ), + ], + "green": [ + "G", + dict( + metavar="GREEN_BAND", + default="R", + nargs="+", + help="Band(s) for the green channel (default: R G B).", + ), + ], + "blue": [ + "B", + dict( + metavar="BLUE_BAND", + default="G", + nargs="+", + help="Band(s) for the blue channel (default: G B I).", + ), + ], + "field_name": [ + "F", + dict(metavar="FIELD", help='Field identifier, e.g., "SPLUS-n01s10"'), + ], + "size_unit": [ + "u", + dict( + metavar="SIZE_UNIT", + choices=["arcsec", "arcmin", "pixels"], + default="pix", + help="Unit of the size parameter (arcsec, arcmin, pixels).", + ), + ], + "data_release": ["D", dict(default="DR4", help="S-PLUS data release version.")], + "outfile": [ + "o", + dict( + default="trilogy_image.png", + help="Output filename for the Lupton RGB image.", + ), + ], + "noiselum": [ + "n", + dict(type=float, default=0.15, help="Controls noise luminance suppression"), + ], + "satpercent": [ + "s", + dict( + type=float, + default=0.15, + help="Percentile value for saturation clipping (default: 0.15).", + ), + ], + "colorsatfac": [ + "c", + dict(type=float, default=2, help="Factor for color saturation (default: 2)."), + ], # positional arguments - 'ra': ['pos', dict(metavar='RA', help="Object's right ascension")], - 'dec': ['pos', dict(metavar='DEC', help="Object's declination")], - 'size': ['pos', dict(metavar='SIZE', help='Size of the stamp in size_unit (defaults to pixels) (see --size_unit).')], + "ra": ["pos", dict(metavar="RA", help="Object's right ascension")], + "dec": ["pos", dict(metavar="DEC", help="Object's declination")], + "size": [ + "pos", + dict( + metavar="SIZE", + help="Size of the stamp in size_unit (defaults to pixels) (see --size_unit).", + ), + ], } -TRILOGY_IMAGE_PROG_DESC = f''' +TRILOGY_IMAGE_PROG_DESC = f""" Create a Trilogy RGB composite (multi-filter blend). -''' +""" + def trilogy_image(): - ''' + """ Script function for downloading a S-PLUS Trilogy RGB image. Raises @@ -195,47 +368,93 @@ def trilogy_image(): Returns ------- None - ''' - parser = create_parser(args_dict=TRILOGY_IMAGE_ARGS, program_description=TRILOGY_IMAGE_PROG_DESC) + """ + parser = create_parser( + args_dict=TRILOGY_IMAGE_ARGS, program_description=TRILOGY_IMAGE_PROG_DESC + ) args = parser.parse_args(args=sys.argv[1:]) conn = Core(args.username, args.password) - print_level(f'Downloading Trilogy RGB image at RA: {args.ra}, DEC: {args.dec}', 1, args.verbose) - conn.trilogy_image(args.ra, args.dec, args.size, R=args.red, G=args.green, B=args.blue, - field_name=args.field_name, size_unit=args.size_unit, - data_release=args.data_release, outfile=args.outfile, - noiselum=args.noiselum, satpercent=args.satpercent, - colorsatfac=args.colorsatfac, + print_level( + f"Downloading Trilogy RGB image at RA: {args.ra}, DEC: {args.dec}", + 1, + args.verbose, ) + conn.trilogy_image( + args.ra, + args.dec, + args.size, + R=args.red, + G=args.green, + B=args.blue, + field_name=args.field_name, + size_unit=args.size_unit, + data_release=args.data_release, + outfile=args.outfile, + noiselum=args.noiselum, + satpercent=args.satpercent, + colorsatfac=args.colorsatfac, + ) + CALIBRATED_STAMP_ARGS = { # optional arguments - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS cloud password.')], - 'weight': ['w', dict(action='store_true', default=False, help='Download the weight map associated to the stamp.')], - 'field_name': ['F', dict(metavar='FIELD', help='Field identifier, e.g., "SPLUS-n01s10"')], - 'size_unit': ['u', dict(metavar='SIZE_UNIT', choices=['arcsec', 'arcmin', 'pixels'], default='pix', help='Unit of the size parameter (arcsec, arcmin, pixels).')], - 'data_release': ['D', dict(default='DR4', help='S-PLUS data release version.')], - 'outfile': ['o', dict(default='stamp.fits.fz', help='Output filename for the stamp FITS image.')], - + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS cloud password.")], + "weight": [ + "w", + dict( + action="store_true", + default=False, + help="Download the weight map associated to the stamp.", + ), + ], + "field_name": [ + "F", + dict(metavar="FIELD", help='Field identifier, e.g., "SPLUS-n01s10"'), + ], + "size_unit": [ + "u", + dict( + metavar="SIZE_UNIT", + choices=["arcsec", "arcmin", "pixels"], + default="pix", + help="Unit of the size parameter (arcsec, arcmin, pixels).", + ), + ], + "data_release": ["D", dict(default="DR4", help="S-PLUS data release version.")], + "outfile": [ + "o", + dict(default="stamp.fits.fz", help="Output filename for the stamp FITS image."), + ], # positional arguments - 'ra': ['pos', dict(metavar='RA', help="Object's right ascension")], - 'dec': ['pos', dict(metavar='DEC', help="Object's declination")], - 'size': ['pos', dict(metavar='SIZE', help='Size of the stamp in size_unit (defaults to pixels) (see --size_unit).')], - 'band': ['pos', dict(metavar='BAND', help='S-PLUS band, e.g., "R", "I", "F660", "U".')], + "ra": ["pos", dict(metavar="RA", help="Object's right ascension")], + "dec": ["pos", dict(metavar="DEC", help="Object's declination")], + "size": [ + "pos", + dict( + metavar="SIZE", + help="Size of the stamp in size_unit (defaults to pixels) (see --size_unit).", + ), + ], + "band": [ + "pos", + dict(metavar="BAND", help='S-PLUS band, e.g., "R", "I", "F660", "U".'), + ], } -CALIBRATED_STAMP_PROG_DESC = f''' +CALIBRATED_STAMP_PROG_DESC = f""" Create a stamp and return a photometrically calibrated PrimaryHDU. This computes a cutout via `stamp(...)`, then loads the appropriate DR6+ per-field zero-point model and applies spatially varying calibration. -''' +""" + def calibrated_stamp(): - ''' + """ Script function for downloading a S-PLUS calibrated stamp FITS image. Raises @@ -246,13 +465,26 @@ def calibrated_stamp(): Returns ------- None - ''' - parser = create_parser(args_dict=CALIBRATED_STAMP_ARGS, program_description=CALIBRATED_STAMP_PROG_DESC) + """ + parser = create_parser( + args_dict=CALIBRATED_STAMP_ARGS, program_description=CALIBRATED_STAMP_PROG_DESC + ) args = parser.parse_args(args=sys.argv[1:]) conn = Core(args.username, args.password) - print_level(f'Downloading calibrated stamp at RA: {args.ra}, DEC: {args.dec} - band {args.band}', 1, args.verbose) - conn.calibrated_stamp(args.ra, args.dec, args.size, args.band, - weight=args.weight, field_name=args.field_name, - size_unit=args.size_unit, data_release=args.data_release, outfile=args.outfile, + print_level( + f"Downloading calibrated stamp at RA: {args.ra}, DEC: {args.dec} - band {args.band}", + 1, + args.verbose, + ) + conn.calibrated_stamp( + args.ra, + args.dec, + args.size, + args.band, + weight=args.weight, + field_name=args.field_name, + size_unit=args.size_unit, + data_release=args.data_release, + outfile=args.outfile, ) diff --git a/splusdata/scubes/__init__.py b/splusdata/scubes/__init__.py index b2fadd3..6600caf 100644 --- a/splusdata/scubes/__init__.py +++ b/splusdata/scubes/__init__.py @@ -1 +1 @@ -from splusdata.core import Core \ No newline at end of file +from splusdata.core import Core diff --git a/splusdata/scubes/core.py b/splusdata/scubes/core.py index c54c66c..18ee410 100644 --- a/splusdata/scubes/core.py +++ b/splusdata/scubes/core.py @@ -16,35 +16,45 @@ from splusdata.vars import BANDS, BANDWAVEINFO, get_band_info from splusdata.features.io import print_level, convert_coord_to_degrees -__scubes_author__ = 'Eduardo A. D. Lacerda ' -__scubes_version__ = '0.1.idr6-beta' +__scubes_author__ = "Eduardo A. D. Lacerda " +__scubes_version__ = "0.1.idr6-beta" + def _getval_array(pathlist, key, ext): return np.array([fits.getval(img, key, ext) for img in pathlist]) + def _getdata_array(pathlist, ext): return np.array([fits.getdata(img, ext=ext) for img in pathlist]) + def _getheader_array(pathlist, ext): return np.array([fits.getheader(img, ext=ext) for img in pathlist]) + def _getval_array_mem(hdull, key, ext): return np.array([hdul[ext].header.get(key) for hdul in hdull]) + def _getdata_array_mem(hdull, ext): return np.array([hdul[ext].data for hdul in hdull]) + def _getheader_array_mem(hdull, ext): return np.array([hdul[ext].header for hdul in hdull]) + def _get_band_info_array(prop): return np.array([get_band_info(band)[prop] for band in BANDS]) + class SCubes: - def __init__(self, ra, dec, field, size=None, username=None, password=None, verbose=1): + def __init__( + self, ra, dec, field, size=None, username=None, password=None, verbose=1 + ): # ignore some warnings without verbosity - warnings.simplefilter('ignore', category=VerifyWarning) - warnings.simplefilter('ignore', category=FITSFixedWarning) + warnings.simplefilter("ignore", category=VerifyWarning) + warnings.simplefilter("ignore", category=FITSFixedWarning) self.conn = Core(username, password, verbose=verbose) self.field = field @@ -59,68 +69,87 @@ def _stamp_config(self, ra, dec, size): self.size = size def _getval(self, obj, key, ext): - return _getval_array_mem(obj, key, ext) if self.mem else _getval_array(obj, key, ext) + return ( + _getval_array_mem(obj, key, ext) + if self.mem + else _getval_array(obj, key, ext) + ) def _getdata(self, obj, ext): return _getdata_array_mem(obj, ext) if self.mem else _getdata_array(obj, ext) def _getheader(self, obj, ext): - return _getheader_array_mem(obj, ext) if self.mem else _getheader_array(obj, ext) + return ( + _getheader_array_mem(obj, ext) if self.mem else _getheader_array(obj, ext) + ) def _download_calibrated_stamps(self, objname, force=False): images = [] wimages = [] _kw_args = dict(ra=self.ra, dec=self.dec, size=self.size, field_name=self.field) - for b in tqdm(BANDS, desc=f'{objname} @ {self.field} - Downloading', leave=True, position=0): - b = b.upper().replace('J0', 'F') + for b in tqdm( + BANDS, + desc=f"{objname} @ {self.field} - Downloading", + leave=True, + position=0, + ): + b = b.upper().replace("J0", "F") kw_args = _kw_args.copy() - kw_args.update(band=b, weight=False) + kw_args.update(band=b, weight=False) if self.mem: x = self.conn.calibrated_stamp(**kw_args) images.append(x) else: - filename = f'{objname}_{self.field}_{b}_{self.size}x{self.size}_swp.fits.fz' + filename = ( + f"{objname}_{self.field}_{b}_{self.size}x{self.size}_swp.fits.fz" + ) kw_args.update(outfile=join(self.outpath, filename)) _ = self.conn.calibrated_stamp(**kw_args) - images.append(kw_args['outfile']) + images.append(kw_args["outfile"]) # wei - kw_args['weight'] = True + kw_args["weight"] = True if self.mem: wimages.append(self.conn.calibrated_stamp(**kw_args)) else: - kw_args['outfile'] = join(self.outpath, filename.replace('swp', 'swpweight')) + kw_args["outfile"] = join( + self.outpath, filename.replace("swp", "swpweight") + ) _ = self.conn.calibrated_stamp(**kw_args) - wimages.append(kw_args['outfile']) + wimages.append(kw_args["outfile"]) self.images = images self.wimages = wimages def _photospectra(self, flam_scale=None, ext=1): flam_scale = 1e-19 if flam_scale is None else flam_scale _c = const.c - scale = (1/flam_scale) + scale = 1 / flam_scale Jy2fnu = 3631e-23 - self.wl__b = _get_band_info_array('pivot_wave')*u.Angstrom + self.wl__b = _get_band_info_array("pivot_wave") * u.Angstrom self.flam_unit = u.erg / u.s / u.cm / u.cm / u.AA self.fnu_unit = u.erg / u.s / u.cm / u.cm / u.Hz # flux calib_data__byx = self._getdata(self.images, ext) - fnu__byx = calib_data__byx*Jy2fnu*self.fnu_unit - flam__byx = scale*(fnu__byx*_c/self.wl__b[:, None, None]**2).to(self.flam_unit) + fnu__byx = calib_data__byx * Jy2fnu * self.fnu_unit + flam__byx = scale * (fnu__byx * _c / self.wl__b[:, None, None] ** 2).to( + self.flam_unit + ) # error in flux zp_factor__byx = self._getdata(self.images, 2) - absdata__byx = np.abs(calib_data__byx/zp_factor__byx) - gain__b = self._getval(self.images, 'GAIN', ext) + absdata__byx = np.abs(calib_data__byx / zp_factor__byx) + gain__b = self._getval(self.images, "GAIN", ext) gain__byx = gain__b[:, None, None] weidata__byx = self._getdata(self.wimages, ext) absweidata__byx = np.abs(self._getdata(self.wimages, ext)) - dataerr__byx = np.sqrt(1/absweidata__byx + absdata__byx/gain__byx) - f0__byx = zp_factor__byx*Jy2fnu - efnu__byx = dataerr__byx*f0__byx*self.fnu_unit - eflam__byx = scale*(efnu__byx*_c/self.wl__b[:, None, None]**2).to(self.flam_unit) + dataerr__byx = np.sqrt(1 / absweidata__byx + absdata__byx / gain__byx) + f0__byx = zp_factor__byx * Jy2fnu + efnu__byx = dataerr__byx * f0__byx * self.fnu_unit + eflam__byx = scale * (efnu__byx * _c / self.wl__b[:, None, None] ** 2).to( + self.flam_unit + ) self.flam__byx = flam__byx self.eflam__byx = eflam__byx @@ -128,7 +157,7 @@ def _photospectra(self, flam_scale=None, ext=1): self.absweidata__byx = absweidata__byx def _stamp_WCS_to_cube_header(self, header): - ''' + """ Convert WCS information from stamp to cube header. Parameters @@ -140,33 +169,36 @@ def _stamp_WCS_to_cube_header(self, header): ------- :class:`~astropy.io.fits.Header` Cube header with updated WCS information. - ''' + """ w = WCS(header) nw = WCS(naxis=3) nw.wcs.cdelt = [w.wcs.cdelt[0], w.wcs.cdelt[1], 1] nw.wcs.crval = [w.wcs.crval[0], w.wcs.crval[1], 0] nw.wcs.crpix = [w.wcs.crpix[0], w.wcs.crpix[1], 0] - nw.wcs.ctype = [w.wcs.ctype[0], w.wcs.ctype[1], ''] + nw.wcs.ctype = [w.wcs.ctype[0], w.wcs.ctype[1], ""] try: nw.wcs.pc[:2, :2] = w.wcs.get_pc() except: pass return nw.to_header() - + def _weights_mask_hdu(self): # WEIGHTS MASK HDU w__byx = self.weidata__byx wmask__byx = np.where(w__byx < 0, 1, 0) wmask__yx = wmask__byx.sum(axis=0) wmask_hdu = fits.ImageHDU(wmask__yx) - wmask_hdu.header['EXTNAME'] = ('WEIMASK', 'Sum of negative weight pixels (from 1 to 12)') + wmask_hdu.header["EXTNAME"] = ( + "WEIMASK", + "Sum of negative weight pixels (from 1 to 12)", + ) return wmask_hdu - + def _metadata_hdu(self, ext=1): # METADATA HDU tab = [BANDS] - tab.append(_get_band_info_array('central_wave')) - tab.append(_get_band_info_array('pivot_wave')) + tab.append(_get_band_info_array("central_wave")) + tab.append(_get_band_info_array("pivot_wave")) # PSFFWHM psffwhm__b = [] for img in self.images: @@ -174,63 +206,82 @@ def _metadata_hdu(self, ext=1): hdr = img[ext].header else: hdr = fits.getheader(img, ext=ext) - key = [k for k in hdr.keys() if 'FWHMMEAN' in k] + key = [k for k in hdr.keys() if "FWHMMEAN" in k] if len(key) == 1: psffwhm__b.append(hdr.get(key[0])) tab.append(psffwhm__b) - tab = Table(tab, names=['FILTER', 'CENTWAVE', 'PIVOTWAVE', 'PSFFWHM']) - return fits.BinTableHDU(tab, name='METADATA') + tab = Table(tab, names=["FILTER", "CENTWAVE", "PIVOTWAVE", "PSFFWHM"]) + return fits.BinTableHDU(tab, name="METADATA") def _create_cube_hdulist(self, objname, ext=1, overwrite=False): cube_prim_hdu = fits.PrimaryHDU() - cube_prim_hdu.header['TILE'] = self.field - cube_prim_hdu.header['OBJECT'] = objname - cube_prim_hdu.header['SIZE'] = (self.size, 'Side of the stamp in pixels') - cube_prim_hdu.header['RA'] = self.ra - cube_prim_hdu.header['DEC'] = self.dec - hdr = self.images[0][ext].header if self.mem else fits.getheader(self.images[0], ext=ext) + cube_prim_hdu.header["TILE"] = self.field + cube_prim_hdu.header["OBJECT"] = objname + cube_prim_hdu.header["SIZE"] = (self.size, "Side of the stamp in pixels") + cube_prim_hdu.header["RA"] = self.ra + cube_prim_hdu.header["DEC"] = self.dec + hdr = ( + self.images[0][ext].header + if self.mem + else fits.getheader(self.images[0], ext=ext) + ) cube_prim_hdu.header.update(self._stamp_WCS_to_cube_header(hdr)) - for key in ['X0TILE', 'X01TILE', 'Y0TILE', 'Y01TILE']: + for key in ["X0TILE", "X01TILE", "Y0TILE", "Y01TILE"]: cube_prim_hdu.header[key] = hdr.get(key) # DATA HDU flam_hdu = fits.ImageHDU(self.flam__byx.value, cube_prim_hdu.header) - flam_hdu.header['EXTNAME'] = ('DATA', 'Name of the extension') + flam_hdu.header["EXTNAME"] = ("DATA", "Name of the extension") # ERRORS HDU eflam_hdu = fits.ImageHDU(self.eflam__byx.value, cube_prim_hdu.header) - eflam_hdu.header['EXTNAME'] = ('ERRORS', 'Name of the extension') + eflam_hdu.header["EXTNAME"] = ("ERRORS", "Name of the extension") hdul = [cube_prim_hdu, flam_hdu, eflam_hdu] # INFO TO HEADERS for hdu in hdul[1:]: - hdu.header['BSCALE'] = (self.flam_scale, 'Linear factor in scaling equation') - hdu.header['BZERO'] = (0, 'Zero point in scaling equation') - hdu.header['BUNIT'] = (f'{self.flam_unit}', 'Physical units of the array values') + hdu.header["BSCALE"] = ( + self.flam_scale, + "Linear factor in scaling equation", + ) + hdu.header["BZERO"] = (0, "Zero point in scaling equation") + hdu.header["BUNIT"] = ( + f"{self.flam_unit}", + "Physical units of the array values", + ) hdul.append(self._weights_mask_hdu()) hdul.append(self._metadata_hdu(ext)) - print_level(f'writting cube {self.cubepath}', 1, self.verbose) + print_level(f"writting cube {self.cubepath}", 1, self.verbose) fits.HDUList(hdul).writeto(self.cubepath, overwrite=overwrite) - print_level(f'Cube successfully created!', 1, self.verbose) + print_level(f"Cube successfully created!", 1, self.verbose) return fits.open(self.cubepath) - def create_cube(self, flam_scale=None, objname=None, outpath=None, force=False, data_ext=1, return_scube=False, force_mem=False): + def create_cube( + self, + flam_scale=None, + objname=None, + outpath=None, + force=False, + data_ext=1, + return_scube=False, + force_mem=False, + ): self.flam_scale = 1e-19 if flam_scale is None else flam_scale - self.objname = 'myobj' if objname is None else objname - self.outpath = '.' if outpath is None else outpath + self.objname = "myobj" if objname is None else objname + self.outpath = "." if outpath is None else outpath mkcube = True - - try: + + try: makedirs(self.outpath) except FileExistsError: - print_level(f'{self.outpath}: directory already exists', 1, self.verbose) + print_level(f"{self.outpath}: directory already exists", 1, self.verbose) - self.cubepath = join(self.outpath, f'{self.objname}_cube.fits') + self.cubepath = join(self.outpath, f"{self.objname}_cube.fits") if exists(self.cubepath) and not force: mkcube = False - print_level(f'{self.cubepath}: cube already exists', 1, self.verbose) - + print_level(f"{self.cubepath}: cube already exists", 1, self.verbose) + if not self.mem and force_mem: self.mem = True - + cube = None if mkcube: self._download_calibrated_stamps(objname, force=force) @@ -240,4 +291,4 @@ def create_cube(self, flam_scale=None, objname=None, outpath=None, force=False, if return_scube: return read_scube(cube) - return cube \ No newline at end of file + return cube diff --git a/splusdata/scubes/read.py b/splusdata/scubes/read.py index ffc481c..ae0c780 100644 --- a/splusdata/scubes/read.py +++ b/splusdata/scubes/read.py @@ -6,6 +6,7 @@ from astropy.coordinates import SkyCoord from astropy.visualization import make_lupton_rgb + class tupperware_none(Namespace): def __init__(self): pass @@ -13,7 +14,8 @@ def __init__(self): def __getattr__(self, attr): r = self.__dict__.get(attr, None) return r - + + # create filters indexes list def _parse_filters_tuple(f_tup, filters): if isinstance(f_tup, list): @@ -32,11 +34,18 @@ def _parse_filters_tuple(f_tup, filters): i_f.append(f_tup) return i_f -def make_RGB_tom(flux__byx, - rgb=(7, 5, 9), rgb_f=(1, 1, 1), - pminmax=(5, 95), im_max=255, - # astropy.visualization.make_lupton_rgb() input vars - minimum=(0, 0, 0), Q=0, stretch=10): + +def make_RGB_tom( + flux__byx, + rgb=(7, 5, 9), + rgb_f=(1, 1, 1), + pminmax=(5, 95), + im_max=255, + # astropy.visualization.make_lupton_rgb() input vars + minimum=(0, 0, 0), + Q=0, + stretch=10, +): # get filters index(es) pmin, pmax = pminmax RGB = [] @@ -54,15 +63,18 @@ def make_RGB_tom(flux__byx, # percentiles Cmin, Cmax = np.nanpercentile(C, pmin), np.nanpercentile(C, pmax) # calc color intensities - RGB.append(im_max*(C - Cmin)/(Cmax - Cmin)) + RGB.append(im_max * (C - Cmin) / (Cmax - Cmin)) R, G, B = RGB # filters factors fR, fG, fB = rgb_f # make RGB image - return make_lupton_rgb(fR*R, fG*G, fB*B, Q=Q, minimum=minimum, stretch=stretch) + return make_lupton_rgb( + fR * R, fG * G, fB * B, Q=Q, minimum=minimum, stretch=stretch + ) + def get_distance(x, y, x0, y0, pa=0.0, ba=1.0): - ''' + """ Return an image (:class:`numpy.ndarray`) of the distance from the center ``(x0, y0)`` in pixels, assuming a projected disk. @@ -93,25 +105,26 @@ def get_distance(x, y, x0, y0, pa=0.0, ba=1.0): ------- pixel_distance : array Pixel distances. - ''' + """ y = np.asarray(y) - y0 x = np.asarray(x) - x0 x2 = x**2 y2 = y**2 - xy = x*y + xy = x * y - a_b = 1/ba + a_b = 1 / ba cos_th = np.cos(pa) sin_th = np.sin(pa) - A1 = cos_th**2 + a_b**2*sin_th**2 - A2 = -2*cos_th*sin_th*(a_b**2 - 1) - A3 = sin_th**2 + a_b**2*cos_th* 2 + A1 = cos_th**2 + a_b**2 * sin_th**2 + A2 = -2 * cos_th * sin_th * (a_b**2 - 1) + A3 = sin_th**2 + a_b**2 * cos_th * 2 + + return np.sqrt(A1 * x2 + A2 * xy + A3 * y2) - return np.sqrt(A1*x2 + A2*xy + A3*y2) def get_image_distance(shape, x0, y0, pa=0.0, ba=1.0): - ''' + """ Return an image (:class:`numpy.ndarray`) of the distance from the center ``(x0, y0)`` in pixels, assuming a projected disk. @@ -144,12 +157,24 @@ def get_image_distance(shape, x0, y0, pa=0.0, ba=1.0): -------- :func:`get_distance` - ''' + """ y, x = np.indices(shape) return get_distance(x, y, x0, y0, pa, ba) -def radial_profile(prop, bin_r, x0, y0, pa=0.0, ba=1.0, rad_scale=1.0, mask=None, mode='mean', return_npts=False): - ''' + +def radial_profile( + prop, + bin_r, + x0, + y0, + pa=0.0, + ba=1.0, + rad_scale=1.0, + mask=None, + mode="mean", + return_npts=False, +): + """ Calculate the radial profile of an N-D image. Parameters @@ -211,10 +236,13 @@ def radial_profile(prop, bin_r, x0, y0, pa=0.0, ba=1.0, rad_scale=1.0, mask=None See also -------- :func:`get_image_distance` - ''' + """ + def red(func, x, fill_value): - if x.size == 0: return fill_value, fill_value - if x.ndim == 1: return func(x), len(x) + if x.size == 0: + return fill_value, fill_value + if x.ndim == 1: + return func(x), len(x) return func(x, axis=-1), x.shape[-1] imshape = prop.shape[-2:] @@ -223,21 +251,21 @@ def red(func, x, fill_value): r__yx = get_image_distance(imshape, x0, y0, pa, ba) / rad_scale if mask is None: mask = np.ones(imshape, dtype=bool) - if mode == 'mean': + if mode == "mean": reduce_func = np.mean - elif mode == 'median': + elif mode == "median": reduce_func = np.median - elif mode == 'sum': + elif mode == "sum": reduce_func = np.sum - elif mode == 'var': + elif mode == "var": reduce_func = np.var - elif mode == 'std': + elif mode == "std": reduce_func = np.std else: - raise ValueError('Invalid mode: %s' % mode) + raise ValueError("Invalid mode: %s" % mode) if isinstance(prop, np.ma.MaskedArray): - n_bad = prop.mask.astype('int') + n_bad = prop.mask.astype("int") max_bad = 1.0 while n_bad.ndim > 2: max_bad *= n_bad.shape[0] @@ -254,94 +282,107 @@ def red(func, x, fill_value): if mask.any(): dist_flat = r__yx[mask] dist_idx = np.digitize(dist_flat, bin_r) - prop_flat = prop[...,mask] + prop_flat = prop[..., mask] for i in range(0, nbins): - prop_profile[..., i], npts[i] = red(reduce_func, prop_flat[..., dist_idx == i+1], reduce_fill_value) + prop_profile[..., i], npts[i] = red( + reduce_func, prop_flat[..., dist_idx == i + 1], reduce_fill_value + ) if return_npts: return prop_profile, npts return prop_profile + class read_scube: def __init__(self, data): self._read_data(data) self._init() - + def _read_data(self, data): self._hdulist = fits.open(data) if isinstance(data, str) else data def _init_wcs(self): - ''' + """ Initializes the World Coordinate System (WCS) from the data header. - ''' + """ self.wcs = WCS(self.data_header, naxis=2) def _init_centre(self): - ''' + """ Calculates the central coordinates of the object in pixel space using WCS. - ''' - self.cent_coord = SkyCoord(self.ra, self.dec, unit=('deg', 'deg'), frame='icrs') + """ + self.cent_coord = SkyCoord(self.ra, self.dec, unit=("deg", "deg"), frame="icrs") self.x0, self.y0 = self.wcs.world_to_pixel(self.cent_coord) self.i_x0 = int(self.x0) self.i_y0 = int(self.y0) def _mag_values(self): - ''' + """ Computes the magnitude per square arcsecond and corresponding errors from the flux values. - ''' - a = 1/(2.997925e18*3631.0e-23*self.pixscale**2) - x = a*(self.flux__byx*self.pivot_wave[:, np.newaxis, np.newaxis]**2) - self.mag_arcsec2__byx = -2.5*np.log10(x) - self.emag_arcsec2__byx = (2.5*np.log10(np.exp(1)))*self.eflux__byx/self.flux__byx + """ + a = 1 / (2.997925e18 * 3631.0e-23 * self.pixscale**2) + x = a * (self.flux__byx * self.pivot_wave[:, np.newaxis, np.newaxis] ** 2) + self.mag_arcsec2__byx = -2.5 * np.log10(x) + self.emag_arcsec2__byx = ( + (2.5 * np.log10(np.exp(1))) * self.eflux__byx / self.flux__byx + ) def _init(self): - ''' + """ Initializes the class by setting WCS, central coordinates, and magnitude values. Also calculates the pixel distance from the central coordinates. - ''' + """ self._init_wcs() self._init_centre() self._mag_values() self.pa, self.ba = 0, 1 self.set_ellipticity(self.pa, self.ba) - #self.pixel_distance__yx = get_image_distance(self.weimask__yx.shape, x0=self.x0, y0=self.y0, pa=self.pa, ba=self.ba) + # self.pixel_distance__yx = get_image_distance(self.weimask__yx.shape, x0=self.x0, y0=self.y0, pa=self.pa, ba=self.ba) def set_ellipticity(self, pa, ba): self.pa = pa self.ba = ba - self.pixel_distance__yx = get_image_distance(self.weimask__yx.shape, x0=self.x0, y0=self.y0, pa=self.pa, ba=self.ba) + self.pixel_distance__yx = get_image_distance( + self.weimask__yx.shape, x0=self.x0, y0=self.y0, pa=self.pa, ba=self.ba + ) def get_filter_i(self, filt): return self.filters.index(filt) def lRGB_image( - self, rgb=('r', 'g', 'i'), rgb_f=(1, 1, 1), - pminmax=(5, 95), im_max=255, + self, + rgb=("r", "g", "i"), + rgb_f=(1, 1, 1), + pminmax=(5, 95), + im_max=255, # astropy.visualization.make_lupton_rgb() input vars - minimum=(0, 0, 0), Q=0, stretch=3): - ''' + minimum=(0, 0, 0), + Q=0, + stretch=3, + ): + """ Creates an RGB image from the data cube using specified filters. Parameters ---------- rgb : tuple of str or tuple of int, optional Tuple specifying the filters to use for the red, green, and blue channels (default is ('rSDSS', 'gSDSS', 'iSDSS')). - + rgb_f : tuple of float, optional Scaling factors for the red, green, and blue channels (default is (1, 1, 1)). - + pminmax : tuple of int, optional Percentiles for scaling the RGB intensities (default is (5, 95)). - + im_max : int, optional Maximum intensity value for the RGB image (default is 255). - + minimum : tuple of float, optional Minimum values for scaling the RGB channels (default is (0, 0, 0)). - + Q : float, optional Parameter for controlling the contrast in the Lupton RGB scaling (default is 0). - + stretch : float, optional Stretch factor for enhancing the RGB intensities (default is 3). @@ -349,7 +390,7 @@ def lRGB_image( ------- np.ndarray 3D array representing the RGB image with shape (height, width, 3). - ''' + """ # check filters if len(rgb) != 3: return None @@ -370,102 +411,109 @@ def lRGB_image( i_rgb = [_parse_filters_tuple(x, self.filters) for x in rgb] return make_RGB_tom( - self.flux__byx, rgb=i_rgb, rgb_f=rgb_f, - pminmax=pminmax, im_max=im_max, - minimum=minimum, Q=Q, stretch=stretch + self.flux__byx, + rgb=i_rgb, + rgb_f=rgb_f, + pminmax=pminmax, + im_max=im_max, + minimum=minimum, + Q=Q, + stretch=stretch, ) @property def weimask__byx(self): - return np.broadcast_to(self.weimask__yx, (len(self.filters), self.size, self.size)) + return np.broadcast_to( + self.weimask__yx, (len(self.filters), self.size, self.size) + ) @property def primary_header(self): - return self._hdulist['PRIMARY'].header + return self._hdulist["PRIMARY"].header @property def data_header(self): - return self._hdulist['DATA'].header + return self._hdulist["DATA"].header @property def metadata(self): - return self._hdulist['METADATA'].data - + return self._hdulist["METADATA"].data + @property def filters(self): - return self.metadata['FILTER'].tolist() + return self.metadata["FILTER"].tolist() @property def central_wave(self): - return self.metadata['CENTRALWAVE'] + return self.metadata["CENTRALWAVE"] @property def pivot_wave(self): - return self.metadata['PIVOTWAVE'] + return self.metadata["PIVOTWAVE"] @property def psf_fwhm(self): - return self.metadata['PSFFWHM'] + return self.metadata["PSFFWHM"] @property def tile(self): - return self.primary_header.get('TILE', None) - + return self.primary_header.get("TILE", None) + @property def galaxy(self): - return self.primary_header.get('GALAXY', None) - + return self.primary_header.get("GALAXY", None) + @property def size(self): - return self.primary_header.get('SIZE', None) - + return self.primary_header.get("SIZE", None) + @property def ra(self): - return self.primary_header.get('RA', None) - + return self.primary_header.get("RA", None) + @property def dec(self): - return self.primary_header.get('DEC', None) + return self.primary_header.get("DEC", None) @property def x0tile(self): - return self.primary_header.get('X0TILE', None) + return self.primary_header.get("X0TILE", None) @property def y0tile(self): - return self.primary_header.get(['Y0TILE'], None) + return self.primary_header.get(["Y0TILE"], None) @property def pixscale(self): - return self.data_header.get('PIXSCALE', 0.55) + return self.data_header.get("PIXSCALE", 0.55) - @property + @property def weimask__yx(self): - return self._hdulist['WEIMASK'].data + return self._hdulist["WEIMASK"].data - @property + @property def flux__byx(self): - return self._hdulist['DATA'].data + return self._hdulist["DATA"].data - @property + @property def eflux__byx(self): - return self._hdulist['ERRORS'].data - + return self._hdulist["ERRORS"].data + @property def n_x(self): - return self.data_header.get('NAXIS1', None) - + return self.data_header.get("NAXIS1", None) + @property def n_y(self): - return self.data_header.get('NAXIS2', None) - + return self.data_header.get("NAXIS2", None) + @property def n_filters(self): - return self.data_header['NAXIS3'] - + return self.data_header["NAXIS3"] + @property def SN__byx(self): - return self.flux__byx/self.eflux__byx + return self.flux__byx / self.eflux__byx @property def mag__byx(self): @@ -473,4 +521,4 @@ def mag__byx(self): @property def emag__byx(self): - return self.emag_arcsec2__byx \ No newline at end of file + return self.emag_arcsec2__byx diff --git a/splusdata/scubes/scripts.py b/splusdata/scubes/scripts.py index c83877f..d547edf 100644 --- a/splusdata/scubes/scripts.py +++ b/splusdata/scubes/scripts.py @@ -8,58 +8,60 @@ from splusdata.scripts.args import create_parser from splusdata.scubes.core import __scubes_author__, __scubes_version__ + def ml2header_updheader(cube_filename, ml_table, force=False): - ''' - Updates a S-CUBES raw cube primary header with the masterlist + """ + Updates a S-CUBES raw cube primary header with the masterlist information. Parameters ---------- cube_filename : str Path to S-CUBES raw cube. - + ml_table : :class:`astropy.table.table.Table` Masterlist read using :meth:`astropy.io.ascii.read` - + force : bool, optional - Force the update the key value is the key is existent at the - S-CUBES header. - ''' - with fits.open(cube_filename, 'update') as hdul: - hdu = hdul['PRIMARY'] + Force the update the key value is the key is existent at the + S-CUBES header. + """ + with fits.open(cube_filename, "update") as hdul: + hdu = hdul["PRIMARY"] # SNAME CONTROL - sname = hdu.header.get('OBJECT', None) + sname = hdu.header.get("OBJECT", None) if sname is None: - print_level('header: missing SNAME information') + print_level("header: missing SNAME information") sys.exit(1) - if sname not in ml_table['SNAME']: - print_level(f'masterlist: {sname}: missing SNAME information') + if sname not in ml_table["SNAME"]: + print_level(f"masterlist: {sname}: missing SNAME information") sys.exit(1) - mlcut = ml_table[ml_table['SNAME'] == sname] + mlcut = ml_table[ml_table["SNAME"] == sname] for col in ml_table.colnames: v = mlcut[col][0] desc = None if v is np.ma.masked: v = None - if '__' in col: - col, desc = col.split('__') + if "__" in col: + col, desc = col.split("__") if not force and (col in hdu.header): - continue - if col == 'FIELD' or col == 'SNAME': continue - if col == 'SIZE': - col = 'SIZE_ML' - desc = 'SIZE masterlist' + if col == "FIELD" or col == "SNAME": + continue + if col == "SIZE": + col = "SIZE_ML" + desc = "SIZE masterlist" hdu.header.set(col, value=v, comment=desc) -SPLUS_MOTD_TOP = '┌─┐ ┌─┐┬ ┬┌┐ ┌─┐┌─┐ ' -SPLUS_MOTD_MID = '└─┐───│ │ │├┴┐├┤ └─┐ ' -SPLUS_MOTD_BOT = '└─┘ └─┘└─┘└─┘└─┘└─┘ ' -SPLUS_MOTD_SEP = '----------------------' -SCUBES_PROG_DESC = f''' +SPLUS_MOTD_TOP = "┌─┐ ┌─┐┬ ┬┌┐ ┌─┐┌─┐ " +SPLUS_MOTD_MID = "└─┐───│ │ │├┴┐├┤ └─┐ " +SPLUS_MOTD_BOT = "└─┘ └─┘└─┘└─┘└─┘└─┘ " +SPLUS_MOTD_SEP = "----------------------" + +SCUBES_PROG_DESC = f""" {SPLUS_MOTD_TOP} | Create S-PLUS galaxies data cubes, a.k.a. S-CUBES. {SPLUS_MOTD_MID} | S-CUBES is an organized FITS file with data, errors, {SPLUS_MOTD_BOT} | mask and metadata about some galaxy present on any @@ -77,31 +79,45 @@ def ml2header_updheader(cube_filename, ml_table, force=False): Note that *10h37m2.5s* is a totally different angle from *10:37:2.5* (*159.26 deg* and *10.62 deg* respectively). -''' +""" -size_help = 'Size of the cube in pixels. ' -size_help += 'If size is a odd number, the program ' -size_help += 'will choose the closest even integer.' +size_help = "Size of the cube in pixels. " +size_help += "If size is a odd number, the program " +size_help += "will choose the closest even integer." SCUBES_ARGS = { # optional arguments - 'force': ['f', dict(action='store_true', default=False, help='Force overwrite of existing files.')], - 'size': ['l', dict(default=500, type=int, help=size_help)], - 'workdir': ['w', dict(default=getcwd(), help='Working directory.')], - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS cloud password.')], - 'force_mem': ['F', dict(action='store_true', default=False, help='Force memory mapping of downloaded input files.')], - + "force": [ + "f", + dict( + action="store_true", + default=False, + help="Force overwrite of existing files.", + ), + ], + "size": ["l", dict(default=500, type=int, help=size_help)], + "workdir": ["w", dict(default=getcwd(), help="Working directory.")], + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS cloud password.")], + "force_mem": [ + "F", + dict( + action="store_true", + default=False, + help="Force memory mapping of downloaded input files.", + ), + ], # positional arguments - 'field': ['pos', dict(metavar='SPLUS_TILE', help='Name of the S-PLUS field')], - 'ra': ['pos', dict(metavar='RA', help="Object's right ascension")], - 'dec': ['pos', dict(metavar='DEC', help="Object's declination")], - 'object': ['pos', dict(metavar='OBJECT_NAME', help="Object's name")], + "field": ["pos", dict(metavar="SPLUS_TILE", help="Name of the S-PLUS field")], + "ra": ["pos", dict(metavar="RA", help="Object's right ascension")], + "dec": ["pos", dict(metavar="DEC", help="Object's declination")], + "object": ["pos", dict(metavar="OBJECT_NAME", help="Object's name")], } - + + def scubes_argparse(args): - ''' + """ A particular parser of the command-line arguments for `scubes` script. Parameters @@ -113,14 +129,15 @@ def scubes_argparse(args): ------- :class:`argparse.Namespace` Command-line arguments parsed. - ''' + """ # closest even - args.size = round(float(args.size)/2)*2 + args.size = round(float(args.size) / 2) * 2 return args + def scubes(): - ''' + """ Script function for creating S-PLUS galaxy data cubes (S-CUBES). Raises @@ -131,50 +148,91 @@ def scubes(): Returns ------- None - ''' + """ from splusdata.scubes.core import SCubes parser = create_parser(args_dict=SCUBES_ARGS, program_description=SCUBES_PROG_DESC) # ADD VERSION OPTION - parser.add_argument('--version', action='version', version='%(prog)s {version}'.format(version=__scubes_version__)) + parser.add_argument( + "--version", + action="version", + version="%(prog)s {version}".format(version=__scubes_version__), + ) args = scubes_argparse(parser.parse_args(args=sys.argv[1:])) SCubes( - ra=args.ra, dec=args.dec, field=args.field, - size=args.size, username=args.username, password=args.password, + ra=args.ra, + dec=args.dec, + field=args.field, + size=args.size, + username=args.username, + password=args.password, verbose=args.verbose, ).create_cube( - objname=args.object, - outpath=join(args.workdir, args.object), - force=args.force, force_mem=args.force_mem + objname=args.object, + outpath=join(args.workdir, args.object), + force=args.force, + force_mem=args.force_mem, ) -SCUBESML_PROG_DESC = f''' + +SCUBESML_PROG_DESC = f""" {SPLUS_MOTD_TOP} | scubesml script: {SPLUS_MOTD_MID} | Create S-PLUS galaxies data cubes, a.k.a. S-CUBES {SPLUS_MOTD_BOT} | using the masterlist information as input. {SPLUS_MOTD_SEP} + {__scubes_author__} -''' +""" SCUBESML_ARGS = { # optional arguments - 'force': ['f', dict(action='store_true', default=False, help='Force overwrite of existing files.')], - 'size_multiplicator': ['S', dict(default=10, type=float, help='Factor to multiply the SIZE__pix value of the masterlist to create the galaxy size. If size is a odd number, the program will choose the closest even integer.')], - 'min_size': ['m', dict(default=200, type=int, help='Minimal size of the cube in pixels. If size negative, uses the original value of size calculation.')], - 'workdir': ['w', dict(default=getcwd(), help='Working directory.')], - 'verbose': ['v', dict(action='count', default=0, help='Verbosity level.')], - 'username': ['U', dict(default=None, help='S-PLUS Cloud username.')], - 'password': ['P', dict(default=None, help='S-PLUS Cloud password.')], - 'sname': ['O', dict(default=None, metavar='OBJECT_SNAME', help="Object's masterlist SNAME")], - 'force_mem': ['F', dict(action='store_true', default=False, help='Force memory mapping of downloaded input files.')], - - 'masterlist': ['pos', dict(metavar='MASTERLIST', help='Path to masterlist file')] + "force": [ + "f", + dict( + action="store_true", + default=False, + help="Force overwrite of existing files.", + ), + ], + "size_multiplicator": [ + "S", + dict( + default=10, + type=float, + help="Factor to multiply the SIZE__pix value of the masterlist to create the galaxy size. If size is a odd number, the program will choose the closest even integer.", + ), + ], + "min_size": [ + "m", + dict( + default=200, + type=int, + help="Minimal size of the cube in pixels. If size negative, uses the original value of size calculation.", + ), + ], + "workdir": ["w", dict(default=getcwd(), help="Working directory.")], + "verbose": ["v", dict(action="count", default=0, help="Verbosity level.")], + "username": ["U", dict(default=None, help="S-PLUS Cloud username.")], + "password": ["P", dict(default=None, help="S-PLUS Cloud password.")], + "sname": [ + "O", + dict(default=None, metavar="OBJECT_SNAME", help="Object's masterlist SNAME"), + ], + "force_mem": [ + "F", + dict( + action="store_true", + default=False, + help="Force memory mapping of downloaded input files.", + ), + ], + "masterlist": ["pos", dict(metavar="MASTERLIST", help="Path to masterlist file")], } + def scubesml_argparse(args): - ''' + """ A particular parser of the command-line arguments for `scubesml` script. Parameters @@ -186,21 +244,22 @@ def scubesml_argparse(args): ------- :class:`argparse.Namespace` Command-line arguments parsed. - ''' + """ # closest even from astropy.io import ascii try: args.ml = ascii.read(args.masterlist) except: - print_level(f'{args.masterlist}: unable to read file') + print_level(f"{args.masterlist}: unable to read file") sys.exit(1) - + return args + def scubesml(): - ''' - Script for creating S-PLUS galaxy data cubes (S-CUBES) using the + """ + Script for creating S-PLUS galaxy data cubes (S-CUBES) using the masterlist for the input arguments. Raises @@ -211,33 +270,51 @@ def scubesml(): Returns ------- None - ''' + """ from splusdata.scubes.core import SCubes - - parser = create_parser(args_dict=SCUBESML_ARGS, program_description=SCUBESML_PROG_DESC) + + parser = create_parser( + args_dict=SCUBESML_ARGS, program_description=SCUBESML_PROG_DESC + ) # ADD VERSION OPTION - parser.add_argument('--version', action='version', version='%(prog)s {version}'.format(version=__scubes_version__)) + parser.add_argument( + "--version", + action="version", + version="%(prog)s {version}".format(version=__scubes_version__), + ) args = scubesml_argparse(parser.parse_args(args=sys.argv[1:])) - sname_list = args.ml['SNAME'] if args.sname is None else [args.sname] + sname_list = args.ml["SNAME"] if args.sname is None else [args.sname] for sname in sname_list: - mlcut = args.ml[args.ml['SNAME'] == sname] - ra = mlcut['RA__deg'][0] - dec = mlcut['DEC__deg'][0] - field = mlcut['FIELD'][0] - size_pix = max(round(args.size_multiplicator*float(mlcut['SIZE__pix'][0])/2)*2, args.min_size) - #print(size_pix) - #print(ra, dec, field, size_pix) + mlcut = args.ml[args.ml["SNAME"] == sname] + ra = mlcut["RA__deg"][0] + dec = mlcut["DEC__deg"][0] + field = mlcut["FIELD"][0] + size_pix = max( + round(args.size_multiplicator * float(mlcut["SIZE__pix"][0]) / 2) * 2, + args.min_size, + ) + # print(size_pix) + # print(ra, dec, field, size_pix) creator = SCubes( - ra=ra, dec=dec, field=field, - size=size_pix, - username=args.username, password=args.password, verbose=args.verbose + ra=ra, + dec=dec, + field=field, + size=size_pix, + username=args.username, + password=args.password, + verbose=args.verbose, + ) + # try: + creator.create_cube( + objname=sname, + outpath=args.workdir, + force=args.force, + force_mem=args.force_mem, ) - #try: - creator.create_cube(objname=sname, outpath=args.workdir, force=args.force, force_mem=args.force_mem) if creator.cubepath is not None: ml2header_updheader(creator.cubepath, args.ml) - #except Exception as e: - # print_level(f'{sname}: Cube creation failed: {e}') \ No newline at end of file + # except Exception as e: + # print_level(f'{sname}: Cube creation failed: {e}') diff --git a/splusdata/vacs/__init__.py b/splusdata/vacs/__init__.py index aa8c2b7..253fd2f 100644 --- a/splusdata/vacs/__init__.py +++ b/splusdata/vacs/__init__.py @@ -1,2 +1,2 @@ from splusdata.vacs.sqg import * -from splusdata.vacs.pdfs import * \ No newline at end of file +from splusdata.vacs.pdfs import * diff --git a/splusdata/vacs/pdfs.py b/splusdata/vacs/pdfs.py index d91ec6e..c878a1b 100644 --- a/splusdata/vacs/pdfs.py +++ b/splusdata/vacs/pdfs.py @@ -1,17 +1,29 @@ import numpy as np import pandas as pd + # Function that calculates PDFs def Calc_PDF(x, Weights, Means, STDs): - PDF = np.sum(Weights*(1/(STDs*np.sqrt(2*np.pi))) * np.exp((-1/2) * ((x[:,None]-Means)**2)/(STDs)**2), axis=1) - return PDF/np.trapz(PDF, x) + PDF = np.sum( + Weights + * (1 / (STDs * np.sqrt(2 * np.pi))) + * np.exp((-1 / 2) * ((x[:, None] - Means) ** 2) / (STDs) ** 2), + axis=1, + ) + return PDF / np.trapz(PDF, x) + def Calculate_PDFs(PDF_Catalogue, x=np.arange(0, 1.001, 0.001)): Final_PDFs = [] PDF_Catalogue = PDF_Catalogue.reset_index(drop=True) for i in range(len(PDF_Catalogue)): - Final_PDFs.append(Calc_PDF(x, np.fromstring(PDF_Catalogue['PDF_Weights'][i][0:], sep=','), - np.fromstring(PDF_Catalogue['PDF_Means'][i][0:], sep=','), - np.fromstring(PDF_Catalogue['PDF_STDs'][i][0:], sep=','))) + Final_PDFs.append( + Calc_PDF( + x, + np.fromstring(PDF_Catalogue["PDF_Weights"][i][0:], sep=","), + np.fromstring(PDF_Catalogue["PDF_Means"][i][0:], sep=","), + np.fromstring(PDF_Catalogue["PDF_STDs"][i][0:], sep=","), + ) + ) - return x, Final_PDFs \ No newline at end of file + return x, Final_PDFs diff --git a/splusdata/vacs/sqg.py b/splusdata/vacs/sqg.py index 8a80d6d..cce1835 100644 --- a/splusdata/vacs/sqg.py +++ b/splusdata/vacs/sqg.py @@ -6,47 +6,58 @@ from astropy.coordinates import SkyCoord import requests + class SQGClass: __author__ = "Lilianne Nakazono" __version__ = "0.1.0" _pos = ["RA", "DEC"] - _morph = ['FWHM_n', 'A', 'B', 'KRON_RADIUS'] - _feat = ['u_iso', - 'J0378_iso', - 'J0395_iso', - 'J0410_iso', - 'J0430_iso', - 'g_iso', - 'J0515_iso', - 'r_iso', - 'J0660_iso', - 'i_iso', - 'J0861_iso', - 'z_iso'] - - def __init__(self, model='auto', verbose=False): + _morph = ["FWHM_n", "A", "B", "KRON_RADIUS"] + _feat = [ + "u_iso", + "J0378_iso", + "J0395_iso", + "J0410_iso", + "J0430_iso", + "g_iso", + "J0515_iso", + "r_iso", + "J0660_iso", + "i_iso", + "J0861_iso", + "z_iso", + ] + + def __init__(self, model="auto", verbose=False): self.model = model self.verbose = verbose try: - if self.model == 'RF16': + if self.model == "RF16": if self.verbose: print("Loading model...") - self.model_splus = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S4M") + self.model_splus = requests.get( + "https://splus.cloud/files/models/iDR3n4_RF_12S4M" + ) self.model_splus = self.model_splus.content self.model_splus = pickle.loads(self.model_splus) - elif self.model == 'RF18': + elif self.model == "RF18": if self.verbose: print("Loading model...") - self.model_wise = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S2W4M") + self.model_wise = requests.get( + "https://splus.cloud/files/models/iDR3n4_RF_12S2W4M" + ) self.model_wise = self.model_wise.content self.model_wise = pickle.loads(self.model_wise) else: if self.verbose: print("Loading both models...") - self.model_splus = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S4M") + self.model_splus = requests.get( + "https://splus.cloud/files/models/iDR3n4_RF_12S4M" + ) self.model_splus = self.model_splus.content self.model_splus = pickle.loads(self.model_splus) - self.model_wise = requests.get("https://splus.cloud/files/models/iDR3n4_RF_12S2W4M") + self.model_wise = requests.get( + "https://splus.cloud/files/models/iDR3n4_RF_12S2W4M" + ) self.model_wise = self.model_wise.content self.model_wise = pickle.loads(self.model_wise) except: @@ -54,7 +65,7 @@ def __init__(self, model='auto', verbose=False): @staticmethod def get_wise(df, list_names): - ''' + """ Deletes objects without WISE counterpart (> 2 arcsecond) from data Keywords arguments: @@ -62,7 +73,7 @@ def get_wise(df, list_names): list_names -- names of the columns corresponding to W1 and W2 magnitudes/error/snr from WISE returns a copy of a subset from data with objects that have WISE counterpart - ''' + """ copy = df.copy(deep=True) copy = copy.query("d2d <= 2") @@ -71,21 +82,32 @@ def get_wise(df, list_names): return copy def irsa_query(self, data): - ''' + """ Query ALLWISE catalogue from IRSA database. Note that this is optimized for querying per fields (i.e., the dataframe input must correspond to a unique field from S-PLUS) Keywords arguments: data -- dataframe containing information of RA and DEC for each object return dataframe containing WISE sources - ''' + """ if self.verbose: - print("Note that irsa_query() is not yet optimized for large data inputs. Please consider running classify() per S-PLUS field if irsa_query()=True.") + print( + "Note that irsa_query() is not yet optimized for large data inputs. Please consider running classify() per S-PLUS field if irsa_query()=True." + ) ramin = np.min(data["RA"]) ramax = np.max(data["RA"]) decmin = np.min(data["DEC"]) decmax = np.max(data["DEC"]) - col = ['ra', 'dec', 'w1mpro', 'w2mpro', 'w1snr', 'w2snr', 'w1sigmpro', 'w2sigmpro'] + col = [ + "ra", + "dec", + "w1mpro", + "w2mpro", + "w1snr", + "w2snr", + "w1sigmpro", + "w2sigmpro", + ] if ramin < 2 and ramax > 358: ramax = np.max(data.query("RA<2").RA) ramin = np.min(data.query("RA>358").RA) @@ -113,15 +135,23 @@ def irsa_query(self, data): data.rename(columns={feat: feat + "_" + str(0)}, inplace=True) df_wise = df_wise.rename( - columns={"col_0": col[0], "col_1": col[1], "col_2": col[2], - "col_3": col[3], "col_4": col[4], "col_5": col[5], - "col_6": col[6], "col_7": col[7]}) + columns={ + "col_0": col[0], + "col_1": col[1], + "col_2": col[2], + "col_3": col[3], + "col_4": col[4], + "col_5": col[5], + "col_6": col[6], + "col_7": col[7], + } + ) return df_wise @staticmethod def crossmatch(data, df_wise): - ''' + """ Perform cross match between two catalogues based on RA and DEC data_ra -- RA from S-PLUS catalogue @@ -130,12 +160,12 @@ def crossmatch(data, df_wise): df_wise_dec -- DEC from ALLWISE catalogue returns dataframe containing all original columns from data cross-matched with the ALLWISE query catalogue - ''' + """ coo_splus = SkyCoord(data["RA"] * u.deg, data["DEC"] * u.deg) coo_wise = SkyCoord(df_wise["ra"] * u.deg, df_wise["dec"] * u.deg) idx, d2d, d3d = coo_splus.match_to_catalog_sky(coo_wise) - data['d2d'] = pd.Series(d2d.arcsec) + data["d2d"] = pd.Series(d2d.arcsec) data = pd.concat([data, df_wise.iloc[idx].reset_index()], axis=1) return data @@ -146,30 +176,42 @@ def check_columns(data, list_columns): for element in list_columns: data[element] except ValueError: - print("Please ensure that ", element, - "column is in data. Please provide a ALLWISE cross-matched data or set match_irsa == True. Use model == 'opt' if the provided data do not have any source with WISE counteerpart.") + print( + "Please ensure that ", + element, + "column is in data. Please provide a ALLWISE cross-matched data or set match_irsa == True. Use model == 'opt' if the provided data do not have any source with WISE counteerpart.", + ) return - def check_match_irsa(self, data, match_irsa): - ''' + def check_match_irsa(self, data, match_irsa): + """ Check statement of match_irsa and performs cross-match if True data -- dataframe containing S-PLUS information match_irsa -- If true, query ALLWISE catalogue and performs cross-match. If false, checks if all necessary columns are in data - ''' + """ if match_irsa == False: self.check_columns(data, self._feat_wise) self.check_columns(data, self._error_wise) else: try: - for element in self._pos: # check if RA and DEC are in data to perfom the cross-match + for ( + element + ) in ( + self._pos + ): # check if RA and DEC are in data to perfom the cross-match data[element] except ValueError: - print("Please ensure that ", element, - "column is in data, otherwise the cross-match with ALLWISE cannot be done.") + print( + "Please ensure that ", + element, + "column is in data, otherwise the cross-match with ALLWISE cannot be done.", + ) if self.verbose: - print("Querying ALLWISE catalogue via TapPlus(url='https://irsa.ipac.caltech.edu/TAP')...") + print( + "Querying ALLWISE catalogue via TapPlus(url='https://irsa.ipac.caltech.edu/TAP')..." + ) df_wise = self.irsa_query(data) @@ -178,13 +220,22 @@ def check_match_irsa(self, data, match_irsa): data = self.crossmatch(data, df_wise) return data - def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro': 'w1mpro', - 'w2mpro': 'w2mpro', 'w1snr': 'w1snr', - 'w2snr': 'w2snr', - 'w1sigmpro': 'w1sigmpro', - 'w2sigmpro': 'w2sigmpro'}, verbose=False): - - ''' + def classify( + self, + df, + return_prob=True, + match_irsa=False, + columns_wise={ + "w1mpro": "w1mpro", + "w2mpro": "w2mpro", + "w1snr": "w1snr", + "w2snr": "w2snr", + "w1sigmpro": "w1sigmpro", + "w2sigmpro": "w2sigmpro", + }, + verbose=False, + ): + """ Create classifications for sources with or without counterpart Keywords arguments: @@ -196,13 +247,17 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro the source has WISE counterpart, otherwise returns classification from model "RF16" (flagged as 1). match_irsa -- determines if cross-match with ALLWISE catalogue will be performed. If model == "RF16", match_irsa == False. returns a dataframe with classes - ''' + """ data = df.copy(deep=True) self._feat_wise = [columns_wise["w1mpro"], columns_wise["w2mpro"]] - self._error_wise = [columns_wise["w1snr"], columns_wise["w2snr"], - columns_wise["w1sigmpro"], columns_wise["w2sigmpro"]] + self._error_wise = [ + columns_wise["w1snr"], + columns_wise["w2snr"], + columns_wise["w1sigmpro"], + columns_wise["w2sigmpro"], + ] if self.model == "RF16": if match_irsa: @@ -218,17 +273,23 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro try: data[element] except: - raise (KeyError)("Please ensure that ", element, - "column is in data. Use splusdata.connect.query() to retrieve it.") + raise (KeyError)( + "Please ensure that ", + element, + "column is in data. Use splusdata.connect.query() to retrieve it.", + ) self.results = pd.DataFrame() if verbose: print( - "Note that this function assumes that the input data were previously extinction corrected. Starting classification... ") + "Note that this function assumes that the input data were previously extinction corrected. Starting classification... " + ) if self.model == "RF16": - y_pred = pd.DataFrame(self.model_splus.predict(data[self._morph + self._feat])) + y_pred = pd.DataFrame( + self.model_splus.predict(data[self._morph + self._feat]) + ) y_pred = y_pred.astype("float64") y_pred.index = data.index @@ -238,7 +299,9 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro if return_prob: if verbose: print("Calculating probabilities...") - prob_df = pd.DataFrame(self.model_splus.predict_proba(data[self._morph + self._feat])) + prob_df = pd.DataFrame( + self.model_splus.predict_proba(data[self._morph + self._feat]) + ) prob_df.index = data.index results = pd.concat([results, prob_df], axis=1) @@ -247,8 +310,11 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro elif self.model == "RF18": data = self.check_match_irsa(data, match_irsa) data_wise = self.get_wise(data, self._feat_wise + self._error_wise) - ypred_wise = pd.DataFrame(self.model_wise.predict( - data_wise[self._morph + self._feat + self._feat_wise])) # classify sources with WISE + ypred_wise = pd.DataFrame( + self.model_wise.predict( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) # classify sources with WISE ypred_wise.index = data_wise.index results = ypred_wise results.columns = ["CLASS"] @@ -258,7 +324,10 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro print("Calculating probabilities...") prob_wise_df = pd.DataFrame( - self.model_wise.predict_proba(data_wise[self._morph + self._feat + self._feat_wise])) + self.model_wise.predict_proba( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) prob_wise_df.index = data_wise.index results = pd.concat([results, prob_wise_df], axis=1) results.columns = ["CLASS", "PROB_QSO", "PROB_STAR", "PROB_GAL"] @@ -268,13 +337,18 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro data_wise = self.get_wise(data, self._feat_wise + self._error_wise) data_nowise = data.drop(data_wise.index) - y_pred = pd.DataFrame(self.model_splus.predict(data_nowise[self._morph + self._feat])) + y_pred = pd.DataFrame( + self.model_splus.predict(data_nowise[self._morph + self._feat]) + ) y_pred = y_pred.astype("float64") y_pred.index = data_nowise.index y_pred["model_flag"] = "1" # S-PLUS - ypred_wise = pd.DataFrame(self.model_wise.predict( - data_wise[self._morph + self._feat + self._feat_wise])) # classify sources with WISE + ypred_wise = pd.DataFrame( + self.model_wise.predict( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) # classify sources with WISE ypred_wise.index = data_wise.index ypred_wise["model_flag"] = "0" # S-PLUS + WISE results = pd.concat([y_pred, ypred_wise], axis=0) @@ -283,14 +357,27 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro if return_prob: if verbose: print("Calculating probabilities...") - prob_df = pd.DataFrame(self.model_splus.predict_proba(data_nowise[self._morph + self._feat])) + prob_df = pd.DataFrame( + self.model_splus.predict_proba( + data_nowise[self._morph + self._feat] + ) + ) prob_wise_df = pd.DataFrame( - self.model_wise.predict_proba(data_wise[self._morph + self._feat + self._feat_wise])) + self.model_wise.predict_proba( + data_wise[self._morph + self._feat + self._feat_wise] + ) + ) prob_df.index = data_nowise.index prob_wise_df.index = data_wise.index prob_results = pd.concat([prob_df, prob_wise_df], axis=0) results = pd.concat([results, prob_results], axis=1) - results.columns = ["CLASS", "model_flag", "PROB_QSO", "PROB_STAR", "PROB_GAL"] + results.columns = [ + "CLASS", + "model_flag", + "PROB_QSO", + "PROB_STAR", + "PROB_GAL", + ] else: raise (ValueError)("Parameter 'model' should be 'RF16', 'RF18' or 'auto'.") @@ -298,4 +385,4 @@ def classify(self, df, return_prob=True, match_irsa=False, columns_wise={'w1mpro if verbose: print("Finished process.") - return results.sort_index(axis=0) \ No newline at end of file + return results.sort_index(axis=0) diff --git a/splusdata/vars.py b/splusdata/vars.py index c823a04..986e7f9 100644 --- a/splusdata/vars.py +++ b/splusdata/vars.py @@ -1,66 +1,220 @@ -BANDS = ['u', 'J0378', 'J0395', 'J0410', 'J0430', 'g', 'J0515', 'r', 'J0660', 'i', 'J0861', 'z'] +BANDS = [ + "u", + "J0378", + "J0395", + "J0410", + "J0430", + "g", + "J0515", + "r", + "J0660", + "i", + "J0861", + "z", +] WAVELENGHTS = [3536, 3770, 3940, 4094, 4292, 4751, 5133, 6258, 6614, 7690, 8611, 8831] DR_POINTINGS = { - "dr4": - { - "link": "https://splus.cloud/files/documentation/iDR4/tabelas/iDR4_pointings.csv", - "ra_col": "RA", - "dec_col": "DEC", - "field_col": "Field" - }, - "dr5": - { - "link": "https://splus.cloud/files/documentation/idr5/auxiliary_tables/iDR5_pointings.csv", - "ra_col": "RA_d", - "dec_col": "DEC_d", - "field_col": "iDR5_Field_Name" - }, - "dr6": - { - "link": "https://splus.cloud/files/dr6_docs/dr6_list.csv", - "ra_col": "ra", - "dec_col": "dec", - "field_col": "field" - } + "dr4": { + "link": "https://splus.cloud/files/documentation/iDR4/tabelas/iDR4_pointings.csv", + "ra_col": "RA", + "dec_col": "DEC", + "field_col": "Field", + }, + "dr5": { + "link": "https://splus.cloud/files/documentation/idr5/auxiliary_tables/iDR5_pointings.csv", + "ra_col": "RA_d", + "dec_col": "DEC_d", + "field_col": "iDR5_Field_Name", + }, + "dr6": { + "link": "https://splus.cloud/files/dr6_docs/dr6_list.csv", + "ra_col": "ra", + "dec_col": "dec", + "field_col": "field", + }, } -BANDWAVEINFO_COLS = ['central_wave','delta_wave','trapz_wave','trapz_width','mean_wave','mean_width','mean_1_wave','mean_1_width','pivot_wave','alambda_av'] +BANDWAVEINFO_COLS = [ + "central_wave", + "delta_wave", + "trapz_wave", + "trapz_width", + "mean_wave", + "mean_width", + "mean_1_wave", + "mean_1_width", + "pivot_wave", + "alambda_av", +] BANDWAVEINFO = { - 'u': [3576.5900319003185,324.89284892848946,3542.1443522227532,322.8263652334996,3542.1443530585752,322.8263652334996,3541.969648229418,322.47994720424987,3533.2815060278176,1.6104999642584739], - 'J0378': [3770.667656676567,150.9900099000988,3774.011802885624,135.9577193603005,3774.0118028527427,135.9577193603005,3773.981717868756,135.7529151577093,3773.1649561852896,1.5176985654998583], - 'J0395': [3940.669006690067,102.79782797828011,3941.090848532596,100.77648141472922,3941.0908486234594,100.77648141472922,3941.074611723043,100.664883796743,3940.6981217186826,1.4589092149273093], - 'J0410': [4094.0795907959077,200.31290312903093,4096.244672289782,193.3458020035327,4096.244672355809,193.3458020035327,4096.2061617095505,193.11785361385836,4094.928007328532,1.4033758391668887], - 'J0430': [4292.020120201202,200.16240162401664,4293.375187582308,195.06906616375664,4293.375187489372,195.06906616375664,4293.3383586839045,194.82176337031825,4292.105790056625,1.333792079417452], - 'g': [4774.026040260403,1505.4610546105455,4821.096812428365,1312.4365137224963,4821.096812419468,1312.4365137224963,4821.067328442389,1312.1660034147683,4758.487858695381,1.1988078567002278], - 'J0515': [5132.820973209732,207.0569705697053,5134.215757880079,203.61492665730344,5134.215757789297,203.61492665730344,5134.203017754601,203.4775265629296,5133.132479753686,1.0983855043535193], - 'r': [6274.743347433474,1436.69456694567,6295.687061301507,1274.085080552066,6295.687061340775,1274.085080552066,6295.670469955384,1273.717585280557,6251.8309742925085,0.8639150819807985], - 'J0660': [6613.9931899319,147.28497284972855,6614.316749499244,146.66560010379152,6614.316749615694,146.66560010379152,6614.295914036835,146.52087279655584,6613.875560388883,0.7979103155197371], - 'i': [7702.49932499325,1506.845468454685,7709.955146600505,1438.1013394782774,7709.955151443986,1438.1013394782774,7709.81234140586,1437.2077440992825,7670.614459826697,0.6479238989570093], - 'J0861': [8611.481664816649,409.6853968539672,8609.86480549928,401.9084864881284,8609.866138424737,401.90954637204015,8609.836152859618,401.52828207867316,8607.254217015387,0.5391299762216387], - 'z': [8881.70071700717,1270.4971049710512,8985.814083997144,1308.4931400034807,8985.814681949101,1308.4931400034807,8986.541255819202,1307.6954289945693,8941.476066226343,0.5121720871113746] + "u": [ + 3576.5900319003185, + 324.89284892848946, + 3542.1443522227532, + 322.8263652334996, + 3542.1443530585752, + 322.8263652334996, + 3541.969648229418, + 322.47994720424987, + 3533.2815060278176, + 1.6104999642584739, + ], + "J0378": [ + 3770.667656676567, + 150.9900099000988, + 3774.011802885624, + 135.9577193603005, + 3774.0118028527427, + 135.9577193603005, + 3773.981717868756, + 135.7529151577093, + 3773.1649561852896, + 1.5176985654998583, + ], + "J0395": [ + 3940.669006690067, + 102.79782797828011, + 3941.090848532596, + 100.77648141472922, + 3941.0908486234594, + 100.77648141472922, + 3941.074611723043, + 100.664883796743, + 3940.6981217186826, + 1.4589092149273093, + ], + "J0410": [ + 4094.0795907959077, + 200.31290312903093, + 4096.244672289782, + 193.3458020035327, + 4096.244672355809, + 193.3458020035327, + 4096.2061617095505, + 193.11785361385836, + 4094.928007328532, + 1.4033758391668887, + ], + "J0430": [ + 4292.020120201202, + 200.16240162401664, + 4293.375187582308, + 195.06906616375664, + 4293.375187489372, + 195.06906616375664, + 4293.3383586839045, + 194.82176337031825, + 4292.105790056625, + 1.333792079417452, + ], + "g": [ + 4774.026040260403, + 1505.4610546105455, + 4821.096812428365, + 1312.4365137224963, + 4821.096812419468, + 1312.4365137224963, + 4821.067328442389, + 1312.1660034147683, + 4758.487858695381, + 1.1988078567002278, + ], + "J0515": [ + 5132.820973209732, + 207.0569705697053, + 5134.215757880079, + 203.61492665730344, + 5134.215757789297, + 203.61492665730344, + 5134.203017754601, + 203.4775265629296, + 5133.132479753686, + 1.0983855043535193, + ], + "r": [ + 6274.743347433474, + 1436.69456694567, + 6295.687061301507, + 1274.085080552066, + 6295.687061340775, + 1274.085080552066, + 6295.670469955384, + 1273.717585280557, + 6251.8309742925085, + 0.8639150819807985, + ], + "J0660": [ + 6613.9931899319, + 147.28497284972855, + 6614.316749499244, + 146.66560010379152, + 6614.316749615694, + 146.66560010379152, + 6614.295914036835, + 146.52087279655584, + 6613.875560388883, + 0.7979103155197371, + ], + "i": [ + 7702.49932499325, + 1506.845468454685, + 7709.955146600505, + 1438.1013394782774, + 7709.955151443986, + 1438.1013394782774, + 7709.81234140586, + 1437.2077440992825, + 7670.614459826697, + 0.6479238989570093, + ], + "J0861": [ + 8611.481664816649, + 409.6853968539672, + 8609.86480549928, + 401.9084864881284, + 8609.866138424737, + 401.90954637204015, + 8609.836152859618, + 401.52828207867316, + 8607.254217015387, + 0.5391299762216387, + ], + "z": [ + 8881.70071700717, + 1270.4971049710512, + 8985.814083997144, + 1308.4931400034807, + 8985.814681949101, + 1308.4931400034807, + 8986.541255819202, + 1307.6954289945693, + 8941.476066226343, + 0.5121720871113746, + ], } + def get_band_info(band): - ''' + """ Retreive T80-South filter information. Parameters ---------- band : str - The name of the T80-S filter. + The name of the T80-S filter. Returns ------- None or dict - If the information for the required band is found, + If the information for the required band is found, returns a dictionary, otherwise None. See Also -------- https://github.com/splus-collab/splus_filters for more info about T80-South filters. - ''' + """ retval = None band = band.lower() # solve for U or u, G or g, etc... @@ -68,4 +222,4 @@ def get_band_info(band): if len(key) == 1: key = key[0] retval = {k: v for (k, v) in zip(BANDWAVEINFO_COLS, BANDWAVEINFO[key])} - return retval \ No newline at end of file + return retval diff --git a/tests/test_core.py b/tests/test_core.py new file mode 100644 index 0000000..9cdc47b --- /dev/null +++ b/tests/test_core.py @@ -0,0 +1,80 @@ +import io + +from PIL import Image + +import pytest + +from splusdata.core import open_image, save_image, Core, SplusdataError + + +def make_png_bytes(size=(8, 8), color=(255, 0, 0)): + im = Image.new("RGB", size, color) + buf = io.BytesIO() + im.save(buf, format="PNG") + return buf.getvalue() + + +def test_open_image_and_save_image(tmp_path): + data = make_png_bytes((16, 16)) + im = open_image(data) + assert isinstance(im, Image.Image) + assert im.size == (16, 16) + + out = tmp_path / "out.png" + save_image(data, str(out)) + assert out.exists() + + +class FakeClient: + def __init__(self): + self._collections = [ + {"id": 1, "name": "dr4", "patterns": {"": ""}}, + {"id": 2, "name": "dr6", "patterns": {"": ""}}, + ] + + def get_collections(self): + return self._collections + + def list_files(self, collection_id, filter_str=None, filter_name=None, **kwargs): + # Return a few fake files with filenames that allow testing include/exclude + files = [ + {"id": 10, "filename": f"{filter_str}_a.fz", "file_type": "fz"}, + {"id": 11, "filename": f"{filter_str}_b.fz", "file_type": "fz"}, + {"id": 12, "filename": f"{filter_str}_c.txt", "file_type": "txt"}, + ] + return files + + def download_file(self, file_id, output_path=None, timeout=None, **kwargs): + # Return a tiny fits-like header bytes for tests (not used in these tests) + return b"FAKEBYTES" + + +def make_core_with_fake_client(): + # Create Core instance without calling __init__ to avoid prompts / network + core = Core.__new__(Core) + core.client = FakeClient() + core.collections = [] + core.verbose = 0 + return core + + +def test_get_collection_id_by_pattern_found(): + core = make_core_with_fake_client() + col = core.get_collection_id_by_pattern("dr4") + assert col["name"] == "dr4" + + +def test_get_collection_id_by_pattern_not_found_raises(): + core = make_core_with_fake_client() + with pytest.raises(SplusdataError): + core.get_collection_id_by_pattern("no-such-dr") + + +def test_get_file_metadata_basic(): + core = make_core_with_fake_client() + # Add a pattern to the dr4 collection so pattern check passes + core.client._collections[0]["patterns"] = {"": ""} + # Should return one of the fake files + res = core.get_file_metadata("SPLUS-0001", "R", pattern="", data_release="dr4") + assert isinstance(res, dict) + assert "filename" in res