PyOrbital example

[Rutger]

Hey all,

Its been a while since i used PyOrbital for anything serious but last week i decided to give my earlier attempt of predicting the VIIRS granules another go. I have attached an example result to this post.
The starting time of the granules are determined by peaking at the ladsweb ftp and grab it from the filenames present. The code to generate the granules footprint is still the same as in my post earlier this year. For plotting i only kept the outer edge and converted it to a Matplotlib patch, this allows an intersection test with my AOI (Egypt in this case) which is used for highlighting the granules covering the AOI.

Note how the left and right edges of the granule dont align very well with the next. I mentioned it before but this images shows it more clearly, im still a bit clueless what this is. For what i understand the same code is used the determine the orbit of the satellite, which looks fine. But even with these glitches its already very usable.


Regards,
Rutger

[Martin Raspaud]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hi Rutger,

Funny you bring that up now again, since last week I looked again at
you posts from last time...

I didn't have much time to find out what the problem was, but I guess
we have to look at the geolocation code.

Right now I'm using the satellite track as a rotation vector for the
pixel vectors (vectors pointing at each pixel), and using quaternion
rotation to get the pixels. (We're setting the look angles in the
instrument description).

So two questions:
- - Is that way of doing correct in principle?
- - If it is, is there a mistake in the code?

What do you think?

I'll try to have a look at this later today...

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJQjlTYAAoJEBdvyODiyJI4BcMIAOEOFHcv56pMKOIeOpe3SYRJ
OQfwGOrMWMpAQRKLI1QoyS1/Luq0x0n9wggRaGAa0MoDkdsnMVWJuYAWe8F7ZpE4
SCMI+L9bkNsts1K4UuCejAvhJSJGlEdKrDTJhQNQUAc/eYPNnF4811LAQtPx5K4z
8AAHchacS67ZTrbakOZY5ikycSWp0pJKSRJ8iieyXdmGZMZDtJ4FjRJBlgsk8dDs
pt3T4A5thL3A+gp49cp9YFL9/vqWWegaD2ueW3HBow3twdmm/xSqh+KnVWEAIEo6
Zvk9OPxPrNfysj24WsBFA6XvlOzF1KoGcfVQwBxuF9BUuMGGs0ZoErposliwPU0=
=8ysI
-----END PGP SIGNATURE-----

[Hróbjartur Þorsteinsson]

That issue seems to only affect the calculation of the
swath edges.

I plotted some 5 minute granules with our SPP.py code,
see attachment. I recall I did something rather crude
when calculating the swath edge,

Perhaps digging into "rightSwathEdge" and "leftSwathEdge"
subroutines in SPP can help, I calculated the edges using
great circle arcs based on perpendicularity to the ground
motion of the satellite.

Note the ground motion and satellite-bearing was calculated crudely,
by taking a small step forward, instead of evaluating it from
earth rotation and orbital parameters.

Cheers,
Hrob.

[Martin Raspaud]

On 29/10/12 09:18, Rutger wrote:
> Hey all,

[...]

> Note how the left and right edges of the granule dont align very well with
> the next. I mentioned it before but this images shows it more clearly, im
> still a bit clueless what this is. For what i understand the same code is
> used the determine the orbit of the satellite, which looks fine. But even
> with these glitches its already very usable.

Hi Rutger,

I think I fixed it! The last line of the computation should be:
pos_time = get_lonlatalt(pixels_pos, sgeom.times(t))

instead of:
pos_time = get_lonlatalt(pixels_pos, t)

So instead of using a fixed time (and thus a fixed gmst) for all the
pixels, we use a different time for every pixel...

Best regards,
Martin

[Rutger]

Hey,

You're spot on! This makes my day. :) I spent some time last night reading up on the quaternion rotation, i had never heard of it before. But Hrob's comment on the Earth rotation already gave a hint that it should be some time related issue. The Earth rotates about .3 decimal degrees during the recording of a granule which is also about the deviation at the end of the recording.

Excellent, now the next step on my todo list is validation with the GITCO geolocation product. Which should be pretty straightforward. Every scanline is 32 pixels high so reading out the 32th and -32th pixels at the beginning and the end should be very close to the predicted first and last scanline.

I have attached the improved map.

Regards,
Rutger

[Martin Raspaud]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 30/10/12 08:28, Rutger wrote:
> Hey,

Hi Rutger,

> You're spot on! This makes my day. :) I spent some time last night
> reading up on the quaternion rotation, i had never heard of it
> before. But Hrob's comment on the Earth rotation already gave a
> hint that it should be some time related issue. The Earth rotates
> about .3 decimal degrees during the recording of a granule which is
> also about the deviation at the end of the recording.

> Excellent, now the next step on my todo list is validation with the
> GITCO geolocation product. Which should be pretty straightforward.
> Every scanline is 32 pixels high so reading out the 32th and -32th
> pixels at the beginning and the end should be very close to the
> predicted first and last scanline.

