Python中的Globcolour数据和投影错误

Question

我在显示来自Globcolour(1)，由于图像的matplotlib和cartopy定义使用的投影。

我下载了一个NetCDF格式的总悬浮物图像（这里是数据）。请在此输入链接描述)，当我试图将其与cartopy包中的海岸线一起显示时，海岸线和数据之间有一个众所周知的差距。正如你所看到的那样，像素应该在海岸线旁边（黑线），而不是超过陆地（旗帜图像中的黄色像素）。

这不应该发生。我使用QGIS检查并直接加载netcdf文件，海岸线设置正确。

最初我使用了PlateeCarrer投影的图像，考虑到如果图像是在WGS84中，它们会匹配，但显然它们不匹配。我试过使用matplotlib函数中的变换选项，但没有成功。要么是差距依然存在，要么是图的坐标变成了投影坐标，我的数据（是地理坐标）消失了。

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数据框架的形式接收数据。打开NetCDF的方法是 xarray.open_dataset 然后将其转化为熊猫与 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

会不会影响到？

Answer 1

我们回答费利佩通过电子邮件，我复制粘贴在这里。

一个小的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')

我认为你正面临2个问题。

1）我们的产品可能会重叠一些陆地区域，因为在GlobColour处理过程中的Level-3 rebinning：如果一个4公里的像素只有水面上的角落，我们将填补完整的像素。我们保留这些像素是因为它们可能对某些需求有用（例如陆地水位变化的区域），但在质量标志中，我们提供了一个LAND掩码，可以用来删除这些像素。如果你愿意，你也可以使用你自己的LAND掩码。下面的Python例子展示了如何使用LAND掩码。

2) 我怀疑你的Python代码引入了至少半个像素的东西南移，也许是因为latlon数组是针对每个像素的中心，但cartopy需要的范围是外部限制。

GlobColour标志在产品用户指南中定义了 http:/www.globcolour.infoCDR_DocsGlobCOLOUR_PUG.pdf 第76页。

GlobColour团队

Answer 2

你确定你的数据是WGS84的吗？在元数据中，我只看到：

'earth_radius': 6378.137

我的意思是假设一个半径为6378. 137公里的球形地球。我无法访问你的数据，但我会尝试设置一个新的数据。cartopy.crs.Globe 该半径的实例。