1

I am analyzing GOES-16 Satelite data and need to calculate the underlying lat, lon of each grid coordinate as part of further analysis. I am currently attempting to use pyproj to do this and am running into longitude values outside of the expected extent for the image, the longitude values are within the expected extent.

Downloading the image:

import boto3
import xarray as xr
s3 = boto3.client('s3')
s3.download_file('noaa-goes17', 'ABI-L1b-RadC/2019/022/21/OR_ABI-L1b-RadC-M3C02_G17_s20190222102189_e20190222104562_c20190222104588.nc', 
'test.nc')

ds = xr.open_dataset('test.nc')

Converting the projection:

from pyproj import Proj
p = Proj(proj='geos', h='35786023.0', lon_0='-137.0', sweep='x')

dat = ds.metpy.parse_cf('Rad')
cord_grid = np.meshgrid(dat['x'].values, dat['y'].values)
lons, lats = p(cord_grid[0], cord_grid[1], inverse=True)

the min and max lon value here is -179.99 and 179.9 here respectively. I would only expect this image to contain data between -89.6 and 175.6 as shown below.

When I check the expected geographical extent for this image I get the following:

ds['geospatial_lat_lon_extent']


Out[61]:
<xarray.DataArray 'geospatial_lat_lon_extent' ()>
array(9.96921e+36, dtype=float32)
Coordinates:
    t        datetime64[ns] ...
    y_image  float32 ...
    x_image  float32 ...
Attributes:
    long_name:                       geospatial latitude and longitude refere...
    geospatial_westbound_longitude:  175.62358
    geospatial_northbound_latitude:  53.50006
    geospatial_eastbound_longitude:  -89.62357
    geospatial_southbound_latitude:  14.57134
    geospatial_lat_center:           29.967
    geospatial_lon_center:           -137.0
    geospatial_lat_nadir:            0.0
    geospatial_lon_nadir:            -137.0
    geospatial_lat_units:            degrees_north
    geospatial_lon_units:            degrees_east

I am still relatively new to geospatial data manipulation. What am I doing wrong here?

skuzmier
  • 23
  • 5

2 Answers2

0

Looking at your data, it looks like it crosses the antimeridian (-180/180 longitude). Using rioxarray to inspect the data.

The usual inspection may not reveal this:

import numpy as np
from pyproj import Transformer
import rioxarray

rds = rioxarray.open_rasterio("test.nc", masked=True)
<xarray.Dataset>
Dimensions:      (band: 1, x: 10000, y: 6000)
Coordinates:
  * y            (y) float64 1.583e+06 1.584e+06 ... 4.588e+06 4.589e+06
  * x            (x) float64 -2.505e+06 -2.504e+06 ... 2.504e+06 2.505e+06
  * band         (band) int64 1
    spatial_ref  int64 0
Data variables:
    Rad          (band, y, x) float64 ...
    DQF          (band, y, x) float64 ...
rds.rio.crs
CRS.from_wkt('PROJCS["unnamed",GEOGCS["unknown",DATUM["unnamed",SPHEROID["Spheroid",6378137,298.2572221]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Geostationary_Satellite"],PARAMETER["central_meridian",-137],PARAMETER["satellite_height",35786023],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=geos +lon_0=-137 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs +sweep=x"]]')
rds.Rad.rio.transform_bounds("epsg:4326")
(-177.77133129663548,
 14.561801004390198,
 179.99792429593583,
 53.52729472471354)

But, once you transform the data: (Side note, here is why I use the Transformer instead of Proj https://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-longitude-to-projection-converter)

transformer = Transformer.from_crs(rds.rio.crs, "epsg:4326", always_xy=True)
x_coords, y_coords = np.meshgrid(
    rds.coords[rds.rio.x_dim], rds.coords[rds.rio.y_dim]
)
lon, lat = transformer.transform(x_coords, y_coords)

You will see that the corners are not all on the same side of the line:

(lon[0][0], lat[0][0]), (lon[-1][-1], lat[-1][-1]), (lon[0][-1], lat[0][-1]), (lon[-1][0], lat[-1][0])```
((-161.57648657675344, 14.798043104222199),
 (-89.56850957728933, 53.5204809865526),
 (-112.4235101487757, 14.798043165552787),
 (175.56851955360526, 53.52047990993389))

Due to this, when you do the min and max to get the bounds, it will be misleading as the min and max longitude will all be next to the antimeridian (-180/180) for lat/lon.

lon.min(), lat.min(), lon.max(), lat.max()
(-179.99994221231273, 14.564185533235145, 179.9999815862486, 53.5204809865526)
snowman2
  • 646
  • 4
  • 11
0

With the GOES package that I created, you can get the Lon/Lat of ABI pixels easily. Many examples are shown in this tutorial. I hope you find it useful.

joaohenry23
  • 145
  • 5