Yes, I checked quickly yesterday, the bow-tie effect is clearly seen
when you take in the first and last pixel of each scan... However, I
had a hard time to find precise values on the internet for the
backward and forward fov's... From what I gathered, the width of the
scan along-track is 11.87km, so if the satellite is about 824km above
the earth, we can approximate the along-track fov phi as:
tan(phi) = 11.87/824
which gives about 14.4 mrad...

Did you find any better value for that?

Here are the changes I made to the code to accommodate for wider scans
(more than one scanline):

#taking just borders and middle for now
scan_pixels = 3
#scan_pixels = 32

across_track = (np.arange(6400) - 3199.5) / 3200 * np.deg2rad(-55.84)

# pitch rotation: np.arctan2(11.87/2, 824.0)

y_max_angle = np.arctan2(11.87/2, 824.0)
along_track = np.array([-y_max_angle, 0, y_max_angle])

scan = np.vstack((np.tile(across_track, scan_pixels),
np.repeat(along_track, 6400))).T

npp = np.tile(scan, [scanline_nb, 1])

[...]

offset = np.arange(scanline_nb) * 1.779166667
times = (np.tile(np.arange(6400) * 0.0002779947917, [scanline_nb,
scan_pixels]) + np.expand_dims(offset, 1))


About the validation against GITCO, I would expect discrepancies due to:
- - the vertical axis points to the center of the earth at the moment,
not the subsatellite point
- - the attitude correction is not applied
- - the start time of the granules has to be more precise than the
second... I think the file name provide tenth of seconds, right?
- - no terrain correction is done

But tell me how it goes!
I started validating against HRPT data yesterday, and with all the
shortages mentionned above, I still got "only" 6.5km mean error
compared to the result of AAPP :)

> I have attached the improved map.

Really nice! I'm no basemap expert, so I'd gladly have a look at how
you did this :)

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJQj5J6AAoJEBdvyODiyJI4aakH/2j4k5Nt86cRKIYQzCjNB8FU
eWfnnJn9F9nrOhAB8aXGKIcYliaNoI3CR8kN4sRnuyY4krfL1yuuhtZU0haV8tPn
0//RMzIILKARo/B5XytSDj0Fddpvgb+gpSQRBotVlqpMCY1SGzL7Lcm48OJ320j1
8lE94o3w6dODmN0CYW60ex1R+jQZ8GcsGMTbyVE8UC6Aj9JbXxrid9t4oQwqCamU
LZZ+ESLmDRsFxdPVryg2EY0IGfzXyRwSLgbhXcg/Y4qR2ZfSFukeL5XGkBXlmZI8
nVub4ErZ8949l5ZxNuv6IE20cRPIuTyMbi9VzS2ohox6nqMSY7nminz1frGvWHk=
=+y3C
-----END PGP SIGNATURE-----

[Rutger]

Hey,

I dont take the tenth of a second into account right now, but i should indeed use it.

Regarding the sensor; VIIRS has a trick to minimize pixel deformation as a result of the bow-tie effect. From what i understand, the imagery resolution across track is 371m (so indeed 11.872km for the full array). But the resolution across scan is 95m, at nadir however the pixels are aggregated using a 3x1 ratio creating a resolution of 285 meter. Aggregation takes place only along the scan, so the upper and lower edges of the swath can be treated as normal i guess. The left and right edges might be impacted by this, because at the edge there is no aggregation performed.

Source: http://eospso.gsfc.nasa.gov/eos_homepage/mission_profiles/docs/NPP.pdf

Ill try to clear it up a bit more, perhaps plotting all geolocation points from the non-terrain corrected will reveal how they handle these aggregations.

Attached are some pictures about the NPP sensor.

Regards,
Rutger

[Rutger]

Hey,

> Really nice! I'm no basemap expert, so I'd gladly have a look at how
> you did this :)

As promised, ive attached the code im using. Its still a bit of a work in progress but it should be clean enough to see how im handling the plotting.
It basically comes down to setting up a Basemap object, and get it as an Matplotlib axes with plt.gc() and then add patches (polygons in this case) to it with ax.add_patch().

The creation of polygons is done with the functions makePatch() and geom2path().

I've made two npp_geometry functions, yours and an 'edge' only version i made. It wasnt really clear to me what direction the y_angle of the scanline was pointing. I figured a negative angle looks back with respect to the satellites direction, but with trial and error i find the opposite looking better (slight overlap due to the bowtie effect as opposed to a small gap between granules).

Let me know if you have any questions.

Regards,
Rutger

On Tuesday, October 30, 2012 9:40:28 AM UTC+1, Martin Raspaud wrote:

[Rutger]

Hey,

> I've made two npp_geometry functions, yours and an 'edge' only version i made. It wasnt really clear to me what direction the y_angle of the scanline was pointing.
> I figured a negative angle looks back with respect to the satellites direction, but with trial and error i find the opposite looking better (slight overlap due to the bowtie effect as opposed to a small gap between granules).

I made a mistake there, accidently reversed the entire array instead of the angle to be appended, it should be:
y_angle = np.tile(-y_max_angle, xpixels)
y_angle = np.append(y_angle, np.tile(0, scanline_nb-2))
y_angle = np.append(y_angle, np.tile(y_max_angle, xpixels))
y_angle = np.append(y_angle, np.tile(0, scanline_nb-2))


