GOES卫星产品没有提供经纬度，需要单独计算，这个代码提供了GOES-16卫星获取经纬度的方法uba/goes-latlon: Python script that can be used to generate latitude/longitude coordinates for GOES-16 full-disk extent. (github.com)
但该方法仅支持GOES-16卫星，笔者以GOES-18为例进行了扩展。经纬度获取是通过投影参数获得，知道其他卫星的投影参数（下面是GOES-18卫星）可以修改以下变量仿照扩展：
G16Proj4String = '+proj=geos +h=35786023.0 \
+a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 \
+lat_0=0.0 +lon_0=-137.0 +sweep=x +no_defs'扩展后的goesproj.py代码是：
# FROM: https://github.com/uba/goes-latlon/blob/main/goesproj.py
# MIT License
# Copyright (c) 2022 Douglas Uba
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
__author__ = 'Douglas Uba'
__email__ = 'douglas.uba@inpe.br'
from netCDF4 import Dataset
import numpy as np
import pyproj
# Define KM_PER_DEGREE
KM_PER_DEGREE = 40075.16/360.0
# GOES-16 viewing point (satellite position) height above the Earth
H = 35786023.0 # nf.variables['nominal_satellite_height'][:]
# GOES-16 Spatial Reference System (proj4 format)
G16Proj4String = '+proj=geos +h=35786023.0 \
+a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 \
+lat_0=0.0 +lon_0=-75.0 +sweep=x +no_defs'
# Lat/Lon WGS84 Spatial Reference System (proj4 string format)
LatLonWGS84Proj4String = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'
# GOES-16 Spatial Reference System
G16Proj = pyproj.Proj(G16Proj4String)
# Lat/Lon WSG84 Spatial Reference System
LatLonWGS84Proj = pyproj.Proj(LatLonWGS84Proj4String)
# Full-disk (FD) properties by resolution (GOES-16)
G16FDNLinesDic = {'0.5' : 21696, '1.0' : 10848, '2.0' : 5424, '4.0' : 2712, '10.0' : 1086}
G16FDNColsDic = {'0.5' : 21696, '1.0' : 10848, '2.0' : 5424, '4.0' : 2712, '10.0' : 1086}
G16FDScaleDic = {'0.5' : 0.000014, '1.0' : 0.000028, '2.0' : 0.000056, '4.0' : 0.000112, '10.0' : 0.000280}
G16FDOffsetDic = {'0.5' : 0.151865, '1.0' : 0.151858, '2.0' : 0.151844, '4.0' : 0.151816, '10.0' : 0.151900}
# update projection params
def update_projection_params(satellite='G16'):
if satellite == 'G16':
# pass # default
G16Proj4String = '+proj=geos +h=35786023.0 \
+a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 \
+lat_0=0.0 +lon_0=-75.0 +sweep=x +no_defs'
elif satellite == 'G18':
# nf.variables['goes_imager_projection']
G16Proj4String = '+proj=geos +h=35786023.0 \
+a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 \
+lat_0=0.0 +lon_0=-137.0 +sweep=x +no_defs'
else:
raise NotImplementedError
# GOES-16 Spatial Reference System
G16Proj = pyproj.Proj(G16Proj4String)
return G16Proj, G16Proj4String
def GetFullDiskInfos(resolution):
return G16FDNLinesDic[resolution], G16FDNColsDic[resolution], \
G16FDScaleDic[resolution], G16FDOffsetDic[resolution]
def BuildLatLonGrid(resolution, path=None, satellite='G16'):
# update projection params by satellite
G16Proj, G16Proj4String = update_projection_params(satellite=satellite)
# Get full-disk properties
nlines, ncols, scale, offset = GetFullDiskInfos(resolution)
# Create satellite-projection coordinates array
x = np.arange(0, ncols)
y = np.arange(0, nlines)
# Apply scale and offset (scanning angle) + H factor = projection coordinates
x = ((x * scale) - offset) * H
y = ((y * (-1 * scale)) + offset) * H
# Create matrix with values
x, y = np.meshgrid(x, y)
# Reshape to vector (1-d)
x = x.reshape(nlines * ncols)
y = y.reshape(nlines * ncols)
# Transform
transformer = pyproj.Transformer.from_proj(G16Proj, LatLonWGS84Proj) # source -> from
lon, lat = transformer.transform(x, y)
# Reshape to matrix format
lat = lat.reshape(nlines, ncols)
lon = lon.reshape(nlines, ncols)
lat = np.ma.masked_invalid(lat)
lon = np.ma.masked_invalid(lon)
# Export [lat, lon] to file, if requested
if(path):
export2file(lat, lon, path)
return lat, lon
def getScaleOffset(values, n):
# From: http://james.hiebert.name/blog/work/2015/04/18/NetCDF-Scale-Factors.html
# stretch/compress data to the available packed range
max = np.max(values); min = np.min(values)
scale = (max - min)/(2**n - 1)
# translate the range to be symmetric about zero
offset = min + 2**(n - 1) * scale
return (scale, offset)
def pack(values, scale, offset):
return np.floor((values - offset)/scale)
def export2file(lat, lon, path):
nc = Dataset(path, 'w', format='NETCDF4');
nc.createDimension('size', lat.shape[0])
latvar = nc.createVariable('lat', 'i2', ('size', 'size'), zlib=True) # i2 -> int16
lonvar = nc.createVariable('lon', 'i2', ('size', 'size'), zlib=True) # i2 -> int16
# Pack data using scale and offset
lat_scale, lat_offset = getScaleOffset(lat, 16) # 16 -> int16 = i2
lon_scale, lon_offset = getScaleOffset(lon, 16) # 16 -> int16 = i2
lat_packed = pack(lat, lat_scale, lat_offset)
lon_packed = pack(lon, lon_scale, lon_offset)
latvar[:] = lat_packed.astype(np.int16)
lonvar[:] = lon_packed.astype(np.int16)
# Adjust meta-data info
latvar.scale_factor = lat_scale; latvar.add_offset = lat_offset
lonvar.scale_factor = lon_scale; lonvar.add_offset = lon_offset
nc.close()
