Python Apex.convert方法代码示例

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def func4():
    """
    Test (lat, lt) to (mlat, mlt) conversions.

    1, Use a = Apex(date=...); "date" determines which IGRF coefficients are
       used in conversions. Uses current date as default.
    2, height is needed for better conversion.
       champ and grace files use qd coordinates.
    3, mlt in qd and apex coordinates are the same.
    """
    import champ_grace as cg
    from apexpy import Apex as Apex
    import matplotlib.pyplot as plt
    a = cg.ChampDensity('2005-1-1', '2005-1-2')
    b = Apex(date=2005)
    mlatt, mlt = b.convert(
            lat=a.lat, lon=a.long, source='geo', dest='mlt',
            datetime=a.index, height=a.height)
    mlat, mlongt = b.convert(
            lat=a.lat, lon=a.long, source='geo', dest='qd',
            height=a.height)
    mlat2 = np.array(a.Mlat)
    mlt2 = np.array(a.MLT)
    plt.plot(mlat-mlat2)
    plt.plot(abs(mlt-mlt2) % 24, '.')
    plt.show()
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def test_convert_invalid_lat():
    A = Apex(date=2000, refh=300)
    with pytest.raises(ValueError):
        A.convert(91, 0, 'geo', 'geo')
    with pytest.raises(ValueError):
        A.convert(-91, 0, 'geo', 'geo')
    A.convert(90, 0, 'geo', 'geo')
    A.convert(-90, 0, 'geo', 'geo')

    assert_allclose(A.convert(90+1e-5, 0, 'geo', 'apex'), A.convert(90, 0, 'geo', 'apex'), rtol=0, atol=1e-8)

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_temperature():
    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    g = g2a
    alts = [130, 400]
    for ialt, alt in enumerate(alts):
        alt_ind = np.argmin(np.abs(g['Altitude'][0, 0, 2:-2]/1000-alt))+2
        alt_str = '%6.2f' % (g['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
                1, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
                centrallon=g3ca.calculate_centrallon(g, 'polar',  useLT=True),
                coastlines=False)
        lon0, lat0, zdata0 = g3ca.contour_data('Temperature', g, alt=alt)
        fp = (lat0[:,0]>slat) & (lat0[:,0]<nlat)
        lon0, lat0, zdata0 = lon0[fp, :], lat0[fp,:], zdata0[fp,:]
        hc = ax.contourf(lon0, lat0, zdata0, 21,
                         transform=ccrs.PlateCarree(), cmap='jet',
                         extend='both')
        hc = plt.colorbar(hc, pad=0.17)
        hc.set_label('Temperature (K)')
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+titletime, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    plt.show()
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_ion_drift(show=True,save=True):
    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)
    g = g2a

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    for ialt, alt in enumerate([130, 200, 300, 400, 500, 600]):
        alt_ind = np.argmin(np.abs(g['Altitude'][0, 0, :]/1000-alt))
        alt_str = '%6.2f' % (g['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
                3, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
                centrallon=g3ca.calculate_centrallon(g, 'polar',  useLT=True),
                coastlines=False)
        lon0, lat0, ewind, nwind = g3ca.vector_data(g, 'ion', alt=alt)
        lon0, lat0, ewind, nwind = g3ca.convert_vector(
                lon0, lat0, ewind, nwind, plot_type='polar',
                projection=projection)
        hq = ax.quiver(
                lon0, lat0, ewind, nwind, scale=1500, scale_units='inches',
                regrid_shape=20, headwidth=5)
        ax.quiverkey(hq, 0.93, -0.1, 1000, '1000 m/s')
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'03_ion_drift_%s%s.pdf' % (tstrday,tstring))
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
    def satellite_position_lt_lat(self, mag=False, ns='N'):
        """ Show the lt and lat positions of the satellite in a polar
        coordinate.

        Input:
            mag: if True, for MLT and Mlat position
            ns: N or S for North and South hemispheres, respectively

        Output:
            hcup, hcdown: scatter handles for the up and down orbits,
                respectively.
        """
        if self.empty:
            return
        lt = self['LT']
        lat = self['lat']
        if mag:
            from apexpy import Apex
            gm = Apex()
            mlat,mlt = gm.convert(
                    self['lat'], self['long'], 'geo', 'mlt', datetime=self.index)
            lt, lat = mlt, mlat
        ct = lat>0 if ns is 'N' else lat<0
        theta = lt[ct]/12*np.pi
        r = 90 - abs(lat[ct])
        hc = plt.scatter(theta, r, linewidths=0)
        return hc

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_vert_vgradrho_rho_diff(show=True, save=True):
    rho1 = np.array(g1a['Rho'])
    vgradrho1 = \
        g1a['V!Dn!N (up)'] \
        * cr.calc_rusanov_alts_ausm(g1a['Altitude'],rho1)/g1a['Rho']

    rho2 = np.array(g2a['Rho'])
    vgradrho2 = \
        g2a['V!Dn!N (up)']\
        * cr.calc_rusanov_alts_ausm(g2a['Altitude'],rho2)/g2a['Rho']

    g1a['vgradrho_diff'] = vgradrho1-vgradrho2

    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    for ialt, alt in enumerate([130, 200, 300, 400, 500, 600]):
        alt_ind = np.argmin(np.abs(g1a['Altitude'][0, 0, :]/1000-alt))
        alt_str = '%6.2f' % (g1a['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
            3, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
            centrallon=g3ca.calculate_centrallon(g1a, 'polar',  useLT=True),
            coastlines=False)
        lon0, lat0, zdata0 = g3ca.contour_data('vgradrho_diff', g1a, alt=alt)
        fp = (lat0[:,0]>slat) & (lat0[:,0]<nlat)
        lon0, lat0, zdata0 = lon0[fp, :], lat0[fp,:], zdata0[fp,:]
        hc = ax.contourf(
            lon0, lat0, zdata0, np.linspace(-1,1,21)*1e-4,
            transform=ccrs.PlateCarree(), cmap='seismic', extend='both')
        hcb = plt.colorbar(hc, pad=0.17)
        hcb.formatter.set_powerlimits((0,0))
        hcb.update_ticks()
        hcb.set_label(r'$-\vec{u}\cdot\frac{\nabla\rho}{\rho}$ (up)')
        # wind difference
        lon1, lat1, ewind1, nwind1 = g3ca.vector_data(g1a, 'neutral', alt=alt)
        lon2, lat2, ewind2, nwind2 = g3ca.vector_data(g2a, 'neutral', alt=alt)
        lon0, lat0 = lon1, lat1
        lon0, lat0, ewind0, nwind0 = g3ca.convert_vector(
            lon0, lat0, ewind2-ewind1, nwind2-nwind1, plot_type='polar',
            projection=projection)
        hq = ax.quiver(
            lon0, lat0, ewind0, nwind0, scale=1000, scale_units='inches',
            regrid_shape=20, headwidth=5)
        ax.quiverkey(hq, 0.93, -0.1, 500, '500 m/s')
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'06_vert_vgradrho_rho_diff_%s%s.pdf' % (tstrday,tstring))
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_temperature_diff(show=True, save=True):
    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    for ialt, alt in enumerate([130, 200, 300, 400, 500, 600]):
        alt_ind = np.argmin(np.abs(g1a['Altitude'][0, 0, :]/1000-alt))
        alt_str = '%6.2f' % (g1a['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
            3, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
            centrallon=g3ca.calculate_centrallon(g1a, 'polar',  useLT=True),
            coastlines=False)
        # temperature diff
        lon1, lat1, zdata1 = g3ca.contour_data('Temperature', g1a, alt=alt)
        lon2, lat2, zdata2 = g3ca.contour_data('Temperature', g2a, alt=alt)
        fp = (lat1[:,0]>slat) & (lat1[:,0]<nlat)
        lon0, lat0, zdata0 = lon1[fp, :], lat1[fp,:], (zdata2-zdata1)[fp,:]
        hc = ax.contourf(
            lon0, lat0, zdata0, np.linspace(-80,80,21),
            transform=ccrs.PlateCarree(), cmap='seismic', extend='both')
        hc = plt.colorbar(hc, pad=0.17)
        hc.set_label(r'$T_2-T_1$ (K)')
        # geomagnetic pole
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())

        # wind difference
        lon1, lat1, ewind1, nwind1 = g3ca.vector_data(g1a, 'neutral', alt=alt)
        lon2, lat2, ewind2, nwind2 = g3ca.vector_data(g2a, 'neutral', alt=alt)
        lon0, lat0 = lon1, lat1
        lon0, lat0, ewind0, nwind0 = g3ca.convert_vector(
            lon0, lat0, ewind2-ewind1, nwind2-nwind1, plot_type='polar',
            projection=projection)
        hq = ax.quiver(
            lon0, lat0, ewind0, nwind0, scale=1000, scale_units='inches',
            regrid_shape=20, headwidth=5)
        ax.quiverkey(hq, 0.93, -0.1, 500, '500 m/s')

        # rho difference
        #lon3, lat3, zdata3 = g3ca.contour_data('Rho', g1a, alt=alt)
        #lon4, lat4, zdata4 = g3ca.contour_data('Rho', g2a, alt=alt)
        #diffrho = 100*(zdata4-zdata3)/zdata3
        #hc = ax.contour(
        #        lon4, lat4, diffrho, [-10], transform=ccrs.PlateCarree(),
        #        colors='g',linestyles='-')
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'02_temperature_diff_%s%s.pdf' % (tstrday,tstring))
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_den_win(show=True, save=True):
    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    g = g2a
    for ialt, alt in enumerate([130, 200, 300, 400, 500, 600]):
        alt_ind = np.argmin(np.abs(g['Altitude'][0, 0, 2:-2]/1000-alt))+2
        alt_str = '%6.2f' % (g['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
                3, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
                centrallon=g3ca.calculate_centrallon(g, 'polar',  useLT=True),
                coastlines=False)
        # density
        lon0, lat0, zdata0 = g3ca.contour_data('Rho', g, alt=alt)
        fp = (lat0[:,0]>slat) & (lat0[:,0]<nlat)
        lon0, lat0, zdata0 = lon0[fp, :], lat0[fp,:], zdata0[fp,:]
        hc = ax.contourf(lon0, lat0, zdata0, 21,
                         transform=ccrs.PlateCarree(), cmap='jet',
                         extend='both')
        hc = plt.colorbar(hc, pad=0.17)
        hc.set_label(r'$\rho$ (kg/m$^3$)')

        # wind
        lon0, lat0, ewind, nwind = g3ca.vector_data(g, 'neutral', alt=alt)
        lon0, lat0, ewind, nwind = g3ca.convert_vector(
                lon0, lat0, ewind, nwind, plot_type='polar',
                projection=projection)
        hq = ax.quiver(
                lon0, lat0, ewind, nwind, scale=1500, scale_units='inches',
                regrid_shape=20)
        ax.quiverkey(hq, 0.93, -0.1, 1000, '1000 m/s')
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())

        # rho difference
        # lon3, lat3, zdata3 = g3ca.contour_data('Rho', g1a, alt=alt)
        # lon4, lat4, zdata4 = g3ca.contour_data('Rho', g2a, alt=alt)
        # diffrho = 100*(zdata4-zdata3)/zdata3
        # hc = ax.contour(
        #         lon4, lat4, diffrho, [-10], transform=ccrs.PlateCarree(),
        #         colors='g',linestyles='-')
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'01_den_win_run2_%s%s.pdf' %(tstrday,tstring))
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
    def polar_quiver_wind(self, ax, ns='N'):
        # Wind vector in lat-long coordinates.
        # For different map projections, the arithmetics to calculate xywind
        # are different
        if self.empty:
            return
        from mpl_toolkits.basemap import Basemap
        from apexpy import Apex
        # Creat polar coordinates
        projection,fc = ('npstere',1) if ns=='N' else ('spstere',-1)
        m = Basemap(projection=projection,boundinglat=fc*40,lon_0=0,resolution='l')
        m.drawcoastlines(color='gray',zorder=1)
        m.fillcontinents(color='lightgray',zorder=0)
        dt = self.index.min() + (self.index.max()-self.index.min())/2
        m.nightshade(dt,zorder=2)
        #m.drawparallels(np.arange(-80,81,20))
        #m.drawmeridians(np.arange(-180,181,60),labels=[1,1,1,1])

        # Calculate mlat and mlon
        lat_grid = np.arange(-90,91,10)
        lon_grid = np.arange(-180,181,10)
        lon_grid, lat_grid = np.meshgrid(lon_grid, lat_grid)
        gm = Apex(date=2005)
        mlat,mlon = gm.convert(lat_grid,lon_grid,'geo','apex')
        hc1 = m.contour(lon_grid,lat_grid,mlat,levels=np.arange(-90,91,10),
                        colors='k', zorder=3, linestyles='dashed',
                        linewidths=1, latlon=True)
        # hc2 = m.contour(lon_grid,lat_grid,mlon,levels=np.arange(-180,181,45),
        #                 colors='k', zorder=3, linestyles='dashed', latlon=True)
        plt.clabel(hc1,inline=True,colors='k',fmt='%d')
        # plt.clabel(hc2,inline=True,colors='k',fmt='%d')

        # Calculate and plot x and y winds
        lat = self.lat
        lon = self.long
        wind = self.wind
        winde1 = self.winde
        winde = winde1*wind
        windn1 = self.windn
        windn = windn1*wind
        # only appropriate for the npstere and spstere
        xwind = fc*winde*np.cos(lon/180*np.pi)-windn*np.sin(lon/180*np.pi)
        ywind = winde*np.sin(lon/180*np.pi)+fc*windn*np.cos(lon/180*np.pi)
        hq = m.quiver(np.array(lon),np.array(lat),xwind,ywind,color='blue',
                      scale=300, scale_units='inches',zorder=3, latlon=True)
        #plt.quiverkey(hq,1.05,1.05,100,'100 m/s',coordinates='axes',labelpos='E')
        #m.scatter(np.array(lon),np.array(lat),
        #          s=50, c=self.index.to_julian_date(),linewidths=0, zorder=4,latlon=True)
        return m

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def mag_parallels(date, parallels=range(-75, 76, 15), height=350, N=1000):
    """
    Return a mapping between magnetic latitudes specified by
    *parallels* to the tuple of mapped geographic latitudes and
    longitudes. The mapping is made across *N* uniformly spaced
    geographic longitudes, on :class:`datetime` *date*, and at
    *height* (in [km]) in apex geomagnetic coordinates. If *date* is
    None, use the current date and time in the coordinate
    transformation.
    """
    apex = Apex(date=date)
    parallel_map = OrderedDict()
    lons = NP.linspace(-180, 180, N)
    for parallel in parallels:
        glat, glon = apex.convert(parallel, lons, source="apex", dest="geo")
        # sort by geographic longitude
        glat, glon = zip(*sorted(zip(glat, glon), key=lambda x: x[1]))
        parallel_map[parallel] = glat, glon
    return parallel_map

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
    def satellite_position_lt_lat(self, ax, nlat=90, slat=0,
                                  mag=False, plot_type='polar',
                                  *args, **kwargs):
        """ Show the (M)lt and (M)lat positions of the satellite in a polar
        or rectangular coordinate axis.

        Input:
            ax: which axis to draw orbit on
            nlat: northward latitude limit (default 90)
            slat: southward latitude limit (default 0)
            mag: if True, for MLT and Mlat position
        Return:
            hc: handle of the scatter
        """
        if self.empty:
            return
        tmp = self
        if mag:
            from apexpy import Apex
            import datetime as dt
            a = Apex(date=self.index.year.mean())
            mlat, mlt = a.convert(tmp.lat, tmp.long, 'geo','mlt',
                                  height=tmp.height, datetime=tmp.index)
            tmp['MLT'] = mlt
            tmp['Mlat'] = mlat
        ltp='MLT' if mag else 'LT'
        latp='Mlat' if mag else 'lat'
        ct = (self[latp]>slat) & (self[latp]<nlat)
        if 'pol' in plot_type.lower():
            if nlat*slat<0:
                print('Error: For polar plot, nlat and slat'
                      ' should have the same signs')
            csign = 1 if nlat>0 else -1
            theta = tmp.loc[ct,ltp]/12*np.pi
            r = 90-csign*tmp.loc[ct,latp]
            hc = ax.scatter(theta, r, linewidths=0, *args, **kwargs)
        if 'rec' in plot_type.lower():
            hc = ax.scatter(tmp.loc[ct, ltp], tmp.loc[ct, latp], linewidths=0,
                            *args, **kwargs)
        return hc

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_den_win():
    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    g = g2a
    alts = [130, 400]
    for ialt, alt in enumerate(alts):
        alt_ind = np.argmin(np.abs(g['Altitude'][0, 0, 2:-2]/1000-alt))+2
        alt_str = '%6.2f' % (g['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
                1, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
                centrallon=g3ca.calculate_centrallon(g, 'polar',  useLT=True),
                coastlines=False)
        # density
        lon0, lat0, zdata0 = g3ca.contour_data('Rho', g, alt=alt)
        fp = (lat0[:,0]>slat) & (lat0[:,0]<nlat)
        lon0, lat0, zdata0 = lon0[fp, :], lat0[fp,:], zdata0[fp,:]
        hc = ax.contourf(lon0, lat0, zdata0, 21,
                         transform=ccrs.PlateCarree(), cmap='jet',
                         extend='both')
        hc = plt.colorbar(hc, pad=0.17)
        hc.set_label(r'$\rho$ (kg/m$^3$)')
        # wind
        lon0, lat0, ewind, nwind = g3ca.vector_data(g, 'neutral', alt=alt)
        lon0, lat0, ewind, nwind = g3ca.convert_vector(
                lon0, lat0, ewind, nwind, plot_type='polar',
                projection=projection)
        hq = ax.quiver(
                lon0, lat0, ewind, nwind, scale=1500, scale_units='inches',
                regrid_shape=20)
        ax.quiverkey(hq, 0.93, -0.1, 1000, '1000 m/s')
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+titletime, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    plt.show()
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def test_convert_geo2mlt():
    datetime = dt.datetime(2000, 3, 9, 14, 25, 58)
    A = Apex(date=2000, refh=300)
    assert_allclose(A.convert(60, 15, 'geo', 'mlt', height=100, ssheight=2e5, datetime=datetime)[1], A.mlon2mlt(A.geo2apex(60, 15, 100)[1], datetime, ssheight=2e5))

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_divrhov_rho_diff(show=True, save=True):
    # density change (shrink)
    lon1 = np.array(g1a['Longitude'])
    lat1 = np.array(g1a['Latitude'])
    alt1 = np.array(g1a['Altitude'])
    Re = 6371*1000 # Earth radius, unit: m
    RR = Re+alt1
    omega = 2*np.pi/(24*3600)
    rho1 = np.array(g1a['Rho'])
    nwind1 = np.array(g1a['V!Dn!N (north)'])
    ewind1 = np.array(g1a['V!Dn!N (east)']) + omega*RR*np.cos(lat1)
    uwind1 = np.array(g1a['V!Dn!N (up)'])
    div_rhov1 = \
        (cr.calc_div_hozt(lon1, lat1, alt1, rho1*nwind1, rho1*ewind1)\
        +cr.calc_div_vert(alt1, rho1*uwind1))/rho1

    # density change (no shrink)
    lon2 = np.array(g2a['Longitude'])
    lat2 = np.array(g2a['Latitude'])
    alt2 = np.array(g2a['Altitude'])
    Re = 6371*1000 # Earth radius, unit: m
    RR = Re+alt2
    omega = 2*np.pi/(24*3600)
    rho2 = np.array(g2a['Rho'])
    nwind2 = np.array(g2a['V!Dn!N (north)'])
    ewind2 = np.array(g2a['V!Dn!N (east)']) + omega*RR*np.cos(lat2)
    uwind2 = np.array(g2a['V!Dn!N (up)'])
    div_rhov2 = \
        (cr.calc_div_hozt(lon2, lat2, alt2, rho2*nwind2, rho2*ewind2)\
        +cr.calc_div_vert(alt2, rho2*uwind2))/rho2
    g1a['divrhov_diff'] = div_rhov1-div_rhov2

    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)

    plt.close('all')
    plt.figure(figsize=(7.26, 9.25))
    for ialt, alt in enumerate([130, 200, 300, 400, 500, 600]):
        alt_ind = np.argmin(np.abs(g1a['Altitude'][0, 0, :]/1000-alt))
        alt_str = '%6.2f' % (g1a['Altitude'][0, 0, alt_ind]/1000)
        ax, projection = gcc.create_map(
            3, 2, ialt+1, 'polar', nlat=nlat, slat=slat, dlat=10,
            centrallon=g3ca.calculate_centrallon(g1a, 'polar',  useLT=True),
            coastlines=False)
        lon0, lat0, zdata0 = g3ca.contour_data('divrhov_diff', g1a, alt=alt)
        fp = (lat0[:,0]>slat) & (lat0[:,0]<nlat)
        lon0, lat0, zdata0 = lon0[fp, :], lat0[fp,:], zdata0[fp,:]
        hc = ax.contourf(
            lon0, lat0, zdata0, np.linspace(-1,1,21)*1e-4,
            transform=ccrs.PlateCarree(), cmap='seismic', extend='both')
        hc = plt.colorbar(hc, pad=0.17)
        hc.formatter.set_powerlimits((0,0))
        hc.update_ticks()
        hc.set_label(r'$-\frac{\nabla\cdot(\rho\vec{u})}{\rho}$')
        # wind difference
        lon1, lat1, ewind1, nwind1 = g3ca.vector_data(g1a, 'neutral', alt=alt)
        lon2, lat2, ewind2, nwind2 = g3ca.vector_data(g2a, 'neutral', alt=alt)
        lon0, lat0 = lon1, lat1
        lon0, lat0, ewind0, nwind0 = g3ca.convert_vector(
                lon0, lat0, ewind2-ewind1, nwind2-nwind1, plot_type='polar',
                projection=projection)
        hq = ax.quiver(
                lon0, lat0, ewind0, nwind0, scale=1000, scale_units='inches',
                regrid_shape=20, headwidth=5)
        ax.quiverkey(hq, 0.93, -0.1, 500, '500 m/s')
        ax.scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())
        plt.title('%s km' % alt_str, y=1.05)
    plt.text(0.5, 0.95, 'Time: '+tstring, fontsize=15,
            horizontalalignment='center', transform=plt.gcf().transFigure)
    if show:
        plt.show()
    if save:
        plt.savefig(spath+'08_divrhov_rho_diff_%s%s.pdf' % (tstrday,tstring))
    return