Regards,
Rutger

[Martin Raspaud]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 05/11/12 08:18, Rutger wrote:
> Hey,
>
>> I've made two npp_geometry functions, yours and an 'edge' only
>> version i
> made. It wasnt really clear to me what direction the y_angle of the
> scanline was pointing.
>> I figured a negative angle looks back with respect to the
>> satellites
> direction, but with trial and error i find the opposite looking
> better (slight overlap due to the bowtie effect as opposed to a
> small gap between granules).
>
> I made a mistake there, accidently reversed the entire array
> instead of the angle to be appended, it should be: y_angle =
> np.tile(-y_max_angle, xpixels) y_angle = np.append(y_angle,
> np.tile(0, scanline_nb-2)) y_angle = np.append(y_angle,
> np.tile(y_max_angle, xpixels)) y_angle = np.append(y_angle,
> np.tile(0, scanline_nb-2))

Hi Rutger,

Thanks for the correction! I was getting confused with you last mail... ;)

Keep me updated if you do some validation of the geolocations against
non-terrain corrected data.

We should also have a look at attitude correction probably...

Best regards,
Martin

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJQl6NmAAoJEBdvyODiyJI42kAIAIx5aYw2Lqa2/3/2kMPrZB0D
zknHLfvzMWgwjubkYQFWSwD/8XD2AAnvvKQb3TnDNRgh8qXir9lWGeP1V/MS+F2g
2NPO1hszm2KQsVEoKt+XBmajPkjV9XiwtDFgYPVTW7Zo24j5hz+O40z1v/J4wQ2G
atP0TEAngi0V0zsX7NEK4ZOL23/BrTJZ09c67Obvf6jSslwi5fuNgNsH+BMJ0tqa
u7wnI66mm2xxsnGjxM0imssDDAClsUI7E1GCA/+taTK3XWgU2XG1wvHsDl7NHab9
Mkhdw38sCZfow0F1Evql3ohkYF8K3xqWLqPWejrZsfTAZGH0FczBJWuxQvQGRFE=
=ELdM
-----END PGP SIGNATURE-----

[Rutger]

Hey,

> We should also have a look at attitude correction probably...

You have mentioned it before, i noticed you already added the pitch/roll/yaw parameters in the geoloc file. I guess the hard thing is getting those values? Im planning to read up on this stuff because its all very new for me. Information on getting the attitude isnt very widespread, or im looking in the wrong place. Would it be something which can be derived from the TLE (or SGP4 parameters)?

NASA also provides attitude parameters in the VIIRS hdf's, which might be nice for checking. They provide three per scanline.

The non-terrain corrected geolocation files are not produced anymore, at least not on the ladsweb ftp where im getting my data from. But because the files are almost equal i'll start of creating a workflow for the TC ones. Another detail of interest, every 31 or 32 granules there is one which only has 47 scanlines. There is a n_scanlines tag in the hdf and it can also be seen at the duration of the granule. I havent previously taken that into account.


Regards,
Rutger

[Martin Raspaud]

On 09/11/12 09:05, Rutger wrote:
> Hey,

Hi Rutger,

>

>> We should also have a look at attitude correction probably...
>
> You have mentioned it before, i noticed you already added the
> pitch/roll/yaw parameters in the geoloc file. I guess the hard thing is
> getting those values? Im planning to read up on this stuff because its all
> very new for me. Information on getting the attitude isnt very widespread,
> or im looking in the wrong place. Would it be something which can be
> derived from the TLE (or SGP4 parameters)?

The actual problem is to get attitude correction right. The code I added
is really straightforward, and if I look in the NWP AAPP documentation,
I see that there are already two different ways to do it with NOAA or
METOP satellites. So it needs to be checked.

The second point on estimating the attitude correction is a bit harder.
I know Thomas Lavergne had some ideas about that (Thomas, feel free to
chip in here), but if I remember correctly, the only good way to do it
is to detect landmarks on the ground and estimate the attitude from there.

I don't think it's possible to use SGP parameters to compute the
attitude, because the attitude is probably corrected quite often... But
I'm no expert here.

> NASA also provides attitude parameters in the VIIRS hdf's, which might be
> nice for checking. They provide three per scanline.
>
> The non-terrain corrected geolocation files are not produced anymore, at
> least not on the ladsweb ftp where im getting my data from. But because the
> files are almost equal i'll start of creating a workflow for the TC ones.
> Another detail of interest, every 31 or 32 granules there is one which only
> has 47 scanlines. There is a n_scanlines tag in the hdf and it can also be
> seen at the duration of the granule. I havent previously taken that into
> account.

I'll ask Adam, maybe we can generate non terrain corrected geolocations
in CSPP, so I can send you the file if you're willing to work on that.

Otherwise, I guess pixels over sea should not have that much terrain
correction, so indeed, you can try with that.

About the number of scanlines, yes, we should definitely read it from
the hdf file.

Best regards,
Martin