Latitude and Longitude for GK2A AMI products

[Frankie91]

Hello everyone,

Hope this is the right place to ask this question. I am currently trying to obtain estimates of the latitude and longitude for each pixel of GK2A AMI AOD L2 products for its 'East Asia' observation area. I first wanted to test whether I could, with some adaptation, use the L1 reader provided by this package, but I don't see it when running available_readers().

My second attempt was to adapt the following calc_latlon function (original source: https://github.com/lsterzinger/goes-l2-latlon-tutorial/blob/main/tutorial.ipynb) which I already use to make identical estimates for identical GOES ABI L2 products:

def calc_latlon(ds):
x = ds.x
y = ds.y
goes_imager_projection = ds.goes_imager_projection

x,y = np.meshgrid(x,y)

r_eq = goes_imager_projection.attrs["semi_major_axis"]
r_pol = goes_imager_projection.attrs["semi_minor_axis"]
l_0 = goes_imager_projection.attrs["longitude_of_projection_origin"] * (np.pi/180)
h_sat = goes_imager_projection.attrs["perspective_point_height"]
H = r_eq + h_sat

a = np.sin(x)**2 + (np.cos(x)**2 * (np.cos(y)**2 + (r_eq**2 / r_pol**2) * np.sin(y)**2))
b = -2 * H * np.cos(x) * np.cos(y)
c = H**2 - r_eq**2

r_s = (-b - np.sqrt(b**2 - 4*a*c))/(2*a)

s_x = r_s * np.cos(x) * np.cos(y)
s_y = -r_s * np.sin(x)
s_z = r_s * np.cos(x) * np.sin(y)

lat = np.arctan((r_eq**2 / r_pol**2) * (s_z / np.sqrt((H-s_x)**2 +s_y**2))) * (180/np.pi)
lon = ((l_0 - np.arctan(s_y / (H-s_x))) * (180/np.pi))
ds = ds.assign_coords({
    "Latitude":(["y","x"],lat),
    "Longitude":(["y","x"],lon)
})
ds.Latitude.attrs["units"] = "degrees_north"
ds.Longitude.attrs["units"] = "degrees_west"
return ds

The Netcdf files provided by the South Korean Meteorological Agency, however, do not include the required information about the horizontal and vertical scan angles (in radians), used by the function to define x and y at the beginning. Everything else is provided, as the metadata for an example file show:

<xarray.DataArray 'gk2a_imager_projection' ()>
[1 values with dtype=int32]
Attributes: (12/13)
    long_name:                       AMI fixed grid projection
    grid_mapping_name:               geostationary
    latitude_of_projection_origin:   0.0
    longitude_of_projection_origin:  128.1999969482422
    perspective_point_height:        35785864.0
    semi_major_axis:                 6378137.0
                             ...
    column_scale_factor:             20425338.0
    line_scale_factor:               20425338.0
    column_offset:                   1125.5
    line_offset:                     2337.5
    sweep_angle_axis:                x
    spatial_resolution:              2km at nadir

Does anyone know if there is an alternative approach, either in satpy or with other packages, to estimate these scan angles for each GK2A AMI pixel? Or, alternatively, of another way to obtain the same latitude and longitude estimates for GK2A AMI products?

Thank you in advance for your attention!

[mraspaud]

Hi @Frankie91,

I'm not familiar with the AMI data, but satpy is capable of reading the level 1b netcdf files and thus provide the longitudes and latitudes for every pixel.

Something like this should work:

from satpy import Scene
import glob


if __name__ == "__main__":
    curfiles = glob.glob('/path/to/my/ami/files/*.nc')
    scn = Scene(curfiles, reader='ami_l1b')
    channel = 'C05'
    scn.load([channel])
    lons, lats = scn[channel].attrs["area"].get_lonlats()

Tell me if it works!

[Frankie91]

Hi @mraspaud,

Thank you for the prompt answer. I've tried your code snippet with L1b and L2 files both for my region of interest and for the full AMI disk, but I always get the same identical error message:

`INFO:satpy.readers:Cannot use ['C:\\Users\\Dfran\\AppData\\Roaming\\Python\\Python39\\site-packages\\satpy\\etc\\readers\\ami_l1b.yaml']
WARNING:satpy.readers:Don't know how to open the following files: {'C:\\Users\\Dfran\\.spyder-py3\\EastAsia\\gk2a_ami_le1b_nr013_ea020lc_202302270130.nc', 'C:\\Users\\Dfran\\.spyder-py3\\EastAsia\\gk2a_ami_le1b_nr013_ea020lc_202302270120.nc', 'C:\\Users\\Dfran\\.spyder-py3\\EastAsia\\gk2a_ami_le1b_nr013_ea020lc_202302270140.nc'}
Traceback (most recent call last):

  Cell In[72], line 3
    scn = Scene(curfiles, reader='ami_l1b')

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\scene.py:133 in __init__
    self._readers = self._create_reader_instances(filenames=filenames,

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\scene.py:154 in _create_reader_instances
    return load_readers(filenames=filenames,

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\readers\__init__.py:579 in load_readers
    _check_reader_instances(reader_instances)

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\readers\__init__.py:618 in _check_reader_instances
    raise ValueError("No supported files found")

ValueError: No supported files found`

[simonrp84]

Do you have netcdf installed? conda install netcdf

[mraspaud]

Maybe this can help: https://satpy.readthedocs.io/en/stable/quickstart.html#troubleshooting

[Frankie91]

@simonrp84 I assume you meant netcdf4 which yes, I have installed(probably got it with xarray).

@mraspaud I do get the following message with regard to the ami_l1b reader: ami_l1b: cannot find module 'satpy.readers.ami_l1b' (No module named 'pyspectral')

[mraspaud]

@Frankie91 just install the pyspectral package then :) should be available on conda forge and pypi

