Astropy: Python for Astronomy

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 1 of 15. You have probably multiplied a NumPy array by a constant. But if you accidentally mix meters and kilometers in a complex equation, Python will not warn you, and your final array will just contain the wrong numbers. Letting the library handle the physical dimensions automatically is exactly what this episode covers. We are looking at The Heart of Astropy: Units and Quantities. Astropy handles physical units through the Quantity object. A Quantity binds a number, or an array of numbers, to a physical unit. Because a Quantity is a direct subclass of a NumPy ndarray, you do not lose any performance or array functionality. You can slice, broadcast, and run standard NumPy math functions on them natively. They behave like normal arrays that just happen to know what physical property they represent. To create a Quantity, you import the units module from Astropy. Then, you simply multiply your number by the specific unit object. Three multiplied by units dot parsec creates a Quantity of three parsecs. Many developers assume this abstraction is slow, which is a common misconception. It is true that multiplying a massive NumPy array by a unit object creates a new copy of that array in memory. To avoid this overhead, Astropy overrides the bitwise left-shift operator, which looks like two less-than signs. Put your existing array on the left, the left-shift operator in the middle, and the Astropy unit on the right. This attaches the unit to the array in place without copying the underlying data. Here is the key insight. The library evaluates the mathematical relationships between units as you compute. Suppose you want to calculate the travel time of a spacecraft. Your distance is three point zero parsecs. Your speed is one hundred thirty kilometers per second. You create a Quantity for the distance and divide it by the Quantity for the speed. Astropy returns a new Quantity for the time. It handles the division of the values, and it calculates the resulting unit, returning parsec seconds per kilometer. That resulting unit is mathematically correct but not very helpful to read. You can change this using the to method available on every Quantity. You call the to method and pass in your target unit, such as units dot year. Astropy verifies that the source and target dimensions are both units of time, performs the internal conversion, and returns a new Quantity expressed in years. For quick standardizations, you can also access the dot SI or dot CGS properties on any Quantity to immediately convert it to the base units of those systems. When you perform math where the units cancel out entirely, like dividing meters by centimeters, Astropy returns a dimensionless Quantity. The unit still exists under the hood as an unscaled dimensionless type. If you need to pass this result to an external library that expects standard Python floats or plain NumPy arrays, you extract the raw number using the value attribute. Before you extract any final value, remember that you can always simplify messy, combined units down to their absolute fundamental physical dimensions by calling the decompose method. If you would like to support the show, you can search for DevStoriesEU on Patreon—it helps us keep making these. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 2 of 15. A standard sixty-four-bit float cannot represent a timestamp spanning billions of years without destroying its own precision. If you attempt to track a distant celestial event down to the microsecond, standard floating-point arithmetic simply runs out of bits. This exact hardware limitation is why we need Time Scales and Precision: The astropy.time Module. To maintain sub-nanosecond accuracy over cosmic history, this module completely avoids using a single number to represent time. Instead, it relies on a clever two-float architecture. Under the hood, every time value is split into two separate sixty-four-bit floats. The first float stores the integer portion of the Julian Date, which tracks the whole days. The second float stores the fractional part of the day. If Astropy added these together into one variable, the precision would instantly degrade. By keeping them apart, the library can perform highly precise time arithmetic locally, regardless of how deep into the past or future the date sits. When you create a Time object, you pass in the time value alongside two critical arguments: the format and the scale. Here is the key insight. Developers routinely confuse the visual format of a date with its physical time scale. They are entirely separate concepts. The format dictates solely how the data is read or written. Common formats include standard ISO strings, numerical Modified Julian Dates, and specialized representations like FITS headers. Format is just the structural wrapper around the numbers. The time scale, however, defines the physical reference frame of those numbers. Physical time scales account for Earth's rotation, leap seconds, and relativistic effects. Examples include Coordinated Universal Time, known as UTC, Terrestrial Time, or TT, and International Atomic Time, or TAI. A Modified Julian Date is not a physical scale, it is just a format. You can easily have a timestamp formatted as a Modified Julian Date, but physically anchored to the UTC scale. You must explicitly define both the format and the scale when feeding raw data into a Time object. Consider calculating the time difference between two spacecraft telemetry timestamps. You receive the first timestamp as an ISO format string from yesterday, and the second timestamp as an ISO string from today. First, you initialize two Time objects. You pass the raw string into the object, declare the format as ISO, and set the scale to UTC. You do this for both timestamps. Because Astropy supports vectorized operations, you could also pass a massive array of these strings into a single Time object, but the underlying initialization logic remains the same. To find the exact elapsed time, you subtract the earlier Time object from the later one. Astropy automatically executes the internal two-float arithmetic, returning a TimeDelta object that holds the exact duration. Next, you need to map that second telemetry timestamp to an absolute, uniform physical scale, removing the unpredictable jumps caused by UTC leap seconds. You take your second Time object and simply request its equivalent in International Atomic Time. You achieve this by calling the TAI attribute directly on the object. Astropy immediately triggers a recalculation of the internal floats, shifting the reference frame from UTC to TAI while preserving total nanosecond precision. The true power of the astropy.time module is that you never manually manage leap seconds or relativistic corrections; you define your formats and scales upfront, and the two-float architecture guarantees the math stays exact across the universe. I would like to take a moment to thank you for listening — it helps us a lot. Have a great one!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 3 of 15. Manually converting celestial coordinates from one reference system to another involves horrifying spherical trigonometry. A single matrix error can easily leave your telescope pointing at the wrong hemisphere. Navigating the Sky: The SkyCoord Class handles that math for you, turning complex spherical transformations into a simple attribute lookup. In Astropy, SkyCoord is the primary, high-level object for representing celestial positions. It acts as a smart wrapper over lower-level coordinate frames. When you observe a target, you need to record its position. Typically, you use Right Ascension and Declination. To build a SkyCoord, you pass in those two values, specify the units, like degrees or hours, and set the reference frame. If you do not specify a frame, it defaults to the International Celestial Reference System, or ICRS. You can create a single point in the sky easily, but this brings up a frequent trap. Users often parse a data file containing thousands of stars, create a new SkyCoord object for each star, and store them in a standard Python list. When they need to perform calculations, they loop over that list. This is horribly inefficient and bypasses the performance benefits of the library entirely. SkyCoord is built to handle arrays natively. Instead of creating thousands of individual objects, you pass a Numpy array of Right Ascension values and a corresponding array of Declination values into a single SkyCoord initialization. You get back exactly one SkyCoord object that contains your entire dataset. Consider a scenario where you are cross-matching a newly discovered set of one hundred radio transients recorded in ICRS with an existing optical catalog. By loading all one hundred transients into a single array-backed SkyCoord, the library passes the data down to heavily optimized C and Numpy routines. Operations that would take minutes in a Python loop execute in milliseconds. This array structure becomes especially powerful when you need to transform your data. Astronomical catalogs do not always agree on a coordinate system. Your radio transients might be in ICRS, but perhaps you need to analyze their distribution relative to the plane of the Milky Way. For that, you need the Galactic coordinate system. Here is the key insight. You do not need to look up rotation matrices or write conversion functions. Astropy contains a massive, built-in transformation graph that knows how to get from any supported frame to any other supported frame. Because SkyCoord wraps this graph, converting your entire array of one hundred transients takes exactly one step. You take your existing SkyCoord object and call its transform to method, providing the name of the destination frame. Better yet, for common frames, Astropy provides direct attribute access. If you have an ICRS SkyCoord object, you simply request the galactic attribute. The library automatically calculates the transformation and returns a brand new SkyCoord object. This new object contains the exact same physical locations in space, but represented as Galactic longitude and latitude. The true power of this class is that it completely detaches the data representing a celestial position from the underlying spherical geometry required to shift it across the universe of astronomical reference frames. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 4 of 15. You might know exactly where a star is tonight, but without proper motion and Earth's velocity corrections, your telescope will be pointing at empty space next decade. This is where we go Beyond RA and Dec: 3D Tracking and Velocities. By default, coordinates map a flat position on the sky using angles. To model a physical system, like a binary star system, you need the actual physical space between the components. Many people try to compute a three-dimensional separation using just right ascension and declination. That does not work. To use the three-dimensional separation method, you must provide a distance parameter when creating the coordinate object. You pass a distance value with a physical unit, like parsecs, alongside your angular coordinates. Once both objects contain this distance data, you call the three-dimensional separation method on the first coordinate, passing the second coordinate as an argument. Astropy returns the direct physical line between them in three-dimensional space, automatically handling the background trigonometry. Now, those objects are not static. To model their movement over time, you attach velocity data directly to the coordinate object upon creation. You do this by passing proper motion parameters. Specifically, you provide proper motion in right ascension and proper motion in declination. These take angular units per time, like milliarcseconds per year. You also provide a radial velocity parameter, which takes a physical speed, like kilometers per second. When you bundle these together, the coordinate object becomes a full phase-space vector. It holds the object's current location and the exact vector of where it is going. Here is the key insight. Your telescope is also moving through space. Earth rotates on its axis and orbits the sun. If you measure the radial velocity of a star from the ground, your raw measurement is contaminated by Earth's own speed at that exact moment. To share your data with other astronomers, you must strip away Earth's motion. You shift your reference point to a stable frame, usually the center of mass of the solar system. This is called the barycentric radial velocity correction. To calculate this correction in Astropy, you need to define where the observer is using an Earth location object. You create this by passing the exact longitude, latitude, and elevation of your instrument. Alternatively, you can pull an established site directly from Astropy's built-in site registry by querying a known name like the Keck Observatory. Next, you need the exact time the observation took place. With your location and time ready, you call the radial velocity correction method directly on your target's coordinate object. You pass it the observation time and your Earth location. Astropy computes the velocity vector of the Keck Observatory at that specific millisecond, relative to the target, and returns the velocity offset. You add this correction value to your raw measured radial velocity. The result is a clean, standardized measurement referenced to the solar system barycenter. The single most useful takeaway here is that a coordinate is not just a fixed point on a map; it is a full kinematic model that holds the target's position, its trajectory, and the geometry needed to correct for your own motion. Thanks for listening. Take care, everyone.

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 5 of 15. Pandas is amazing for flat data, but what happens when a single column in your dataset needs to hold a time series, a physical unit, and a spatial coordinate? Standard dataframes break down. To solve this, you need Tabular Data Mastery: The QTable Class. If you work in Python, you already know Pandas. Pandas assumes your data is flat, preferring one simple value per cell. In astronomy, data is rarely simple. A single observation might require strict units, a matrix of values, or a complex celestial coordinate. Astropy provides a custom tabular data structure specifically designed for this reality. You will notice Astropy has both a standard Table class and a QTable class. People often confuse the two. The standard Table tracks unit information as background metadata, but it returns plain NumPy arrays when you extract a column. A QTable enforces strict unit-awareness. When you retrieve a column from a QTable, it guarantees you get a native Astropy Quantity object back. The physical unit is permanently attached to the numbers. Let us walk through building a QTable. You start by importing QTable from the astropy dot table module. You can construct it using a dictionary where keys are your column names and values are your data arrays. Suppose you are creating a simple catalog of galaxies. Your first column is plain text. You assign a list of strings representing the galaxy names. Now, you need to record where these galaxies are. You do not want separate, disconnected columns for right ascension and declination. Instead, you create an Astropy SkyCoord object containing the positions for all your galaxies. You pass this single SkyCoord object directly into your table as the second column. This is a mixin column. Mixins allow complex Astropy types, like SkyCoord or Time objects, to act exactly like standard data arrays inside the table. You can slice, sort, and filter the table rows, and the mixin object automatically keeps its internal data perfectly synchronized. Next, you need to store flux measurements across three different wavelength bands for each galaxy. In a typical dataframe, you would be forced to create three separate columns. In a QTable, you use a multidimensional column. You define a single array with a shape of N rows and three columns, attach a physical flux unit to the whole thing, and assign it as your third table column. This is where it gets interesting. The QTable accepts the multidimensional array natively. A single row in that column now holds a fully contained array of three unit-aware measurements. When you print this table, Astropy formats it cleanly. It aligns the nested arrays and displays the units directly under the column headers. Your metadata and your measurements are locked together. The QTable is not just a spreadsheet; it is a specialized container built to respect the physics and geometry of your data. Thanks for hanging out. Hope you picked up something new.

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 6 of 15. Searching a massive catalog for a specific source usually means executing a full linear scan that eats up processing time. But if you structure your data properly, you can reduce those lookups from minutes to milliseconds. That is what we are covering today: Advanced Table Operations: Masking and Joins. Astropy tables are not just arrays. They act as an in-memory relational database that natively understands physical units. Suppose you have a photometric catalog and a spectroscopic catalog, and you need to combine them based on a shared source ID. Before you merge them, you have to handle missing data. Your photometric catalog likely has missing magnitude values because an observation failed or a sensor dropped a frame. You might try to filter out those missing values using standard numpy boolean arrays. Do not do that. Standard numpy masks applied directly to a unit-aware table will silently strip away your physical units and break mixin objects. Instead, Astropy handles this with the MaskedColumn. When you work inside a QTable, which is the unit-aware table class, Astropy automatically uses MaskedColumn operations under the hood for missing data. This ensures your magnitude column keeps its physical units intact while safely marking the missing entries so they do not corrupt your downstream calculations. With your cleanly masked photometric catalog ready, it is time to combine it with your spectroscopic data. Astropy provides a join function that works much like a relational database. You pass the function your left table, your right table, and specify the shared column name as the key, which in this case is your source ID. If you specify an inner join, the function evaluates both datasets and returns a new table containing only the rows where the source ID exists in both catalogs. Here is the key insight. Because you are using QTable objects, the join operation automatically reconciles and preserves the units from both sides. Your merged table will contain the precise physical parameters from the spectroscopic data alongside the properly masked magnitudes from the photometric data, perfectly aligned. Now you have a merged catalog, and your analysis requires pulling the exact parameters for one specific source ID. Without intervention, finding a row in a massive table requires checking every single entry sequentially. Astropy solves this performance bottleneck with table indexing. You can instruct the table to create an index on your source ID column. Behind the scenes, Astropy builds a B-tree data structure mapping your source IDs directly to their underlying row locations in memory. Once the index is built, you retrieve your data using the table location property, passing in your target ID. The engine traverses the B-tree, skips the full scan entirely, and fetches your row in logarithmic time. You can even index multiple columns simultaneously to handle complex lookups, such as finding a row by both source ID and observation date. Building an index takes a fraction of compute time upfront, but it pays off immediately when you need to run repeated lookups across heavy datasets. The true power of Astropy tables is not just storing data, but maintaining the physical integrity of your units across complex relational joins and high-speed B-tree queries. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 7 of 15. You have an old, messy text file of star coordinates, and you need to convert it into a modern, compressed binary standard. In most libraries, that means writing a custom parser, mapping data types, and fighting with binary serializers. In Astropy, you can accomplish this with exactly two lines of code. This is thanks to the Unified I/O Interface. The Unified I/O Interface is a polymorphic data handling layer. Instead of forcing you to import a specific module for every file type, Astropy exposes a generic read and write method directly on its core data structures. Whether you are dealing with a standard Table, a QTable with physical units, or a multidimensional NDData array, the interaction model remains identical across entirely different file standards. There is a common misconception about extracting data. When developers need to read a FITS file, they often go straight for the low-level module, astropy dot io dot fits. That module is strictly for raw file manipulation, like editing individual header cards or inspecting unformatted byte streams. If you just want to extract tabular data from a FITS file, you should not use it. You use the Table dot read method instead. It bypasses the manual header extraction and handles the underlying data mapping automatically. When you call the read method and pass a file path or an open file object, Astropy relies on a robust format inference engine. First, it checks the file extension. If you provide a string ending in dot h5 or dot parquet, it routes the request to the corresponding specialized parser. If the extension is ambiguous, missing, or just a generic text extension, the engine drops down a level and inspects the file contents. It reads the first few bytes of the data looking for known signatures or magic numbers to reliably identify the standard. If the automatic inference fails completely, you bypass it by passing a format argument string, like ascii dot csv, to force the correct parser. Here is how this handles our initial scenario. You have that old ASCII text table of star coordinates. You call QTable dot read and pass it your text file path. Astropy detects the ASCII format, parses the columns, and returns a fully populated QTable object in memory. Next, you want to save this as a compressed FITS binary table. You take that same QTable object in memory and call its write method. You provide a new filename ending in dot fits, and you pass an extra parameter to enable compression. You must also pass an overwrite flag set to true, because the write method safely protects existing files by default. This is where it gets interesting. Astropy dynamically switches context from plain text parsing to complex binary serialization under the hood. The underlying file protocols are completely distinct, but the interface you use remains unchanged. Because the interface is attached to the data structure rather than the file format, your data processing pipelines become decoupled from your storage mediums. The single most useful takeaway here is that changing how your application stores and shares data requires altering a file extension in a single string, rather than rewriting your entire data ingestion pipeline. That is it for today. Thanks for listening — go build something cool.

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 8 of 15. The FITS file standard was created in the 1980s for Fortran systems, storing metadata in rigid eighty-byte blocks modeled directly after physical punched cards. Dealing with punched card structures in modern code is rarely pleasant. Astropy hides this complexity by translating those blocks into native structures, which is exactly why this episode focuses on Demystifying FITS Headers and HDUs. Think of a FITS file as a container. Inside this container are Header Data Units, universally known as HDUs. Every FITS file must contain at least one of these, called the Primary HDU. Often, this primary unit contains nothing but metadata, while the actual images or binary tables follow directly behind it as Extension HDUs. When you open a FITS file using Astropy, it parses this container and gives you back an HDUList object. You navigate this object using standard bracket notation to select individual HDUs. Here is the key insight about indexing. The official FITS standard strictly uses one-based indexing because of its Fortran roots. However, since Astropy is a Python library, it enforces zero-based indexing. The Primary HDU is always at index zero, and the first extension is at index one. It is a simple detail, but forgetting it causes constant off-by-one errors. Every HDU object contains a header attribute. This is where those old eighty-byte blocks live. In FITS terminology, each line of metadata is called a card. A complete card holds a keyword, a value, and sometimes a descriptive comment. Astropy processes these cards and exposes the header to you as an object that acts exactly like a standard Python dictionary. If you need the exposure time, you simply pass the EXPTIME keyword to the header. Modifying values works the same way. You assign a new value to the keyword using an equals sign. If you also need to update the comment on that specific card, you assign a tuple containing both the new value and the new comment string. Astropy automatically repacks this data back into the required eighty-byte standard behind the scenes. This translation works flawlessly until you encounter legacy data. Older FITS files frequently violate the strict standard. They might feature keywords that are too long or contain illegal characters. Because Astropy strictly validates headers by default, attempting to open a non-compliant file will often raise an exception immediately. Instead of abandoning the file, you can manage this using the verify method. Say you open an old telescope image and hit a validation error regarding the formatting. You capture the file object and call its verify method, passing the string argument fix. Astropy will scan all the header cards, automatically repair formatting issues where possible, and suppress errors for the corrupted fields it fixes. Once the header is stabilized, you can safely update a broken OBSERVER card with the correct string, and then write your cleaned-up HDUList back to disk. Always treat your FITS metadata updates as simple dictionary operations, but keep that verify method ready for the moment legacy data breaks the rules. If you want to help keep the show going, search for DevStoriesEU on Patreon. Thanks for tuning in. Until next time!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 9 of 15. You do not need a supercomputer to analyze a one hundred gigabyte image. You just need a way to trick your operating system into reading only the exact pixels you care about. Handling Massive FITS Files and Cloud Storage is the exact mechanism you use to accomplish this. When working with massive astronomical arrays, trying to load an entire file into memory will immediately crash your process. Astropy avoids this using a feature called memory mapping. When you open a FITS file, you can pass an argument called memmap set to True. A memory map does not read the image data into your physical RAM. Instead, it maps the file on your hard drive directly into your computer's virtual memory address space. When you access a specific slice of the array, your operating system fetches only that specific block of data from the disk. You get the convenience of pretending the entire file is loaded, but you only pay the physical memory cost for the pixels you actively process. This behavior introduces a very common trap. You might open a FITS file, assign the image array to a variable, and then call close on the file object. You expect the file handle to be released and the resource freed. It is not. Because memory maps rely on the operating system linking the file to a variable, the file remains open and locked as long as any reference to that array exists in your Python session. To actually free the resource and close the file handle completely, you must explicitly delete your data variable using the standard Python delete command. That handles local files, but modern astronomy datasets often live on remote servers. Downloading hundreds of gigabytes just to inspect a small region is a waste of bandwidth and time. Astropy handles remote data by integrating with the file system spec library. When you open a file, you pass an argument called use fsspec set to True. This tells the library to treat cloud storage locations just like local drives. Consider a two hundred megabyte Hubble Space Telescope image hosted in a public Amazon S3 bucket. You only want a tiny ten by twenty pixel cutout around a specific star. First, you call fits open, passing the S3 uniform resource identifier instead of a local file path, alongside use fsspec set to True. Astropy connects to the bucket and downloads only the FITS header blocks. This lets it understand the structure of the file without touching the actual image data. Next, you navigate to the correct image extension and slice the data array using standard indexing, grabbing just your ten by twenty block. This is where it gets interesting. Astropy translates that array slice into specific HTTP byte range requests. Your machine reaches out to Amazon S3, asks for only the physical bytes corresponding to those two hundred pixels, and downloads exactly that tiny chunk. You get your cutout instantly, while the remaining one hundred and ninety nine megabytes of the file stay on the server untouched. Combining local memory mapping with remote byte range streaming means your analysis scripts are constrained only by the data you actively compute, never by the total size of the file itself. I would like to take a moment to thank you for listening — it helps us a lot. Have a great one!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 10 of 15. If you slice a standard numerical array, you only get raw pixels out. But telescope data is never just a matrix of numbers. When you crop an image to look at a single star, your error bounds, your dead pixel masks, and your sky coordinates must stay perfectly aligned, otherwise your entire analysis breaks. That is exactly why we use Gridded Data: The NDData and CCDData Classes. At its core, Astropy provides the NDData base class as a standard container for gridded astronomical data. While NDData handles the general structure, for actual optical or infrared images, you will typically use its specialized subclass called CCDData. This container exists to tie your raw pixel values directly to their physical properties. Let us walk through initializing a CCDData object for a raw telescope exposure. You do not simply pass in a two-dimensional grid of numbers. You pass the raw data grid, and you must assign a physical unit, such as electrons or analog-to-digital units. Next, you attach a boolean mask. If a high energy cosmic ray strikes your detector during the exposure, you flag those specific ruined pixels as true in a separate mask array, and you hand that directly to the CCDData object. Then comes the error tracking. You attach an uncertainty array. For our raw image, you might calculate the Poisson noise and attach it using either the Standard Deviation Uncertainty class or the Variance Uncertainty class. Now, your brightness values, your units, your bad pixel flags, and your statistical noise limits are locked together as one unified object. This brings up a common confusion. Once you have this fully loaded object, do not try to use standard mathematical operators, like a simple plus or minus sign, to process the image. If you attempt normal array arithmetic to subtract a dark background frame, you will leave your uncertainties behind. Instead, you must use the built-in arithmetic methods provided by the class, such as the add, subtract, multiply, and divide functions. This is where it gets interesting. When you call the subtract method, Astropy does not just subtract the pixel values. It automatically propagates your variance uncertainties through the operation using standard statistical rules, yielding a mathematically sound noise profile in the resulting image. Later in your workflow, you might realize your full telescope image is too large, and you only want to analyze a single galaxy located in the top right corner. You accomplish this using a tool called Cutout2D. You provide your main CCDData object, the exact center coordinates of your target galaxy, and the bounding box size you want to extract. Because you packaged everything in a CCDData container, Cutout2D handles the complex syncing automatically. It returns a new, smaller object. It slices the image pixels, but it also automatically crops the cosmic ray mask and the variance uncertainty array to the exact same dimensions. Crucially, it translates the underlying world coordinate system. The pixel coordinate zero-zero in your new cutout correctly maps to the exact same right ascension and declination as it did in the master image. By treating the pixel data, the physical error margins, and the spatial coordinates as one indivisible unit, these classes prevent silent data corruption as your image moves through complex analysis pipelines. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 11 of 15. Every pixel on a telescope detector is a flat, discrete square, but the sky is a continuous, mathematically complex curved sphere. Mapping a point on that flat detector back to its true position in the universe is a messy geometric problem. World Coordinate Systems, or WCS, is the mechanism that decodes that geometry. WCS acts as the ultimate bridge between hardware and physics. It provides the mathematical transformation required to turn a pixel coordinate into a world coordinate, like Right Ascension and Declination. In astronomical data, these transformation rules are usually stored as standardized metadata inside a FITS file header. The header contains specific keywords defining the reference pixel, the physical coordinates of that pixel, and the rotation or scaling matrix used by the telescope. Astropy reads these keywords and builds an active transformation object that models the entire focal plane. Before doing any math, we must clear up a constant source of indexing confusion. The FITS standard defines pixel coordinates as one-indexed, meaning the very first pixel on the detector is pixel one. Python, and therefore NumPy, is zero-indexed. Astropy resolves this duality through two different layers. The older, low-level API is not a unified indexing system; it requires you to explicitly pass an origin argument of zero or one every single time you transform a coordinate. But the modern, high-level Shared Python Interface operates entirely on Python's zero-indexed convention. It expects zero-based pixels from you and automatically handles the offset to the FITS standard internally. Here is the key insight. The high-level API makes coordinate conversion incredibly concise. Say you have located a bright star in your image array at pixel x equals 30 and y equals 40, and you need to find its true position in space. First, you read your FITS file header using Astropy's input-output tools. Next, you initialize a WCS object by passing that header directly into the WCS class. Finally, you call a method named pixel to world on that WCS object, passing in your x and y pixel values of 30 and 40. The method evaluates the projection math and returns a standard Astropy SkyCoord object. This object contains the physical coordinates and is completely aware of its reference frame, such as ICRS. The beauty of this shared interface is that it shields you from raw numbers. Instead of returning a generic array of floats representing abstract degrees, it gives you rich objects that understand their own physical units and coordinate systems. You can take that SkyCoord and immediately use it to cross-reference a star catalog. The system works exactly the same in reverse. If you have the coordinates of a known galaxy and want to determine exactly which pixels it occupies on your camera sensor, you call the method world to pixel on your WCS object. You pass in your SkyCoord, and it returns the exact x and y floating-point pixel positions. The high-level WCS API isolates your application logic from the underlying projection math, meaning your pixel-to-sky conversion code remains identical whether you are working with a flat radio map or a deeply distorted wide-field optical mosaic. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 12 of 15. In many scientific libraries, the mathematical equation you want to evaluate is tightly coupled to the algorithm used to optimize it. If you want to change how the optimization works, you often have to rewrite the way you defined your math. Astropy eliminates this problem entirely with its approach to Analytical Models and Fitting. The core design philosophy of the astropy modeling subpackage is a strict separation of concerns. The mathematical definition of your equation is one object. The algorithmic engine that fits that equation to your data is a completely separate object. This means you can define your math once, and then test it against several different fitting algorithms simply by passing the model into a different fitter. Astropy provides dozens of predefined mathematical models. If you are modeling a steady trend, you use a Linear1D model. If you are studying an emission line in a spectrum, you reach for a Gaussian1D model. These model objects encapsulate their parameters, such as the amplitude, mean, and standard deviation. Before we look at the fitting process, we need to address a frequent source of confusion. In Astropy, models are callable objects. You pass an array of coordinates directly to an instantiated model, and it evaluates the math to return the corresponding values. Most importantly, when you optimize a model, the fitter does not mutate your original starting object. It always returns a brand new fitted model instance. Your original guess is preserved exactly as you defined it. Consider a scenario where you have noisy one-dimensional array data representing a spectral feature. First, you initialize a Gaussian1D model with your rough initial guesses for the amplitude, mean, and width. Because Astropy separates the math from the optimizer, you attach your physical constraints directly to the model parameters. Suppose you already know the exact center coordinate of this spectral line. You access the mean parameter on your model and set its fixed property to true. The optimization algorithm will now leave that value alone. You can also establish mathematical boundaries. By setting a minimum bound of zero on the amplitude parameter, you force the algorithm to reject any solutions that result in a negative peak. With your starting model fully constrained, you bring in the fitter. For non-linear equations like a Gaussian, you instantiate the Levenberg-Marquardt least squares fitter, known in Astropy as LevMarLSQFitter. You execute the fit by calling this fitter object, passing in your initial Gaussian model along with your horizontal and vertical data arrays. The algorithm runs and outputs a completely new Gaussian1D model containing the optimized parameter values. Because this new model is directly callable, you just pass your original horizontal array into it to generate a smooth mathematical curve that you can plot against your noisy data. Here is the key insight. By storing constraints on the model rather than passing them to the solver, the fitting algorithms remain generic. You spend your effort encoding the physical reality of your problem into the mathematical model, allowing Astropy to swap out the optimization logic in the background without breaking your pipeline. That is all for this one. Thanks for listening, and keep building!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 13 of 15. Instead of writing a massive fifteen parameter mathematical function by hand to fit a messy spectrum, what if you could just add basic shapes together like Lego bricks and let Python handle the calculus? That is exactly what we are covering today with Compound Models and Custom Fits. In Astropy, you rarely need to build complex analytical models from scratch. You take simple, predefined blocks and combine them using standard Python arithmetic operators like addition, subtraction, multiplication, or division. When you add two Astropy models together, you do not get a static array of evaluated numbers. You get a brand new, fully functioning model object. This compound structure knows how to evaluate itself, calculate its derivatives, and track every parameter from the original base models. Consider a specific scenario. You have a spectrum containing two distinct emission lines sitting on top of a flat background continuum. You want to fit this entire structure simultaneously. You start by instantiating a basic polynomial model to act as your flat baseline. Next, you instantiate two separate one dimensional Gaussian models to represent the two emission lines. To create your final spectral profile, you simply define a new variable equal to the polynomial plus the first Gaussian plus the second Gaussian. Astropy automatically merges them into a single compound evaluation tree. This introduces a common trap regarding parameter access. When you create this compound model, you might wonder how to isolate the amplitude of the second emission line. Since both are Gaussians, they both have a parameter named amplitude. Astropy resolves this collision by automatically appending numeric suffixes based on the order the models were combined. The first Gaussian gets its parameters renamed to amplitude zero, mean zero, and standard deviation zero. The second Gaussian gets amplitude one, mean one, and so on. If you ever lose track of what parameter belongs to which component, you can inspect the parameter names attribute of your compound model to see the exact generated list. Now, the second piece of this is fitting data that carries physical units. Astronomical measurements are not just raw floats. Your independent variable might be wavelength in Angstroms, and your dependent variable might be flux density in Janskys. Astropy models handle this elegantly through the Quantity object. When you pass your data arrays into a fitter, you pass the Astropy Quantity objects directly, without stripping the units away. The fitter analyzes the units of your input data and the structure of your compound model. It catches dimensional mismatches before the fit even begins. Once the fit converges, the parameters on the resulting model will automatically hold the correct physical units. Your Gaussian means update to Angstroms, and your amplitudes update to Janskys. You do not have to write custom conversion logic or strip units before fitting only to reattach them later. The compound model enforces dimensional consistency throughout the entire process. Here is the key insight. When you use a plus sign between two models, you are not just chaining functions together. You are constructing a single unified equation that maintains strict unit awareness and parameter traceability from the raw data all the way to the final fit. That is your lot for this one. Catch you next time!

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 14 of 15. You have fourteen thousand unevenly spaced telescope photos of a star. Somewhere in that noise is a zero point zero one percent dip in light caused by an alien planet passing in front of it. Here is the exact algorithm to find it. This episode is about Time Series Analysis and hunting exoplanets. Astropy handles this workflow using the TimeSeries class. Before going further, let us clear up a common misconception. A TimeSeries is not a completely new data structure with its own rules. It is simply a subclass of QTable. This means every table manipulation method, column operation, and masking technique you already know still works perfectly. The only difference is that the first column is strictly enforced as an Astropy Time object. This strictness guarantees your timestamps and your flux measurements stay perfectly aligned, whether your observations are evenly spaced or entirely random. To start hunting an exoplanet, you first load a light curve. Astropy provides built-in readers, allowing you to pass a Kepler FITS file directly to the TimeSeries read function. It automatically sets up the Time column and pulls in the light measurements, usually called flux. Kepler data is extremely dense, sometimes containing thousands of data points over a few months. This raw data often contains high-frequency noise from the spacecraft or stellar activity. You can smooth it out using downsampling via aggregation. You define a time bin size, like ten minutes, and pass an aggregation function such as the mean. Astropy groups the timestamps into these bins and averages the flux, reducing the size of the dataset and suppressing random noise without erasing the actual transit signal. Now you have clean data, but the exoplanet transits are hidden somewhere along a massive timeline. To find the planet, you need a periodogram. You might be familiar with the Lomb-Scargle periodogram, which is excellent for finding smooth, continuous sine waves like the light from a pulsating star. But a planet passing in front of a star does not create a sine wave. It creates a flat line, a sudden drop, a flat bottom, and a sudden rise. The signal is a box. Because of this shape, you use the Box Least Squares periodogram, or BLS. You pass your time column and your flux column into the BLS function. The algorithm then evaluates an entire grid of possible orbital periods, perhaps testing everything from a one-day orbit to a fifty-day orbit. At each tested frequency, it slides a mathematical box across the data, trying to match the depth and duration of a potential planetary transit. It returns a power spectrum showing how well the box fit at every single frequency. You simply extract the period that generated the highest power peak. That is your suspected exoplanet orbit. Once you have that exact period, you need to verify it visually. You take your original TimeSeries and call the fold method, passing in the period you just discovered. Folding takes every single observation across months or years and maps it to a single orbital cycle. Instead of a chronological timeline, your horizontal axis becomes the phase, ranging from zero to one. Every time the planet transited the star over the course of the mission, those individual data points are stacked on top of each other at the exact same phase. This is the part that matters. When you plot that folded TimeSeries, the random timeline scatter vanishes, and a clean, undeniable U-shaped dip appears precisely in the center of your graph. You have just used geometry and time to pull an invisible world out of the noise. Thanks for spending a few minutes with me. Until next time, take it easy.

