1

圆柱投影 Lambert Conformal Conic Projection

我正在尝试绘制 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())
4

1 回答 1

1

我建议使用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+ 时有新的修复。

于 2019-06-10T01:22:45.927 回答