1

我正在分析 GOES-16 卫星数据，需要计算每个网格坐标的基本纬度、经度作为进一步分析的一部分。我目前正在尝试使用 pyproj 来执行此操作，并且遇到图像预期范围之外的经度值，经度值在预期范围内。

下载图像：

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')

转换投影：

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)

这里的最小和最大 lon 值分别是 -179.99 和 179.9。我只希望这张图片包含 -89.6 和 175.6 之间的数据，如下所示。

当我检查此图像的预期地理范围时，我得到以下信息：

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

我对地理空间数据操作还比较陌生。我在这里做错了什么？

4

2 回答 2

0

使用我创建的GOES包，您可以轻松获得 ABI 像素的 Lon/Lat。本教程中显示了许多示例。希望对你有帮助。

于 2021-02-23T09:10:52.640 回答
0

查看您的数据，它看起来像是穿过了反子午线（-180/180 经度）。使用rioxarray检查数据。

通常的检查可能不会揭示这一点：

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)

但是，一旦你转换数据：（旁注，这就是为什么我使用 Transformer 而不是 Proj https://pyproj4.github.io/pyproj/stable/gotchas.html#proj-not-a-generic-latitude-经度投影转换器）

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)

您会看到角并不都在直线的同一侧：

(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))

因此，当您执行 min 和 max 来获得边界时，会产生误导，因为对于纬度/经度，最小和最大经度都将位于反子午线 (-180/180) 旁边。

lon.min(), lat.min(), lon.max(), lat.max()

(-179.99994221231273, 14.564185533235145, 179.9999815862486, 53.5204809865526)
于 2019-09-30T00:54:37.273 回答