Python 中的 Globcolour 数据和投影错误
Globcolour data and projection error in Python
我在显示来自 Globcolour (1) 的一些数据时遇到问题，这是由于 matplotlib 和图像的 cartopy 定义使用了投影。
我下载了 NetCDF 格式的总悬浮物图像（这里是数据 enter link description here），当我尝试显示它时，连同 cartopy 包中的海岸线，两者之间存在臭名昭著的差距海岸线和数据。正如你在下面看到的，像素应该靠近海岸线（黑线），而不是进入陆地（旗帜图像中的黄色像素）
这不应该发生。我检查使用 QGIS 并直接加载海岸线设置正确的 netcdf 文件。
最初我对图像使用了 PlateeCarrer 投影，考虑到如果图像在 WGS84 中，它们会匹配，但显然它们不匹配。我试过在 matplotlib 函数中使用 transform 选项，但没有成功。差距仍然存在，或者图形的坐标变为投影坐标并且我的数据（在地理坐标中）消失了。
NetCDF 文件的属性是：
'grid_type': 'Equirectangular',
'spatial_resolution': 4.6383123,
'nb_equ_bins': 55,
'registration': 5,
'lat_step': 0.041666668,
'lon_step': 0.041666668,
'earth_radius': 6378.137,
'max_north_grid': 11.124998,
'max_south_grid': 9.27,
'max_west_grid': -86.25,
'max_east_grid': -83.97,
'northernmost_latitude': 11.124998,
'southernmost_latitude': 9.249998,
'westernmost_longitude': -86.25,
'easternmost_longitude': -84.0,
'nb_grid_bins': 2475,
'nb_bins': 2475,
'pct_bins': 100.0,
'nb_valid_bins': 1089,
'pct_valid_bins': 44.0,
'netcdf_version': '4.3.3.1 of Jul 8 2016 18:15:50 $',
'DPM_reference': 'GC-UD-ACRI-PUG',
'IODD_reference': 'GC-UD-ACRI-PUG'}
我用来绘制图像的代码是：
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import cartopy.crs as ccrs
import dill as pickel
def paint_maps(df_std=None, fecha=1, attributes=None,
savefol='/media/felipe/TOSHIBA EXT/iMARES/Investigacion/2019_MariculturaPacifico/DB/figures/',
disp_fig=0):
"""Función para dibujar los datos contenidos en los archivos netCDF de SST, Salinidad y propiedad ópticas del agua.
Recibe el dataframe con la información en formato de Pandas Dataframe, y selecciona según una fecha establecida,
el conjunto de datos con coordenadas Lat-Lon que debe dibujar. Esos los dibuja y transforma a formato raster. Unido
se dibuja también la línea de costa proveniente de un archivo shapefile. La función dibuja toda la información
contenida en el dataframe aportado (datos, anomalías, flags, y cualquier otro dato que tenga.
Recibe:
df_std: dataframe con la información a dibujar. Debe venir indexado por fecha, lat y lon.
fecha: día que se elige dibujar. Formato string 'yyyymmdd'. Valor 1 significa que grafica el valor promedio de todas las fechas en cada
píxel. Promedio simple ignorando NaN's
attributes: diccionario con los atributos del netcdf de donde se obtiene nombre de variable y unidades. Creado
con open_netcdf.py
savefol: carpeta donde se guardan las imágenes dibujadas
disp_fig: booleano para imprimir figura en pantalla.
Devuelve:
Nada. Solo crea y guarda figuras"""
# Identifica la fecha solicitada (cuando se ha especificado) y confirma que sea parte del registro. Extrae la
# información del Dataframe en la fecha que se solicitó, o calcula el promedio de todas las fechas para graficar
# el valor promedio.
if fecha != 1:
if isinstance(fecha, str):
fecha = pd.to_datetime(fecha + '120000')
else:
print('La fecha indicada no está en formato String. Reinicie la ejecución.')
try:
idx = pd.IndexSlice
df_map = df_std.loc[idx[:, :, fecha], :]
except:
print('Se generó un error. Posiblemente fecha no está dentro del registro. La fecha debe estar entre el ' + df_std.index[0][-1].strftime('%d/%m/%Y') + ' y el ' + df_std.index[-1][-1].strftime('%d/%m/%Y'))
raise
else:
df_map = df_std.groupby(['lat', 'lon']).mean()
# Reestructura la información para tenerla en forma de matriz y dibujarla de forma más simple. Extrae los valores y
# las latitudes y longitudes correspondientes, así como los valores de la variable y sus flags.
df_map2 = df_map.unstack(level=0)
vari = df_map2['mean_val'].values
flags = df_map2['flag_val'].values
lat = df_map2['mean_val'].columns.get_level_values('lat')
lon = df_map2['mean_val'].index.get_level_values('lon')
# Extrae de los atributos del netcdf el nombre de la variable a graficar y las unidades
variable_str = attributes['variable']['long_name']
variable_units = attributes['variable']['units']
# Dibuja el mapa que se haya seleccionado según fecha (valor promedio del valor o fecha específica)
fig, ax = plt.subplots(1, 2, figsize=(10, 10), subplot_kw={'projection': ccrs.PlateCarree()})
extend = [lon[1], lon[-1], lat[1], lat[-1]]
# Primera figura. Variable a graficar. Usa línea de costa del cartopy y coloca una leyenda abajo
ax[0].set_extent(extend)
ax[0].coastlines(resolution='10m')
#cs = ax[0].pcolormesh(lon, lat, vari.T)
cs = ax[0].pcolormesh(lon, lat, vari.T, transform=ccrs.PlateCarree())
ax[0].set_title(variable_str)
cax, kw = matplotlib.colorbar.make_axes(ax[0], location='bottom', pad=0.05, shrink=0.7)
out = fig.colorbar(cs, cax=cax, extend='both', **kw)
out.set_label('Units: '+variable_units, size=10)
# Segunda figura. Flags de la figura. Usa la leyenda directamente de los datos usados.
ax[1].set_extent(extend)
ax[1].coastlines(resolution='10m')
cs2 = ax[1].pcolormesh(lon, lat, flags.T)
ax[1].set_title('Flags')
cax, kw = matplotlib.colorbar.make_axes(ax[1], location='bottom', pad=0.05, shrink=0.7)
out = fig.colorbar(cs2, cax=cax, extend='both', **kw)
out.set_label('Flags', size=10)
# Salva la figura
plt.savefig(savefol+variable_str+'.jpg', bbox_inches='tight')
with open(savefol+'fig_'+variable_str+'.pickel', 'wb') as f:
pickel.dump(fig, f)
# Imprime figura si se elige opción con disp_fig
if disp_fig == 1:
plt.show()
return
它接收数据作为 Pandas 数据帧。使用 xarray.open_dataset 打开 NetCDF，然后使用 to_dataframe()
将其转换为 Pandas
我在 Ubuntu 中使用 Python 3.7。
最后一件事。加载cartopy.crs包时出现这个错误：
ERROR 1: PROJ: proj_create_from_database: Open of /home/felipe/anaconda3/envs/personal/share/proj failed
会影响吗？
您确定您的数据采用 WGS84 格式吗？查看元数据，我只看到：
'earth_radius': 6378.137
我的意思是假设一个半径为 6378.137 公里的球形地球。我无权访问您的数据，但我会尝试设置一个具有该半径的 cartopy.crs.Globe 实例。
我们通过电子邮件回复了 Felipe，我 copy/paste 在这里：
一个小 Python 脚本，用于从 TSM GlobColour 产品创建您所在区域的地图（我使用月度产品以获得良好的覆盖范围）：
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw=dict(projection=ccrs.PlateCarree()))
# my region of interest
ax.set_extent([-86, -84, 9, 11])
ax.coastlines(resolution='10m', color='red')
nc_dst = nc.Dataset('L3m_20100101-20100131__GLOB_4_AV-MER_TSM_MO_00.nc')
# extent of the product
data_extent = [nc_dst.max_west_grid, nc_dst.max_east_grid,
nc_dst.max_south_grid, nc_dst.max_north_grid]
data = nc_dst.variables['TSM_mean'][:]
flags = nc_dst.variables['TSM_flags'][:]
land = flags & 8 # LAND == 3rd bit == 2^3 == 8
data_noland = np.ma.masked_where(land, data)
ax.imshow(data_noland, origin='upper', extent=data_extent)
plt.savefig('TSM_noland.png')
ax.imshow(data, origin='upper', extent=data_extent)
plt.savefig('TSM.png')
我认为您面临两个问题：
1) 由于 GlobColour 处理过程中的 Level-3 重新分箱，我们的产品可能会与某些陆地区域重叠：如果一个 4km 像素只有水面上的一个角，我们将填充整个像素。我们保留它们是因为它们可能对某些需求有用（例如 land/water 限制变化的区域），但在质量标志中我们提供了一个 LAND 掩码，可用于删除这些像素。如果愿意，您也可以使用自己的 LAND 掩码。下面的 Python 示例显示了如何使用 LAND 掩码。
2) 我怀疑您的 Python 代码引入了至少半个像素的 east/south 偏移，这可能是因为 lat/lon 数组用于每个像素的中心，但范围cartopy需要的是外部限制。
GlobColour 标志在产品用户指南 http://www.globcolour.info/CDR_Docs/GlobCOLOUR_PUG.pdf 第 76 页中定义。
GlobColour 团队
我在显示来自 Globcolour (1) 的一些数据时遇到问题，这是由于 matplotlib 和图像的 cartopy 定义使用了投影。
我下载了 NetCDF 格式的总悬浮物图像（这里是数据 enter link description here），当我尝试显示它时，连同 cartopy 包中的海岸线，两者之间存在臭名昭著的差距海岸线和数据。正如你在下面看到的，像素应该靠近海岸线（黑线），而不是进入陆地（旗帜图像中的黄色像素）
这不应该发生。我检查使用 QGIS 并直接加载海岸线设置正确的 netcdf 文件。
最初我对图像使用了 PlateeCarrer 投影，考虑到如果图像在 WGS84 中，它们会匹配，但显然它们不匹配。我试过在 matplotlib 函数中使用 transform 选项，但没有成功。差距仍然存在，或者图形的坐标变为投影坐标并且我的数据（在地理坐标中）消失了。
NetCDF 文件的属性是：
'grid_type': 'Equirectangular',
'spatial_resolution': 4.6383123,
'nb_equ_bins': 55,
'registration': 5,
'lat_step': 0.041666668,
'lon_step': 0.041666668,
'earth_radius': 6378.137,
'max_north_grid': 11.124998,
'max_south_grid': 9.27,
'max_west_grid': -86.25,
'max_east_grid': -83.97,
'northernmost_latitude': 11.124998,
'southernmost_latitude': 9.249998,
'westernmost_longitude': -86.25,
'easternmost_longitude': -84.0,
'nb_grid_bins': 2475,
'nb_bins': 2475,
'pct_bins': 100.0,
'nb_valid_bins': 1089,
'pct_valid_bins': 44.0,
'netcdf_version': '4.3.3.1 of Jul 8 2016 18:15:50 $',
'DPM_reference': 'GC-UD-ACRI-PUG',
'IODD_reference': 'GC-UD-ACRI-PUG'}
我用来绘制图像的代码是：
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import cartopy.crs as ccrs
import dill as pickel
def paint_maps(df_std=None, fecha=1, attributes=None,
savefol='/media/felipe/TOSHIBA EXT/iMARES/Investigacion/2019_MariculturaPacifico/DB/figures/',
disp_fig=0):
"""Función para dibujar los datos contenidos en los archivos netCDF de SST, Salinidad y propiedad ópticas del agua.
Recibe el dataframe con la información en formato de Pandas Dataframe, y selecciona según una fecha establecida,
el conjunto de datos con coordenadas Lat-Lon que debe dibujar. Esos los dibuja y transforma a formato raster. Unido
se dibuja también la línea de costa proveniente de un archivo shapefile. La función dibuja toda la información
contenida en el dataframe aportado (datos, anomalías, flags, y cualquier otro dato que tenga.
Recibe:
df_std: dataframe con la información a dibujar. Debe venir indexado por fecha, lat y lon.
fecha: día que se elige dibujar. Formato string 'yyyymmdd'. Valor 1 significa que grafica el valor promedio de todas las fechas en cada
píxel. Promedio simple ignorando NaN's
attributes: diccionario con los atributos del netcdf de donde se obtiene nombre de variable y unidades. Creado
con open_netcdf.py
savefol: carpeta donde se guardan las imágenes dibujadas
disp_fig: booleano para imprimir figura en pantalla.
Devuelve:
Nada. Solo crea y guarda figuras"""
# Identifica la fecha solicitada (cuando se ha especificado) y confirma que sea parte del registro. Extrae la
# información del Dataframe en la fecha que se solicitó, o calcula el promedio de todas las fechas para graficar
# el valor promedio.
if fecha != 1:
if isinstance(fecha, str):
fecha = pd.to_datetime(fecha + '120000')
else:
print('La fecha indicada no está en formato String. Reinicie la ejecución.')
try:
idx = pd.IndexSlice
df_map = df_std.loc[idx[:, :, fecha], :]
except:
print('Se generó un error. Posiblemente fecha no está dentro del registro. La fecha debe estar entre el ' + df_std.index[0][-1].strftime('%d/%m/%Y') + ' y el ' + df_std.index[-1][-1].strftime('%d/%m/%Y'))
raise
else:
df_map = df_std.groupby(['lat', 'lon']).mean()
# Reestructura la información para tenerla en forma de matriz y dibujarla de forma más simple. Extrae los valores y
# las latitudes y longitudes correspondientes, así como los valores de la variable y sus flags.
df_map2 = df_map.unstack(level=0)
vari = df_map2['mean_val'].values
flags = df_map2['flag_val'].values
lat = df_map2['mean_val'].columns.get_level_values('lat')
lon = df_map2['mean_val'].index.get_level_values('lon')
# Extrae de los atributos del netcdf el nombre de la variable a graficar y las unidades
variable_str = attributes['variable']['long_name']
variable_units = attributes['variable']['units']
# Dibuja el mapa que se haya seleccionado según fecha (valor promedio del valor o fecha específica)
fig, ax = plt.subplots(1, 2, figsize=(10, 10), subplot_kw={'projection': ccrs.PlateCarree()})
extend = [lon[1], lon[-1], lat[1], lat[-1]]
# Primera figura. Variable a graficar. Usa línea de costa del cartopy y coloca una leyenda abajo
ax[0].set_extent(extend)
ax[0].coastlines(resolution='10m')
#cs = ax[0].pcolormesh(lon, lat, vari.T)
cs = ax[0].pcolormesh(lon, lat, vari.T, transform=ccrs.PlateCarree())
ax[0].set_title(variable_str)
cax, kw = matplotlib.colorbar.make_axes(ax[0], location='bottom', pad=0.05, shrink=0.7)
out = fig.colorbar(cs, cax=cax, extend='both', **kw)
out.set_label('Units: '+variable_units, size=10)
# Segunda figura. Flags de la figura. Usa la leyenda directamente de los datos usados.
ax[1].set_extent(extend)
ax[1].coastlines(resolution='10m')
cs2 = ax[1].pcolormesh(lon, lat, flags.T)
ax[1].set_title('Flags')
cax, kw = matplotlib.colorbar.make_axes(ax[1], location='bottom', pad=0.05, shrink=0.7)
out = fig.colorbar(cs2, cax=cax, extend='both', **kw)
out.set_label('Flags', size=10)
# Salva la figura
plt.savefig(savefol+variable_str+'.jpg', bbox_inches='tight')
with open(savefol+'fig_'+variable_str+'.pickel', 'wb') as f:
pickel.dump(fig, f)
# Imprime figura si se elige opción con disp_fig
if disp_fig == 1:
plt.show()
return
它接收数据作为 Pandas 数据帧。使用 xarray.open_dataset 打开 NetCDF，然后使用 to_dataframe()
我在 Ubuntu 中使用 Python 3.7。
最后一件事。加载cartopy.crs包时出现这个错误：
ERROR 1: PROJ: proj_create_from_database: Open of /home/felipe/anaconda3/envs/personal/share/proj failed
会影响吗？
您确定您的数据采用 WGS84 格式吗？查看元数据，我只看到：
'earth_radius': 6378.137
我的意思是假设一个半径为 6378.137 公里的球形地球。我无权访问您的数据，但我会尝试设置一个具有该半径的 cartopy.crs.Globe 实例。
我们通过电子邮件回复了 Felipe，我 copy/paste 在这里：
一个小 Python 脚本，用于从 TSM GlobColour 产品创建您所在区域的地图（我使用月度产品以获得良好的覆盖范围）：
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig, ax = plt.subplots(figsize=(5, 5), subplot_kw=dict(projection=ccrs.PlateCarree()))
# my region of interest
ax.set_extent([-86, -84, 9, 11])
ax.coastlines(resolution='10m', color='red')
nc_dst = nc.Dataset('L3m_20100101-20100131__GLOB_4_AV-MER_TSM_MO_00.nc')
# extent of the product
data_extent = [nc_dst.max_west_grid, nc_dst.max_east_grid,
nc_dst.max_south_grid, nc_dst.max_north_grid]
data = nc_dst.variables['TSM_mean'][:]
flags = nc_dst.variables['TSM_flags'][:]
land = flags & 8 # LAND == 3rd bit == 2^3 == 8
data_noland = np.ma.masked_where(land, data)
ax.imshow(data_noland, origin='upper', extent=data_extent)
plt.savefig('TSM_noland.png')
ax.imshow(data, origin='upper', extent=data_extent)
plt.savefig('TSM.png')
我认为您面临两个问题：
1) 由于 GlobColour 处理过程中的 Level-3 重新分箱，我们的产品可能会与某些陆地区域重叠：如果一个 4km 像素只有水面上的一个角，我们将填充整个像素。我们保留它们是因为它们可能对某些需求有用（例如 land/water 限制变化的区域），但在质量标志中我们提供了一个 LAND 掩码，可用于删除这些像素。如果愿意，您也可以使用自己的 LAND 掩码。下面的 Python 示例显示了如何使用 LAND 掩码。
2) 我怀疑您的 Python 代码引入了至少半个像素的 east/south 偏移，这可能是因为 lat/lon 数组用于每个像素的中心，但范围cartopy需要的是外部限制。
GlobColour 标志在产品用户指南 http://www.globcolour.info/CDR_Docs/GlobCOLOUR_PUG.pdf 第 76 页中定义。
GlobColour 团队