-3

I'm trying to do a blank Stüve diagram.
My code is:

R=287.04  #Jkg^-1K^-1
cp=1005   #Jkg^-1K^-1
p0=1000 #hPa
L=2.5*10**6 #J kg^-1
temp_celsius = np.array(range(-80,41))
temp_kelvin=temp_celsius+273.15
pressure_hPa=np.array(range(1050,90,-10))

#make a grid of the data
tempdata,pressdata=np.meshgrid(temp_kelvin,pressure_hPa*100.)

#Initialise the arrays
pot_temp_kelvin=np.zeros(tempdata.shape)
es=pot_temp_kelvin
ms=es
pseudo_pot_temp_kelvin=es

#Get the potential temperature                         
pot_temp_kelvin=tempdata*(p0*100/pressdata)**(R/cp)

#Get the saturation mix ratio
#first the saturation vapor pressure after Magnus
#Definition of constants for the Magnus-formula
c1=17.62
c2=243.12

es=6.112*np.exp((17.62*(tempdata-273.15))/(243.12+tempdata-273.15)) #hPa
#Now I'm calculation the saturation mixing ratio
ms=622*(es/(pressdata/100-es)) #g/kg

#At least I need the pseudo-adiabatic potential temperature
pseudo_pot_temp_kelvin=pot_temp_kelvin*np.exp(L*ms/1000./cp/tempdata)

#define the levels which should be plotted in the figure
levels_theta=np.array(range(200,405,5))
levels_ms=np.array([0.1,0.2,0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 8.0, 10.0, 12.0, 15.0, 20.0, 25.0, 30.0])
levels_theta_e=np.array(range(220,410,10))

#The plot
fig = plt.figure(figsize=(15,15))
theta=plt.contour(temp_celsius,pressure_hPa,pot_temp_kelvin,levels_theta,colors='blue')
plt.clabel(theta,levels_theta[0::2], fontsize=10, fmt='%1.0f')

sat_mix_ratio=plt.contour(temp_celsius,pressure_hPa,ms,levels_ms,colors='green')
plt.clabel(sat_mix_ratio,fontsize=10,fmt='%1.1f')
    theta_e=plt.contour(temp_celsius,pressure_hPa,pseudo_pot_temp_kelvin,levels_theta_e,colors='red')
plt.clabel(theta_e,levels_theta_e[1::2],fontsize=10,fmt='%1.0f')

plt.xlabel('Temperature [$^\circ$C]')
plt.ylabel('Pressure [hPa]')
plt.xticks(range(-80,45,5))
plt.xlim((-80,50))
plt.yticks(range(1050,50,-50))
plt.gca().invert_yaxis()
plt.grid(color='black',linestyle='-')
plt.show()

Everything works well but the real Stüve diagram should look like Stüve diagram with sounding

As you can see, the y-axis has a specific scale: p**(R/cp).... (=p**(287.04)/1005)

What do I have to do with my program so that my y-axis looks like the axis in the example?

Schorsch
  • 7,761
  • 6
  • 39
  • 65
  • Can you provide a code? – Farseer Dec 29 '14 at 12:35
  • possible duplicate of [Multiple y-axis conversion scales](http://stackoverflow.com/questions/25159495/multiple-y-axis-conversion-scales) - if this does not help you, please consider clarifying your question. – Schorsch Dec 29 '14 at 13:48

1 Answers1

0

You will have to define your own axis scale for this. My answer is based on this answer and the custom scale example.

This is the custom scaling class:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FormatStrFormatter
from matplotlib import scale as mscale
from matplotlib import transforms as mtransforms

class CustomScale(mscale.ScaleBase):
    name = 'custom'

    def __init__(self, axis, **kwargs):
        mscale.ScaleBase.__init__(self)
        self.thresh = None #thresh

    def get_transform(self):
        return self.CustomTransform(self.thresh)

    def set_default_locators_and_formatters(self, axis):
        pass

    class CustomTransform(mtransforms.Transform):
        input_dims = 1
        output_dims = 1
        is_separable = True

        def __init__(self, thresh):
            mtransforms.Transform.__init__(self)
            self.thresh = thresh

        def transform_non_affine(self, a):
            return a**(R/cp)

        def inverted(self):
            return CustomScale.InvertedCustomTransform(self.thresh)

    class InvertedCustomTransform(mtransforms.Transform):
        input_dims = 1
        output_dims = 1
        is_separable = True

        def __init__(self, thresh):
            mtransforms.Transform.__init__(self)
            self.thresh = thresh

        def transform_non_affine(self, a):
            return a**(cp/R)

        def inverted(self):
            return CustomScale.CustomTransform(self.thresh)

# Now that the Scale class has been defined, it must be registered so
# that ``matplotlib`` can find it.
mscale.register_scale(CustomScale)

Now, you can use the following option to provide your custom scale to your y=axis:

plt.gca().set_yscale('custom')

The following plot compares the custom scale to a log-scale (plt.gca().set_yscale('log')) and the default without scale:

example

Community
  • 1
  • 1
Schorsch
  • 7,761
  • 6
  • 39
  • 65