[Frankie91]

@mraspaud That solved it with the reader, however now I get the following error message:

if __name__ == "__main__":
    curfiles = glob.glob('C:\\Users\\Dfran\\.spyder-py3\\FullDisk\\*.nc')
    scn = Scene(curfiles, reader='ami_l1b')
    channel = 'image_pixel_values'
    scn.load([channel])
    lons, lats = scn.attrs["area"].get_lonlats()
Traceback (most recent call last):

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\scene.py:1360 in _update_dependency_tree
    self._dependency_tree.populate_with_keys(needed_datasets, query)

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\dependency_tree.py:262 in populate_with_keys
    raise MissingDependencies(unknown_datasets, "Unknown datasets:")

MissingDependencies: Unknown datasets: {DataQuery(name='image_pixel_values')}


During handling of the above exception, another exception occurred:

Traceback (most recent call last):

  Cell In[98], line 5
    scn.load([channel])

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\scene.py:1348 in load
    self._update_dependency_tree(needed_datasets, query)

  File ~\AppData\Roaming\Python\Python39\site-packages\satpy\scene.py:1362 in _update_dependency_tree
    raise KeyError(str(err))

KeyError: "Unknown datasets: {DataQuery(name='image_pixel_values')}"

For clarity, I've modified the channel name to 'image_pixel_values' because that is the only band available in the netcdf files that can be downloaded from the korean meteorological agency - they do not provide a single netcdf with all bands, but separate ones for each band.

[djhoese]

@Frankie91 I'm not sure I understand. If you have a single L1b file then can you pass that file to the Scene and load the related band for that file? You can also see scn.available_dataset_names() to see what the band name is. You'll want to make sure whatever band file you use that the resolution matches your L2 products so the lon/lats are the same.

[Frankie91]

@djhoese Yes, that was indeed my plan - work with a band with 2 km resolution so I could get lon-lat arrays that I could then also use for the AOD product. Using scn.available_dataset_names() enabled me to see that I was using the wrong dataset name - I've corrected that to 'NR013' and now I can load the scene. Now I'm trying to solve the following error, related to the last line of the snippet provided by @mraspaud:

Traceback (most recent call last):

  Cell In[121], line 6
    lons, lats = scn.attrs["area"].get_lonlats()

KeyError: 'area'

[mraspaud]

@Frankie91 I think you forgot one piece from the line you copied from my code, it should be
lons, lats = scn[channel].attrs["area"].get_lonlats()
:)

[Frankie91]

@mraspaud correct, my bad :) That finally solved it, Thanks everyone!