# 需要导入模块: from apexpy import Apex [as 别名]
# 或者: from apexpy.Apex import convert [as 别名]
def plot_time_constant(show=False, f107=150, which_alt=600):
    stime1 = pd.Timestamp('2003-03-22 00:10:00')
    etime1 = pd.Timestamp('2003-03-22 06:00:00')
    stime2 = pd.Timestamp('2003-03-23 00:10:00')
    etime2 = pd.Timestamp('2003-03-23 06:00:00')
    timeidx1 = pd.DatetimeIndex(start=stime1, end=etime1, freq='10min')
    timeidx2 = pd.DatetimeIndex(start=stime2, end=etime2, freq='10min')
    if f107==150:
        fp1 = '/media/guod/wd2t/simulation_output/momentum_analysis/'\
              + 'run_no_shrink_iondrift_4_1/data/'
        fn1 = [glob.glob(fp1+'3DALL_t'+k.strftime('%y%m%d_%H%M')+'*.bin')[0]
               for k in timeidx1]
        fp2 = '/media/guod/wd2t/simulation_output/momentum_analysis/'\
              + 'run_no_shrink_iondrift_4_all_day/data/'
        fn2 = [glob.glob(fp2+'3DALL_t'+k.strftime('%y%m%d_%H%M')+'*.bin')[0]
               for k in timeidx2]
        if which_alt==200:
            rholevels=np.linspace(1.8, 2.6, 21)*1e-10
        if which_alt==400:
            rholevels=np.linspace(2.7, 10, 21)*1e-12
        if which_alt==600:
            rholevels=np.linspace(0.1, 0.9, 21)*1e-12
    else:
        fp1 = '/home/guod/simulation_output/momentum_analysis/'\
              + 'run_shrink_70_continue/data/'
        fn1 = [glob.glob(fp1+'3DALL_t'+k.strftime('%y%m%d_%H%M')+'*.bin')[0]
               for k in timeidx]
        fp2 = '/home/guod/simulation_output/momentum_analysis/'\
              + 'run_no_shrink_70/data/'
        fn2 = [glob.glob(fp2+'3DALL_t'+k.strftime('%y%m%d_%H%M')+'*.bin')[0]
               for k in timeidx]
    apex = Apex(date=2003)
    qlat, qlon = apex.convert(-90, 0, source='apex', dest='geo', height=400)
    fig = plt.figure(figsize=[8,9])
    def animate_den_wind(i):
        g1, g2 = [gitm.GitmBin(k[i]) for k in [fn1, fn2]]
        # create axis
        ax = list(range(6))
        projection = ax.copy()
        for ins in range(2):
            nlat, slat = [90, 40] if ins==0 else [-40, -90]
            for irun in range(3):
                ax[ins+irun*2], projection[ins+irun*2] = gcc.create_map(
                        3, 2, 1+ins+irun*2, 'polar', nlat=nlat, slat=slat,
                        dlat=10, centrallon=g3ca.calculate_centrallon(
                            g1, 'polar',  useLT=True),
                        coastlines=False)
        # Density
        lon1, lat1, zdata1 = g3ca.contour_data('Rho', g1, alt=which_alt)
        lon2, lat2, zdata2 = g3ca.contour_data('Rho', g2, alt=which_alt)
        hc = [ax[k].contourf(
                lon1, lat1, zdata1, transform=ccrs.PlateCarree(),
                levels=rholevels,
                cmap='jet', extend='both') for k in [0, 1]]
        ax[0].set_title(g1['time'].strftime('%d-%b-%y %H:%M')+' UT', y=1.05)
        ax[1].set_title(g1['time'].strftime('%d-%b-%y %H:%M')+' UT', y=1.05)
        hc = [ax[k].contourf(
                lon2, lat2, zdata2, transform=ccrs.PlateCarree(),
                levels=rholevels,
                cmap='jet', extend='both') for k in [2, 3]]
        # diff density
        diffzdata = 100*(zdata2-zdata1)/zdata1
        hc = [ax[k].contourf(
                lon2, lat2, diffzdata, 21, transform=ccrs.PlateCarree(),
                levels=np.linspace(-30, 30, 21), cmap='seismic',
                extend='both') for k in [4, 5]]
        hc = [ax[k].contour(
                lon2, lat2, diffzdata, [-10], transform=ccrs.PlateCarree(),
                colors='g',linestyles='-') for k in [4, 5]]

        # wind
        lon1, lat1, ewind1, nwind1 = \
                g3ca.vector_data(g1, 'neutral', alt=which_alt)
        lon2, lat2, ewind2, nwind2 = \
                g3ca.vector_data(g2, 'neutral', alt=which_alt)
        for iax in range(6):
            if iax == 0 or iax == 1:
                lon0, lat0, ewind, nwind = (
                        lon1.copy(), lat1.copy(), ewind1.copy(), nwind1.copy())
                lon0, lat0, ewind, nwind = g3ca.convert_vector(
                        lon0, lat0, ewind, nwind, plot_type='polar',
                        projection=projection[iax])
            elif iax == 2 or iax == 3:
                lon0, lat0, ewind, nwind = (
                        lon2.copy(), lat2.copy(), ewind2.copy(), nwind2.copy())
                lon0, lat0, ewind, nwind = g3ca.convert_vector(
                        lon0, lat0, ewind, nwind, plot_type='polar',
                        projection=projection[iax])
            elif iax == 4 or iax == 5:
                lon0, lat0, ewind, nwind = (
                        lon1.copy(), lat1.copy(), ewind2-ewind1, nwind2-nwind1)
                lon0, lat0, ewind, nwind = g3ca.convert_vector(
                        lon0, lat0, ewind, nwind, plot_type='polar',
                        projection=projection[iax])
            hq = ax[iax].quiver(
                    lon0, lat0, ewind, nwind, scale=1500, scale_units='inches',
                    color='k', regrid_shape=20)
            ax[iax].scatter(qlon, qlat, color='k', transform=ccrs.PlateCarree())

            # ax.quiverkey(hq, 0.93, 0, 1000, '1000 m/s')
#.........这里部分代码省略.........