ABI fixed grid coordinate values

Post any questions, ideas, or topics related to Jython and Python scripting.
Post Reply
User avatar
joleenf
Posts: 1123
Joined: Mon Jan 19, 2009 7:16 pm

ABI fixed grid coordinate values

Post by joleenf »

Is there a way to access the satellite viewing angle of the ABI data or the actual x and y coordinate values?

Thanks,
Joleen
User avatar
bobc
Posts: 988
Joined: Mon Nov 15, 2010 5:57 pm

Re: ABI fixed grid coordinate values

Post by bobc »

Hi Joleen,

I looked into getting at satellite viewing angles in McV. At this point, it appears that the work done in SunRelativePosition.java returns only solar-related angles. I wrote up a script that extracts the angles at each pixel of a local M1 netCDF file, though I haven't verified the accuracy to this point. Basically, the script does the following:

  1. Loads in/displays data from a local netCDF ABI M1 file. I set a pretty high stride for testing purposes.
  2. Obtains the lat/lon of each pixel.
  3. Figures out the distance of each pixel (from step 2 above) from the satellite subpoint.
  4. Determines the satellite viewing angle at each pixel by taking the arctangent of the pixel's distance from the subpoint (from step 3 above) divided by the height of the satellite.
  5. Prints out the lat/lon/sat viewing angle at each pixel.

The script is below. You'd have to change the dataDir and dataFile variables at the beginning of the script to point at your data.

Code: Select all

# setup
import math
from math import sin, cos, sqrt, atan2, radians
dataDir = 'V:/rcarp/Data/ABI/eclipse/M1/'
dataFile = 'OR_ABI-L1b-RadM1-M3C01_G16_s20172331738268_e20172331738325_c20172331738368.nc'

# dictionary for loadGrid
dataParms = dict(
    filename = '%s%s' % (dataDir, dataFile),
    field = 'Rad',
    stride = 50,
    )

# create data object with loadGrid and display data
data = loadGrid(**dataParms)
layer = activeDisplay().createLayer('Image Display', data)

# Use JPythonMethods to extract lists of lat/lons of
# each data point in the data object
ds = getDomainSet(data)
latLons = getLatLons(ds)
lats = latLons[0]
lons = latLons[1]

# This is curtosy of:
# https://stackoverflow.com/questions/19412462/
# getting-distance-between-two-points-based-on-latitude-longitude
# We're figuring out the distance between each data pixel and
# the subpoint of the satellite
distList = []
numPts = len(lats)
R = 6373.0
for x in range(0, numPts):
    # lat1/lon1 are the pixels of data
    lat1 = radians(lats[x])
    lon1 = radians(lons[x])

    # lat2/lon2 are for the satellite subpoint
    lat2 = radians(0)
    lon2 = radians(-75)

    dlat = lat2 - lat1
    dlon = lon2 - lon1

    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
    c = 2 * atan2(sqrt(a), sqrt(1 - a))

    distance = R * c
    distList.append(distance)

# height of satellite in km
satHeight = 35786

# determine the satellite viewing angle by taking the arctangent
# of the opoosite (from distList) divided by satHeight
satAngles = []
for x in range(0, numPts):
    ang = math.degrees(math.atan(distList[x] / satHeight))
    satAngles.append(ang)

# print out lat/lon/sat viewing angles at each pixel
for x in range(0, numPts):
    lat = lats[x]
    lon = lons[x]
    ang = satAngles[x]
    print 'At %s, %s (lat,lon) the viewing angle is %s degrees' % (lat, lon, ang)

Please let me know if you have any questions about this.

Thanks,
Bob Carp
User avatar
bobc
Posts: 988
Joined: Mon Nov 15, 2010 5:57 pm

Re: ABI fixed grid coordinate values

Post by bobc »

As an addendum to my last post, I thought of a better way to specify the subpoint lat/lon. With how GOES-16 has changed positions a couple times since it began sending data, you might not want to have hard-coded subpoint values in the script. You can make use of the nominal_satellite_subpoint_lat and noniminal_satellite_subpoint variables contained in the data:

Code: Select all

# This is curtosy of:
# https://stackoverflow.com/questions/19412462/
# getting-distance-between-two-points-based-on-latitude-longitude
# We're figuring out the distance between each data pixel and
# the subpoint of the satellite
distList = []
numPts = len(lats)
R = 6373.0

# extract subpoint info from metadata
satLat = data['metadataVariables']['nominal_satellite_subpoint_lat']['value']
satLon = data['metadataVariables']['nominal_satellite_subpoint_lon']['value']

for x in range(0, numPts):
    # lat1/lon1 are the pixels of data
    lat1 = radians(lats[x])
    lon1 = radians(lons[x])

    # lat2/lon2 are for the satellite subpoint
    lat2 = radians(float(satLat))
    lon2 = radians(float(satLon))

The lines that changes were:
  • The satLat= and satLon= lines before the for loop were added
  • the lat2= and lon2= lines in the for loop were modified to use the satLat/satLon variables from the previous point

Bob
User avatar
joleenf
Posts: 1123
Joined: Mon Jan 19, 2009 7:16 pm

Re: ABI fixed grid coordinate values

Post by joleenf »

Hi Bob,

Thank-you for the detailed response.

Joleen
Post Reply