Search code examples
pythonpython-xarrayrasterioproj

Regrid geostationary GOES satellite data to a cartesian grid


I am working with GOES satellite data, which is gridded in radians and has the standard projection of geostationary. I am aware of several methods that can reproject the data such that it can be plotted in a cartesian sense. However, I have not been able to find any resources that would allow me to actually regrid the data such that longitude and latitude would actually be coordinates of the dataset. I want to have the data in this format so that I can easily index the data. Any ideas on what to do here?

Sample data can be found here: https://noaa-goes16.s3.amazonaws.com/ABI-L1b-RadC/2020/230/00/OR_ABI-L1b-RadC-M6C07_G16_s20202300001171_e20202300003555_c20202300004038.nc

I am loading and working with the data using xarray.


Solution

  • Disclaim: Super biased answer below. I am one of the core developers on the Satpy and Pyresample packages.

    The Satpy library would be a good place to start as it should make it simple to workaround the stuff you might not care about (conversion from radiance to brightness temperature, conversion from projection space to lon/lat space, etc). See this example jupyter notebook: https://github.com/pytroll/pytroll-examples/blob/master/satpy/GOES-16%20Mosaic.ipynb

    You can also find a tutorial I taught at a SciPy conference a while ago that includes a section on resampling: https://github.com/pytroll/tutorial-satpy-half-day/blob/master/notebooks/04_resampling.ipynb

    With this information you should be able to resample ABI data to your own custom grid (AreaDefinition in pyresample terms). There is more information on making your own pyresample AreaDefinition here: https://pyresample.readthedocs.io/en/latest/geometry_utils.html

    Here is a basic outline of what your code might look like:

    from satpy import Scene
    from pyresample import create_area_def
    
    # analyze filenames and what is in them
    scn = Scene(reader='abi_l1b', filenames=abi_filenames)
    
    # load a specific band
    scn.load(['C07'])
    
    # create a custom area that fits the projection and resolution you want
    your_area_def = create_area_def(...)
    
    # resample all the loaded datasets to your area definition
    new_scn = scn.resample(your_area_def)
    
    # save data as geotiffs to the current directory
    new_scn.save_datasets()
    

    Otherwise, if you are really interested in using xarray and rasterio for this you may want to look at the rioxarray project: https://github.com/corteva/rioxarray