Calculationg Lat / Lon от геостационарной проекции - PullRequest
Новая
       61

Calculationg Lat / Lon от геостационарной проекции

1 голос
/ 27 сентября 2019

Я анализирую данные спутника 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)

минимальное и максимальное значения здесь составляют -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

Я все еще относительно новичок в манипулировании геопространственными данными.Что я тут не так делаю?

1 Ответ

0 голосов
/ 30 сентября 2019

Глядя на ваши данные, похоже, что они пересекают антимеридиан (-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)

Но, как только вы преобразовать данные: (примечание стороны, вот почему я использую трансформатор вместо 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)

вы увидите, что углы не все на той же сторонестрока: 

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

Из-за этого, когда вы делаете мин и макс, чтобы получить оценки, будет вводить в заблуждение, как мин и максимальной долготой все будет рядом с часть небесного меридиана, лежащая ниже горизонта (-180/ 180) для лат / долг. 

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

(-179.99994221231273, 14.564185533235145, 179.9999815862486, 53.5204809865526)
...