diff --git a/pyorbital/geoloc_example.py b/pyorbital/geoloc_example.py index 29d6fcf..c3c1296 100644 --- a/pyorbital/geoloc_example.py +++ b/pyorbital/geoloc_example.py @@ -29,34 +29,42 @@ from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt +# Couple of example Two Line Elements tle1 = "1 33591U 09005A 12345.45213434 .00000391 00000-0 24004-3 0 6113" tle2 = "2 33591 098.8821 283.2036 0013384 242.4835 117.4960 14.11432063197875" +# Choosing a specific time, this should be relatively close to the issue date of the TLE t = datetime(2012, 12, 12, 4, 16, 1, 575000) - -scanline_nb = 351 - +# this is the number of full scan rotations +scans_nb = 10 # we take only every 40th point for plotting clarity scan_points = np.arange(24, 2048, 40) - +# This the maximum scan angle away from nadir for the given TLE that still sees earth. +scan_angle = 55.37 +# period of one full rotation 1/6 s +scan_p = 0.16666667 +# integration time of instrument +int_t = 0.000025 # build the avhrr instrument (scan angles) -avhrr = np.vstack(((scan_points - 1023.5) / 1024 * np.deg2rad(-55.37), - np.zeros((len(scan_points),)))).transpose() -avhrr = np.tile(avhrr, [scanline_nb, 1]) +# creates list of radian angles centered around nadir based on the scan points that see earth +avhrr = np.vstack(((scan_points / 1023.5-1) * np.deg2rad(-scan_angle), + np.zeros((len(scan_points),)))) +avhrr = np.tile( + avhrr[:, np.newaxis, :], [1, scans_nb, 1]) # building the corresponding times array -offset = np.arange(scanline_nb) * 0.1666667 -times = (np.tile(scan_points * 0.000025 + 0.0025415, [scanline_nb, 1]) - + np.expand_dims(offset, 1)) +times = np.tile(scan_points * int_t, [scans_nb, 1]) +offset = np.arange(scans_nb) * scan_p +times += np.expand_dims(offset, 1) # build the scan geometry object -sgeom = ScanGeometry(avhrr, times.ravel()) +sgeom = ScanGeometry(avhrr, times) -# roll, pitch, yaw in radians +# roll, pitch, yaw in radians. This is a static offset. rpy = (0, 0, 0) -# print the lonlats for the pixel positions +# print the longitude and latitude for the pixel positions s_times = sgeom.times(t) pixels_pos = compute_pixels((tle1, tle2), sgeom, s_times, rpy) pos_time = get_lonlatalt(pixels_pos, s_times)