Не удается преобразовать проекцию базовой карты Matplotlib - PullRequest
0 голосов
/ 09 октября 2018

Цилиндрическая проекция Конформная коническая проекция Ламберта

Я пытаюсь построить спутниковое изображение GOES-16 и преобразовать его из проекции GEOS в проекцию LCC.Первоначально я смог преобразовать его в CYL без проблем, и он работал нормально, но, похоже, то же самое не работает для LCC ... Я получаю пустой график.

    nc = Dataset(fname)

      # Subset immediately to get rid of the points with missing values
      # since Python cannot handle them at this huge dimension
    data_subset = nc.variables['CMI'][:][1500:-1,2000:-1]
    #data_subset = nc.variables['CMI'][:]


    print("Running Loop")


    if first:
        bmap = Basemap(projection='cyl', llcrnrlon=-130, llcrnrlat=15, urcrnrlon=-70, urcrnrlat=50,  resolution='l')
        #bmap = Basemap(llcrnrlon=-130, llcrnrlat=15, urcrnrlon=-70, urcrnrlat=50, projection='lcc',  resolution='l')
        #bmap = Basemap(projection='lcc', resolution='c', lat_0=22, lon_0=-83, width=8E4, height=8E4)
        # Create the projection variables
        ori_proj = nc.variables['goes_imager_projection']
        sat_h = ori_proj.perspective_point_height
        sat_lon = ori_proj.longitude_of_projection_origin
        sat_sweep = ori_proj.sweep_angle_axis
        # The projection x and y coordinates equals
        # the scanning angle (in radians) multiplied by the satellite height (http://proj4.org/projections/geos.html)
        X = nc.variables['x'][:] * sat_h
        Y = nc.variables['y'][:] * sat_h
        p = Proj(proj='geos', h=sat_h, lon_0=sat_lon, sweep=sat_sweep)
        # Convert map points to latitude and longitude with the magic provided by Pyproj
        XX, YY = np.meshgrid(X, Y)
        lons, lats = p(XX, YY, inverse=True)
        lons_subset=lons[1500:-1,2000:-1]
        lats_subset=lats[1500:-1,2000:-1]
        levels=np.linspace(-0.1, 1.1, 1000)
        first= False 

    bmap.contourf(lons_subset,lats_subset,data_subset, levels=levels, cmap="gist_gray", extend='both')
    bmap.drawcoastlines(linewidth=0.5, linestyle='solid', color='white')
    bmap.drawcountries(linewidth=0.5, linestyle='solid', color='white')
    bmap.drawparallels(np.arange(-90.0, 90.0, 10.0), linewidth=0.1, color='white', labels=[True, False, False, True])
    bmap.drawmeridians(np.arange(0.0, 360.0, 10.0), linewidth=0.1, color='white', labels=[True, False, False, True])

    date_formatted = datetime.strftime(datetime_start,'%H:%MZ %a %d %b %Y')
    plt.title("GOES-16 ABI Radiance \n Scan from " +date_formatted)
    plt.savefig(image_string, figsize=(8,8), bbox_inches='tight', dpi=900)
    plt.clf()

СледующееПри запуске сценария для LCC получено предупреждение.

  /home/awips/anaconda3/envs/gdal/lib/python3.6/site-packages/matplotlib   /contour.py:1540: UserWarning: Warning: converting a masked element to nan.
  self.zmax = float(z.max())

  /home/awips/anaconda3/envs/gdal/lib/python3.6/site-packages/matplotlib /contour.py:1541: UserWarning: Warning: converting a masked element to nan.
  self.zmin = float(z.min())

1 Ответ

0 голосов
/ 10 июня 2019

Я бы порекомендовал использовать класс pyproj.Transformer:

>>> import pyproj
>>> pyproj.__version__
'2.2.0'
...
>>> crs_geos = pyproj.CRS(proj='geos', h=sat_h, lon_0=sat_lon, sweep=sat_sweep)
>>> crs_lcc = pyproj.CRS(proj='lcc', lat_0=22, lon_0=-83, lat_1=33, lat_2=45)
>>> transformer = pyproj.Transformer.from_crs(crs_geos, crs_lcc)
>>> lons, lats = transformer.transform(XX, YY)

Также вас может заинтересовать: http://pyproj4.github.io/pyproj/v2.2.0rel/api/crs.html#pyproj.crs.CRS.from_cf

Причина этого в том, что он обрабатывает базовые сдвиги,есть новые исправления, так как он использует PROJ 6 +.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...