Hi, this is Alex from DEV STORIES DOT EU. Astropy: Python for Astronomy, episode 15 of 15. Solving the metric of the expanding universe usually requires evaluating massive elliptic integrals by hand. Or, you can just import a single module and let Python do the math. Today we are looking at Cosmological Calculations: Sizing up the Universe. The astropy cosmology subpackage encapsulates general relativity and the Friedmann-Lemaitre-Robertson-Walker metric into effortless function calls. Instead of writing your own integrators to map redshift to physical distance, you define a cosmological model, and the library handles the underlying geometry. A cosmology in Astropy is a Python object. You can build one from scratch using the FlatLambdaCDM class. You pass it a Hubble constant, like seventy, and an initial matter density, like point three. This gives you a spatially flat universe driven by dark energy and cold dark matter. However, you rarely need to type these parameters yourself. Astropy includes standard built-in realizations derived from major surveys. You simply import Planck13 or WMAP9 from the cosmology module, and you instantly have a fully configured model ready to evaluate. Here is the key insight regarding state. The parameters of an Astropy cosmology object are strictly immutable. If you initialize a model and then decide you want to tweak the matter density, you cannot update that attribute in place. Doing so throws an error. Instead, you call the clone method on your existing model, passing the new parameter value as an argument. This returns a brand new cosmology instance with your updated values. Take a concrete scenario using the Planck13 built-in model. Suppose you are observing a galaxy at a redshift of two. You want to know its lookback time, meaning how many years the light took to reach your telescope. You take your Planck13 object, call its lookback time method, and pass it the number two. The method returns a time quantity in gigayears. You use this exact same pattern to calculate luminosity distance. You call the luminosity distance method with your redshift, and Astropy returns the distance in megaparsecs, factoring in the expansion of the universe. That covers inputs to outputs, but what about the inverse? Sometimes you have a physical measurement and you need to find the corresponding redshift. Say you want to know the exact redshift at which the universe was exactly two billion years old. You cannot simply pass an age into the standard methods. Instead, you import a function called z at value. You pass it the method you want to invert, which is the age method of your Planck13 model, and you pass it the target value of two gigayears. Astropy runs a numerical root-finding routine behind the scenes and returns the precise redshift. Converting between redshift and physical time or distance is the backbone of extragalactic astronomy, and this module gives you a tested source of truth so you can focus on data instead of debugging integrals. Since this is the final episode of our Astropy series, I encourage you to read through the official documentation, try these tools hands-on, or visit dev stories dot e u to suggest topics for our next series. Thanks for spending a few minutes with me. Until next time, take it easy.