How compute satellite look angle?

[ahvahsky2008]

Hi. I have in input only tle. How calcullate observation area?

[mraspaud]

Hi,
TLE is a good start, but you need more information than this to compute the observation area, in particular the characteristics of the sensor like the resolution and min/max scanning angle along and across track.

Once you have this, you can try to define a scanning geometry like here: https://github.com/pytroll/pyorbital/blob/main/pyorbital/geoloc_instrument_definitions.py and then use something like this to get the positions of your sensor's pixels on the ground: https://github.com/pytroll/pyorbital/blob/main/pyorbital/geoloc_example.py

Good luck!

[ahvahsky2008]

@mraspaud thx! I try run geoloc_example.py but its raise error
ValueError: operands could not be broadcast together with shapes (3,17901) (1,2)

[mraspaud]

Yes, the script might not be up to date, but if you have the sensor definition set properly, I'm sure you can adapt the example to your use case :)

[ahvahsky2008]

but if you have the sensor definition set properly,

could you pls give some hint? which param
https://i.ibb.co/0hJFHct/xixa.png

[mraspaud]

I have no idea what sensor you are working with, so it's going to be difficult for me to give you a hint 😅

If you paste your sensor definition in the same format as in https://github.com/pytroll/pyorbital/blob/main/pyorbital/geoloc_instrument_definitions.py and a suitable TLE, I'll try to find some time to take a look at this.

[ahvahsky2008]

Thx for reply) Ok. I try work with sat SUOMI NPP All works (may be), but way of satellite incorrect

from geoloc_instrument_definitions import avhrr
import numpy as np
from datetime import datetime
from pyorbital.geoloc import ScanGeometry, compute_pixels, get_lonlatalt
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

from pyorbital.orbital import Orbital

tle1 = "1 37849U 11061A   22207.46972667  .00000025  00000-0  32513-4 0  9995"
tle2 = "2 37849  98.7330 145.4829 0001172  89.7809  20.2392 14.19535383556723"

t = datetime(2022, 7, 26, 19, 4, 1, 575000)

scanline_nb = 10

scan_points = np.arange(24, 2048, 40)

sgeom = avhrr(scanline_nb,scan_points)

rpy = (0, 0, 0)

s_times = sgeom.times(t)
pixels_pos = compute_pixels((tle1, tle2), sgeom, s_times, rpy)
pos_time = get_lonlatalt(pixels_pos, s_times)


m = Basemap(projection='eck4', llcrnrlat=24, urcrnrlat=70, llcrnrlon=-25, urcrnrlon=120,
            lat_ts=58, lat_0=58, lon_0=14, resolution='l')

x, y = m(pos_time[0], pos_time[1])
p1 = m.plot(x, y, marker='+', color='red', markerfacecolor='red', markeredgecolor='red', markersize=1, markevery=1,
            zorder=4, linewidth=0.0)
m.fillcontinents(color='0.85', lake_color=None, zorder=3)
m.drawparallels(np.arange(-90., 90., 5.), labels=[1, 0, 1, 0], fontsize=10, dashes=[1, 0],
                color=[0.8, 0.8, 0.8], zorder=1)
m.drawmeridians(np.arange(-180., 180., 5.), labels=[0, 1, 0, 1], fontsize=10, dashes=[1, 0],
                color=[0.8, 0.8, 0.8], zorder=2)

plt.show(block=True)

[mraspaud]

I'm not sure I understand. The plot shows north America, but the Basemap definition seems to be over Europe?

[ahvahsky2008]

counter question - plot shows position of satellite on proper datetime?

[mraspaud]

yes. Don't forget the time should be given in utc.

Here is the script I used to generate this image:

from pyorbital.geoloc_instrument_definitions import avhrr
import numpy as np
from datetime import datetime
from pyorbital.geoloc import ScanGeometry, compute_pixels, get_lonlatalt
import matplotlib.pyplot as plt

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from pyorbital.orbital import Orbital

def main():

    # Noaa 19, july 27 2022
    tle1 = "1 33591U 09005A   22208.15634740  .00000133  00000+0  97035-4 0  9993"
    tle2 = "2 33591  99.1455 243.1224 0014915  77.9753 282.3088 14.12592326694052"

    t = datetime(2022, 7, 27, 14, 10, 0)

    scanline_nb = 10

    scan_points = np.arange(24, 2048, 40)

    sgeom = avhrr(scanline_nb,scan_points)

    rpy = (0, 0, 0)

    s_times = sgeom.times(t)
    pixels_pos = compute_pixels((tle1, tle2), sgeom, s_times, rpy)
    pos_time = get_lonlatalt(pixels_pos, s_times)

    fig = plt.figure()
    proj4_params = dict(proj='eck4',
                        lat_ts=58, lat_0=58)
    crs = ccrs.PlateCarree(14)
    ax = fig.add_subplot(1, 1, 1, projection=crs)
    ax.set_extent([-25, 120, 24, 90], crs=crs)

    ax.add_feature(cfeature.LAND)
    ax.add_feature(cfeature.OCEAN)
    ax.add_feature(cfeature.COASTLINE)

    x, y = (pos_time[0], pos_time[1])
    p1 = plt.plot(x, y, marker='+', color='red', markerfacecolor='red', markeredgecolor='red', markersize=1, markevery=1,
                zorder=4, linewidth=0.0)

    plt.show()


if __name__ == '__main__':
    main()

[mraspaud]

And what another website reports:

[ahvahsky2008]

Ok, we have satellite trajectory on purpose datetime. Main question - how calc observation area?

[mraspaud]

The script you have provides observation points, right? So from there you can get pixels for given times of interest.

[ahvahsky2008]

I guess this lines are swath ?

[mraspaud]

yes, approximately. You may want to compute some more times in between, a swath is usually more rounded (I'm assuming the satellite goes horizontally in this image)