Add Muneer transposition model

[BernatNicolau]

• [X] Closes to #2117 (this comment)
• [X] I am familiar with the contributing guidelines
• [x] Tests added
• [x] Updates entries in docs/sphinx/source/reference for API changes.
• [x] Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• [x] New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• [x] Pull request is nearly complete and ready for detailed review.
• [x] Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

[BernatNicolau]

It is quite hard to do the testing of the code since PVGIS directly provides the transposed data, without providing GHI, DHI and DNI (if I am wrong please tell me). I wanted to cross check if I obtained the same results by downloading the transposed signal and by applying it with the function, but I think I am reaching a dead end here...

[echedey-ls]

I think I am reaching a dead end here

Not really! Most of the time we don't have original data to work with. It's a good enough test if you can make up some data by using a spreadsheet like Excel and implementing the model there. Just rewriting the equations and making input data up for DHI and b values, then typing both the inputs and expected output in the test should be fine. That way we can confirm this new code complies with the expected behaviour.

[kandersolar]

since PVGIS directly provides the transposed data, without providing GHI, DHI and DNI (if I am wrong please tell me)

It may be possible to recover the original GHI/DHI/DNI components by specifying surface_tilt=0. Of course the poa_direct would need to be corrected for zenith angle to get back to DNI.

[BernatNicolau]

I have come to realize that the Muneer transposition is not as straightforward as I thought. At the beginning I was following this paper which I understand is the original Muneer transposition. However, I have found this article that describes a Muneer model, referencing this book from 2004 which is more complex, but can surely be implemented.

To be able to code the 2004 Muneer model, a inverse_square_spencer71() function needs to be added in pvlib.solarposition, providing the inverse square described in the linked Fourier paper:

$$\frac{1}{r^2} = 1.000110 + 0.034221 \cos(T) + 0.001280 \sin(T) + 0.000719 \cos(2T) + 0.000077 \sin(2T)$$

I am not sure which model should pvlib offer (1990 or 2004). I would gladly accept some guidance here :)

[adriesse]

I am not sure which model should pvlib offer (1990 or 2004). I would gladly accept some guidance here :)

I suggest implementing the version from the book.

PS. Somewhere I have a floppy disk containing the software from the book, in case it is needed. PPS. And somewhere, I also have a floppy disk drive!

[BernatNicolau]

since PVGIS directly provides the transposed data, without providing GHI, DHI and DNI (if I am wrong please tell me)

It may be possible to recover the original GHI/DHI/DNI components by specifying surface_tilt=0. Of course the poa_direct would need to be corrected for zenith angle to get back to DNI.

This could be a good standalone feature to be included!

[cwhanse]

@BernatNicolau I would choose one (my preference is also 2004) and name it muneer2004 or something similar.

[BernatNicolau]

@BernatNicolau I would choose one (my preference is also 2004) and name it muneer2004 or something similar.

After having read the book section, I finally understood that the model is the same as 1990 but explained with greater detail and with some examples that compare different transposition models. So as far as I understand, there is no need to name it muneer2004 and can be kept muneer.

Note that I still need to do testing and add better documentation.

[BernatNicolau]

As a sanity check, I have downloaded data of PVGIS (POA) for a certain location, and have compared it with the Muneer transposition using Solcast data (GHI, DHI and DNI) in 15min, for the same location and resampled in a monthly basis, considering the same structure parameters (tilt=30):

Results are quite aligned (slope of 1.0035), although there is some scatter. (Note that the solar resource is different in each dataset)

I did the same exercise but using the Perez perez transposition to the Solcast data, obtaining bigger differences (slope of 0.9844):

[adriesse]

(Note that the solar resource is different in each dataset)

It's very odd that the scatter patterns are extremely similar--maybe that's a clue to something the comparisons have in common?

[BernatNicolau]

(Note that the solar resource is different in each dataset)

It's very odd that the scatter patterns are extremely similar--maybe that's a clue to something the comparisons have in common?

What I ment by that sentence is that the scatter in the first figure could be due to the intrinsic differences between PVGIS and Solcast, and could be unrelated to the transposition model.

The patterns comparing Muneer and Perez are quite similar, because the correlation between both models is the following (in a monthly basis):

[kandersolar]

(Attempting to) reproduce PVGIS's transposition can be done by retrieving GHI,DHI via the TMY endpoint.

Click to show code

import pvlib
import pandas as pd

surface_tilt = 20
surface_azimuth = 180
latitude = 40
longitude = -80

# get PVGIS GHI, DHI for one month by requesting the TMY dataset:
pvgis_tmy, months, _, meta = pvlib.iotools.get_pvgis_tmy(latitude, longitude)
pvgis_tmy = pvgis_tmy.loc[pvgis_tmy.index.month == 1]

# now get PVGIS POA for the same year/month:
year = months[0]['year']  # the year for the TMY january
pvgis_poa, _, meta = pvlib.iotools.get_pvgis_hourly(latitude, longitude, start=year, end=year, surface_tilt=surface_tilt, surface_azimuth=surface_azimuth, usehorizon=False)
pvgis_poa = pvgis_poa.loc[pvgis_poa.index.month == 1]  # keep only january

# transpose GHI, DHI to POA:
sp = pvlib.solarposition.get_solarposition(pvgis_poa.index + pd.Timedelta('30T'), latitude, longitude)
sp.index = pvgis_poa.index
dni_extra = pvlib.irradiance.get_extra_radiation(pvgis_poa.index)
poa_diffuse_muneer = pvlib.irradiance.muneer(surface_tilt, surface_azimuth, pvgis_tmy['dhi'], pvgis_tmy['ghi'], dni_extra, solar_azimuth=sp['azimuth'], solar_zenith=sp['zenith'])

# compare:
out = pd.DataFrame({
    'PVGIS': pvgis_poa['poa_sky_diffuse'],
    'muneer': poa_diffuse_muneer,
})
out.plot.scatter('PVGIS', 'muneer')

Do we know how PVGIS calculates the b parameter?

[BernatNicolau]

(Attempting to) reproduce PVGIS's transposition can be done by retrieving GHI,DHI via the TMY endpoint.

Click to show code

import pvlib
import pandas as pd

surface_tilt = 20
surface_azimuth = 180
latitude = 40
longitude = -80

# get PVGIS GHI, DHI for one month by requesting the TMY dataset:
pvgis_tmy, months, _, meta = pvlib.iotools.get_pvgis_tmy(latitude, longitude)
pvgis_tmy = pvgis_tmy.loc[pvgis_tmy.index.month == 1]

# now get PVGIS POA for the same year/month:
year = months[0]['year']  # the year for the TMY january
pvgis_poa, _, meta = pvlib.iotools.get_pvgis_hourly(latitude, longitude, start=year, end=year, surface_tilt=surface_tilt, surface_azimuth=surface_azimuth, usehorizon=False)
pvgis_poa = pvgis_poa.loc[pvgis_poa.index.month == 1]  # keep only january

# transpose GHI, DHI to POA:
sp = pvlib.solarposition.get_solarposition(pvgis_poa.index + pd.Timedelta('30T'), latitude, longitude)
sp.index = pvgis_poa.index
dni_extra = pvlib.irradiance.get_extra_radiation(pvgis_poa.index)
poa_diffuse_muneer = pvlib.irradiance.muneer(surface_tilt, surface_azimuth, pvgis_tmy['dhi'], pvgis_tmy['ghi'], dni_extra, solar_azimuth=sp['azimuth'], solar_zenith=sp['zenith'])

# compare:
out = pd.DataFrame({
    'PVGIS': pvgis_poa['poa_sky_diffuse'],
    'muneer': poa_diffuse_muneer,
})
out.plot.scatter('PVGIS', 'muneer')

Do we know how PVGIS calculates the b parameter?

Thanks Kevin for this test, this is a clever way of obtaining poa, ghi and dhi for the same period!

I noted that there is a shift between solar_elevation from PVGIS and the sp['elevation'] of 30min, which if we delete the + pd.Timedelta('30T') when calculating sp, it disappears.

This solves the scatter issue from your figure:

However, I am still trying to figure out which b value is used in PVGIS, which seems to be playing a role in the slope of the correlation.

When comparing overall poa values, the slope between PVGIS (x) and muneer (y), are:

• b=0, y=1.0012x
• b=5.73, y=0.9978x

I think these are pretty good results!

[adriesse]

I think these are pretty good results!

Me too. Plotting individual values rather than monthly averages/sums is also better. You could plot the difference or ratio now to gain insight into the remaining deviations.

[adriesse]

However, I am still trying to figure out which b value is used in PVGIS, which seems to be playing a role in the slope of the correlation.

I think you've implemented something between the Steven-Unsworth and Muneer models. b is not a parameter for the Muneer model.

[BernatNicolau]

I think these are pretty good results!

Me too. Plotting individual values rather than monthly averages/sums is also better. You could plot the difference or ratio now to gain insight into the remaining deviations.

The plots from @kandersolar and from my last comment are hourly values, unlike my first comment comparing PVGIS with Solcast, which were monthly aggregates.

[adriesse]

I think these are pretty good results!

Me too. Plotting individual values rather than monthly averages/sums is also better. You could plot the difference or ratio now to gain insight into the remaining deviations.

The plots from @kandersolar and from my last comment are hourly values, unlike my first comment comparing PVGIS with Solcast, which were monthly aggregates.

Yes, I noticed the evolution in the analysis and was trying to express a thumbs-up.

[BernatNicolau]

You could plot the difference or ratio now to gain insight into the remaining deviations.

I have plotted the difference between PVGIS sky diffuse and the modelled (with b=0), and we can see that most of the differences ocurr during morning/afternoon hours for specific days.

Note that I have only plotted the periods where PVGIS sky diffuse > 0. Also note that the magnitude of the differences is not so high.

I have no idea how to proceed with the validation of the model. It seems that the model is not clearly wrong but we are not obtaining the exact value compared with PVGIS. Do you have any suggestion on how can we validate this model further?

[kandersolar]

It's getting pretty close! I suspect at least some explanation for the remaining differences can be found in the Hofierka document mentioned in #2117. For example I bet the big spikes in that difference plot are related to the change in behavior at solar elevation of 5.7 degrees. The extraterrestrial radiation may also be calculated slightly differently.

[RDaxini]

@echedey-ls good point, the units question is probably out of scope for this PR --- my bad. It'd be good if you started a discussion on the subject.

[BernatNicolau]

I bring good news, I have created the scenario of sky diffuse when the solar altitude is low (<0.1rad) as per suggestion of @kandersolar in this comment and in https://github.com/pvlib/pvlib-python/issues/2117.

This has solved the negative spikes:

I have tried to solve the positive spikes, but I think this is an error of PVGIS. Take a look of this day, where the last hour of the day presents the spike (2014-01-23 22:00:00+00:00):

	PVGIS dhi	PVGIS ghi	PVGIS sky_diffuse	pvlib_muneer sky_diffuse	delta
2014-01-23 13:00:00+00:00	8.00	8.00	7.76	7.76	0.00
2014-01-23 14:00:00+00:00	33.00	33.00	32.00	32.00	0.00
2014-01-23 15:00:00+00:00	77.00	77.00	74.66	74.68	-0.02
2014-01-23 16:00:00+00:00	202.00	264.00	208.53	207.36	1.17
2014-01-23 17:00:00+00:00	91.00	91.00	88.24	88.26	-0.02
2014-01-23 18:00:00+00:00	101.00	101.00	97.93	97.95	-0.02
2014-01-23 19:00:00+00:00	138.00	138.00	133.81	133.84	-0.03
2014-01-23 20:00:00+00:00	155.00	199.00	159.98	159.10	0.88
2014-01-23 21:00:00+00:00	66.00	222.00	94.16	92.03	2.13
2014-01-23 22:00:00+00:00	43.00	43.00	50.88	41.70	9.18

It is a timestamp with 100% diffuse, and the sky diffuse obtained from PVGIS is greater than the dhi. The one obtained from pvlib seems correct to me.

I still haven't figured out which b value is used in PVGIS. I am obtaining best results with b=0.

[kandersolar]

After looking into this some more, it seems to me that we have to regard the PVGIS/Hofierka model as being related to but distinct from the Muneer 1990 model. Compared with Muneer, Hofierka uses different equations for different conditions (aoi < 0, solar_elevation < 5.7) and does away entirely with the b parameter. I don't see how we could have one function that cleanly implements both models.

So, we need to decide which of the two models we are aiming to implement here. If we want to implement Muneer, we cannot use PVGIS to produce validation data for the tests. @AdamRJensen, this was your feature request -- any thoughts on how we should proceed?

[BernatNicolau]

I agree with you @kandersolar . However, either way I don't think we can use PVGIS for validation, I was recently doing more testing and I am obtaining very spurious results while analyzing poa from PVGIS.

I modified your testing code to produce the 12 months of the TMY for that example location:

Click to show code

import pvlib
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import pvlib.tools as tools

surface_tilt = 20
surface_azimuth = 180
latitude = 40
longitude = -80

# get PVGIS GHI, DHI for one month by requesting the TMY dataset:
pvgis_tmy, months, _, meta = pvlib.iotools.get_pvgis_tmy(latitude, longitude)
pvgis_all = pd.DataFrame()

for y in range(1, 13):
  print(y)
  pvgis_tmy_year = pvgis_tmy.loc[pvgis_tmy.index.month == y]

  # now get PVGIS POA for the same year/month:
  year = months[y-1]['year']  # the year for the TMY january
  pvgis_poa, _, meta = pvlib.iotools.get_pvgis_hourly(
      latitude, longitude, start=year, end=year, surface_tilt=surface_tilt, surface_azimuth=surface_azimuth, usehorizon=False)
  pvgis_poa = pvgis_poa.loc[pvgis_poa.index.month == y]  # keep only january

  # transpose GHI, DHI to POA:
  sp = pvlib.solarposition.get_solarposition(
      pvgis_poa.index, latitude, longitude)

  sp.index = pvgis_poa.index
  dni_extra = pvlib.irradiance.get_extra_radiation(pvgis_poa.index)

  poa_diffuse_muneer = pvlib.irradiance.muneer(
      surface_tilt, surface_azimuth, pvgis_tmy_year['dhi'], pvgis_tmy_year['ghi'], dni_extra, solar_azimuth=sp['azimuth'], solar_zenith=sp['zenith'], b=0)
  poa_muneer = pvlib.irradiance.get_total_irradiance(
      surface_tilt=surface_tilt,
      surface_azimuth=surface_azimuth,
      solar_zenith=sp['zenith'],
      solar_azimuth=sp['azimuth'],
      dni=pvgis_tmy_year['dni'],
      ghi=pvgis_tmy_year['ghi'],
      dhi=pvgis_tmy_year['dhi'],
      dni_extra=dni_extra,
      model="muneer",
  )

  aux = pd.DataFrame({
      'sky_diff-PVGIS': pvgis_poa['poa_sky_diffuse'],
      'sky_diff-muneer': poa_diffuse_muneer,
      'poa-PVGIS': pvgis_poa.iloc[:, 0:2].sum(axis=1),
      'poa-muneer': poa_muneer['poa_global']
  })
  out = pd.concat([aux, pvgis_tmy_year, sp, dni_extra], axis=1)
  pvgis_all = pd.concat([pvgis_all, out])



pvgis_all.to_csv('comparison_allyear.csv')

There are some months of the TMY where the poa value seems correct, but there are other months that the poa doesn't make any sense to me. Let me illustrate:

On this day, the poa has a lower value than the ghi, and it seems to be a non-cloudy day.

This day doesn't make sense to me either...

[adriesse]

I thought I had commented on the b parameter earlier, but maybe my comment got lost. In this FORTRAN code that accompanied the 2004 edition of his book, the b parameter is already gone. I suggest implementing this since it is well-defined and documented, but of course if PVGIS does use something different, then that should be implemented as well.

C     PROG4-2.FOR

C     PROGRAM FOR SLOPE IRRADIANCE - OUTPUT FOR SEVEN MODELS 
	Character*1 Anum
    1 format(A1)


      DOUBLE PRECISION EPSBINS(7)
      DOUBLE PRECISION F11R(8), F12R(8), F13R(8)
      DOUBLE PRECISION F21R(8), F22R(8), F23R(8)
      DOUBLE PRECISION B2  
      DOUBLE PRECISION Z,ZH
      
      PI=3.14159
      DTOR=3.14159/180.0
      
      WRITE(*,*)'SLOPE IRRADIANCE - OUTPUT FOR SEVEN MODELS'
      WRITE(*,*)'    '
      WRITE(*,*) 'SELECT THE SYSTEM USED FOR TIME LOG'
      WRITE(*,*) 'INPUT "1" FOR SOLAR, OR "2" FOR LOCAL CLOCK TIME'
      READ(*,*) NTIMES
      
      IF(NTIMES.EQ.2) THEN
      WRITE(*,*) 'INPUT LAT, LONG, STD. TIME MERIDIAN (real values)'
      WRITE(*,*) 'NORTH = +,  WEST = + '
      READ(*,*) YLAT,YLONG,YRLONG
      ELSE                        
      WRITE(*,*) 'INPUT LATITUDE (real value)'
      WRITE(*,*) 'NORTH = +'
      READ(*,*) YLAT
      YRLONG=0.0
      YLONG=0.0
      ENDIF
     
      WRITE(*,*) 'INPUT SURFACE AZIMUTH AND TILT (real values)'
      WRITE(*,*)'AZIMUTH - NORTH=0, EAST=90, SOUTH=180, WEST=270 DEG'
      WRITE(*,*) 'SURFACE TILT - VERTICAL = 90 DEG'
      READ(*,*) AZI,TLT    

      
      WRITE(*,*)'INPUT YEAR (YYYY),MONTH,DAY,HOUR(0-23),MINUTE'
      write(*,*)'only integers to be used'    
      READ(*,*) IYR,IMT,IDY,IHR,IME 

      WRITE(*,*)'ALL IRRADIANCE VALUES ARE TO BE EXPRESSED IN W/m2'
      WRITE(*,*)'NOTE: 1 Wh/m2 = 0.0036 MJ/m2'
      WRITE(*,*)'INPUT HORIZONTAL GLOBAL,DIFFUSE IRRADIANCE , ALBEDO'
      READ(*,*) GRAD,DRAD,RHO
            
C     Calculate GHA360, DEC 
      XLCT=(1.0*IHR)+(1.0*IME/60.0)
      UT=XLCT+(YRLONG/15.0)      
      IF (IMT.GT.2) THEN
      IYR1=IYR
      IMT1=IMT-3
      ELSE
      IYR1=IYR-1
      IMT1=IMT+9
      ENDIF
      INTT1=INT(30.6*IMT1+0.5)
      INTT2=INT(365.25*(IYR1-1976))
      SMLT=((UT/24.0)+IDY+INTT1+INTT2-8707.5)/36525.0
      EPSILN=23.4393-0.013*SMLT
      CAPG=357.528+35999.050*SMLT
      IF(CAPG.GT.360.0) THEN
      G360=CAPG-INT(CAPG/360.0)*360.0
      ELSE
      G360=CAPG
      ENDIF
      CAPC=1.915*SIN(G360*DTOR)+0.020*SIN(2*G360*DTOR)
      CAPL=280.460+36000.770*SMLT+CAPC
      IF(CAPL.GT.360.0) THEN
      XL360=CAPL-INT(CAPL/360.0)*360.0
      ELSE
      XL360=CAPL
      ENDIF
      ALPHA=XL360-2.466*SIN(2*XL360*DTOR)+0.053*SIN(4*XL360*DTOR)
      GHA=15.0*UT-180.0-CAPC+XL360-ALPHA
      IF(GHA.GT.360.0) THEN
      GHA360=GHA-INT(GHA/360.0)*360.0
      ELSE
      GHA360=GHA
      ENDIF
      DEC=ATAN(TAN(EPSILN*DTOR)*SIN(ALPHA*DTOR))/DTOR

C     Calculate Solar Hour Angle 
      IF(NTIMES.EQ.2) THEN
      SHA=GHA360-(YLONG)
      ELSE      
      SHA=GHA360
      ENDIF 
      
C     Calculate Apparent Solar Time
      IF(NTIMES.EQ.2) THEN
      AST=12.0+(SHA/15.0)
      ELSE
      AST=XLCT
      ENDIF

C     Calculate SOLar ALTitude      
      TRM111=SIN(YLAT*DTOR)*SIN(DEC*DTOR)
      TRM112=COS(YLAT*DTOR)*COS(DEC*DTOR)*COS((SHA+180.0)*DTOR)
      TRM11=TRM111-TRM112
      SOLALT=ASIN(TRM11)/DTOR  
      
C     Calculate SOLar AZiMuth 
      TRM121=COS(YLAT*DTOR)*TAN(DEC*DTOR)
      TRM122=SIN(YLAT*DTOR)*COS((SHA+180.0)*DTOR)
      TRM12=TRM121+TRM122
      TRM14=TRM12*COS(DEC*DTOR)/COS(SOLALT*DTOR)
      SOLAZM=ACOS(TRM14)/DTOR
      IF(SHA.GT.0.0) THEN
      SOLAZM=360.0-SOLAZM
      ENDIF

C     Calculate SOLar INCidence angle
      TRM15=COS(TLT*DTOR)*SIN(SOLALT*DTOR)+SIN(TLT*DTOR)*COS(SOLALT
     1*DTOR)*COS((SOLAZM-AZI)*DTOR)
      SOLINC=ACOS(TRM15)/DTOR  
      
C     Calculate Day Number
      DN1=(IDY+INTT1+INTT2)
      IMTX=1
      IYR1=IYR-1
      IMT1=IMTX+9
      INTT1=INT(30.6*IMT1+0.5)
      INTT2=INT(365.25*(IYR1-1976))
      DN2=(INTT1+INTT2)
      DN=DN1-DN2
      
C     Calculate horizontal Extraterrestrial irRADiance & kt
      ERAD=1353.0*(1.+0.033*COS(0.0172024*DN))*SIN(SOLALT*DTOR)  
      TKT=GRAD/ERAD

C     BEAM IRRADIATION      
      IF(SOLINC.LT.90.0) THEN
      BSRAD=(GRAD-DRAD)*COS(SOLINC*DTOR)/SIN(SOLALT*DTOR)
      ELSE
      BSRAD=0.0
      ENDIF

C     GROUND REFLECTED RADIATION
      RSRAD=RHO*GRAD*SIN(TLT*DTOR/2.0)**2

C     SKY IRRADIANCE - ISOTROPIC MODEL
      TLTISO=(COS(TLT*DTOR/2.)**2)
      SISO=TLTISO*DRAD
                               
C     SKY IRRADIANCE - HAY'S MODEL  
      CLRFRA=(GRAD-DRAD)/ERAD
      TLTFAC=(1.0-CLRFRA)*(COS(TLT*DTOR/2.)**2) 
      
      IF(SOLINC.LT.90.0) THEN
      TLTFAC=(CLRFRA*TRM15/TRM11)+TLTFAC
      ENDIF 
      
      SHAY=TLTFAC*DRAD   
C     SKY IRRADIANCE - REINDL'S MODEL
      CLRFRA=(GRAD-DRAD)/ERAD 
      REINDF=SQRT((GRAD-DRAD)/GRAD)
      TLTFAC=(1.0-CLRFRA)*(COS(TLT*DTOR/2.)**2)*
     1(1.+REINDF*SIN(TLT*DTOR/2)**3)
      
      IF(SOLINC.LT.90.0) THEN     
      TLTFAC=(CLRFRA*TRM15/TRM11)+TLTFAC
      ENDIF 
      
      SREIND=TLTFAC*DRAD 
      
C     SKY IRRADIANCE - SKARTVEIT-OLSETH'S MODEL  
      CLRFRA=(GRAD-DRAD)/ERAD  
      BKB=MIN(CLRFRA,1.0)
      BRB=MAX((TRM15/TRM11),0.0)
      CAPBSO=MAX((0.3-2.*BKB),0.0)  
      SSKOL=DRAD*(BKB*BRB+(CAPBSO*COS(TLT*DTOR))+(1.-BKB-CAPBSO)*TLTISO)
      
C     SKY IRRADIANCE - GUEYMARD'S MODEL  
      H01=SOLALT*0.01
      GA0=-0.897-3.364*H01+(3.96*H01**2)-1.909*H01**3
      GA1=4.448-12.962*H01+(34.601*H01**2)-48.784*H01**3+(27.511*H01**4)
      GA2=-2.77+9.164*H01-(18.876*H01**2)+23.776*H01**3-(13.014*H01**4)
      GA3=0.312-0.217*H01-(0.805*H01**2)+0.318*H01**3
      
      CFS=(1.0-(0.2249*SIN(TLT*DTOR)**2)+(0.1231*SIN(2.0*TLT*DTOR))
     1-0.0342*SIN(4.0*TLT*DTOR))/(1.0-0.2249)
      CGH=0.408-0.323*H01+(0.384*H01**2)-0.17*H01**3
      CRD0=EXP(GA0+GA1*TRM15+(GA2*TRM15**2)+GA3*TRM15**3)
      CRD0=CRD0+CFS*CGH  
                                 
      IF((DRAD/GRAD).LE.0.227) THEN 
      CAPY=6.6667*(DRAD/GRAD)-1.4167
      ELSE
      CAPY=1.2121*(DRAD/GRAD)-0.1758
      ENDIF
      CNPT=MAX(MIN(CAPY,1.0),0.0)
      SMLB=0.5+CNPT                          
      CAPBG=2.*SMLB/(PI*(3.+2.*SMLB))
      CRD1=(COS(TLT*DTOR/2.)**2)+CAPBG*(SIN(TLT*DTOR)-(TLT*DTOR)
     1*COS(TLT*DTOR)-PI*SIN(TLT*DTOR/2.)**2)
      SGUEYM=DRAD*((1.0-CNPT)*CRD0+CNPT*CRD1)
                                                        
C     SKY IRRADIANCE - MUNEER'S MODEL
      CLRFRA=(GRAD-DRAD)/ERAD 
      
      IF(SOLINC.GE.90.0) THEN
      CAPB=0.252
      CLRFRA=0.0
      ELSE
c     The user may select one of the following four models:
    
c     CAPB= 0.003 33 -0.415 F -0.698 7 F**2 [for Northern Europe]
c     CAPB= 0.002 63 -0.712 F -0.688 3 F**2 [for Southern Europe]
c     CAPB= 0.080 00 -1.050 F -2.840 0 F**2 [for Japan]
c     CAPB= 0.040 00 -0.820 F -2.026 0 F**2 [for the globe]

c     Model for Northern Europe       
      CAPB=0.00333-0.415*CLRFRA-0.6987*CLRFRA**2 
c     End of model for Northern Europe      
      ENDIF 
      
      TLTFAC=(COS(TLT*DTOR/2.)**2)+CAPB*(SIN(TLT*DTOR)-(TLT*DTOR)
     1*COS(TLT*DTOR)-PI*SIN(TLT*DTOR/2.)**2)    
      DSRAD=TLTFAC*(1.-CLRFRA)+CLRFRA*COS(SOLINC*DTOR)/SIN(SOLALT*DTOR)
      SMUN=DSRAD*DRAD     
      
C     SKY IRRADIANCE - PEREZ MODEL
      DATA F11R / -0.0083117, 0.1299457, 0.3296958, 0.5682053,
     10.8730280, 1.1326077, 1.0601591, 0.6777470 /
      DATA F12R / 0.5877285, 0.6825954, 0.4868735, 0.1874525,
     1-0.3920403, -1.2367284, -1.5999137, -0.3272588 /
      DATA F13R / -0.0620636, -0.1513752, -0.2210958, -0.2951290,
     1-0.3616149, -0.4118494, -0.3589221, -0.2504286 /
      DATA F21R / -0.0596012, -0.0189325, 0.0554140, 0.1088631,
     10.2255647, 0.2877813, 0.2642124, 0.1561313 /
      DATA F22R / 0.0721249, 0.0659650, -0.0639588, -0.1519229,
     1-0.4620442, -0.8230357, -1.1272340, -1.3765031 /
      DATA F23R / -0.0220216, -0.0288748, -0.0260542, -0.0139754,
     10.0012448, 0.0558651, 0.1310694, 0.2506212 /
      DATA EPSBINS / 1.065, 1.23, 1.5, 1.95, 2.8, 4.5, 6.2 /
      B2 = 5.534D-6    
      G=GRAD
      B=(GRAD-DRAD)/TRM11
      Z=(90.0-SOLALT)*DTOR
      SLOPE=TLT*DTOR
      XPINC=SOLINC*DTOR      
      ZENITH =Z/DTOR
      ZH = MAX(TRM11, 0.0871557)
      AIRMASS = 1.0 / (TRM11 + 0.15 * (93.9 - ZENITH)** (-1.253))
      DELTA = DRAD * AIRMASS / 1367.0
      T = ZENITH**3.0
      EPS = (B + DRAD) / DRAD
      EPS = (EPS + T * B2) / (1.0 + T * B2)  

      DO 70001 I = 1,8
      IF ( (I .EQ. 8) .OR. (EPS .LE. EPSBINS(I)) ) GOTO 70002
70001 CONTINUE
70002 F1 = MAX1(0.0, F11R(I) + F12R(I) * DELTA + F13R(I) * Z)
      F2 = F21R(I) + F22R(I) * DELTA + F23R(I) * Z

      A = (1.0-F1)*( 1.0 + COS(SLOPE)) / 2.0
      B1 = F1*MAX(TRM15,0.0)/ZH
      C=F2*SIN(SLOPE)
      SPER = DRAD*(A+B1+C)   
      
C     print results
      WRITE(*,*) '-------------------------------------------------' 
      WRITE(*,*)'ALL IRRADIANCE VALUES ARE EXPRESSED IN W/m2'
      WRITE(*,*) 'EXTRATERRESTRIAL IRRADIANCE =',ERAD
      WRITE(*,*) 'HOURLY CLEARNESS INDEX, KT =', TKT 
      WRITE(*,*) 'SLOPE BEAM IRRADIANCE =',BSRAD       
      WRITE(*,*) 'SLOPE REFLECTED IRRADIANCE =',RSRAD
      WRITE(*,*) '-------------------------------------------------'                        
      WRITE(*,*) 'ISOTROPIC SKY-DIFFUSE IRRADIANCE =',SISO   
      WRITE(*,*) 'HAY SKY-DIFFUSE IRRADIANCE =',SHAY   
      WRITE(*,*) 'SKARTVEIT-OLSETH SKY-DIFFUSE IRRADIANCE =',SSKOL
      WRITE(*,*) 'REINDL SKY-DIFFUSE IRRADIANCE =',SREIND 
      WRITE(*,*) 'GUEYMARD SKY-DIFFUSE IRRADIANCE=', SGUEYM      
      WRITE(*,*) 'MUNEER SKY-DIFFUSE IRRADIANCE=',SMUN     
      WRITE(*,*) 'PEREZ SKY-DIFFUSE IRRADIANCE=',SPER
	write(*,*) 'Key-in any alphanumeric key to exit'
	read(*,1) Anum

      END   

[cwhanse]

Just a comment, this discussion illustrates the original motivation for pvlib - to clarify what each is in each model. Hang in there @BernatNicolau

[BernatNicolau]

I have created the function based on the 2004 book, but I have not pushed it to this PR because I don't know how to proceed:

• Should I push it in this PR?
• Shall we have two muneer functions muneer1990 and muneer2004, or the Muneer from 1990 can be considered as deprecated and only implement the one from 2004?

This article 10.2790/91554 explains the transposition model used in PVGIS which basically is the one from the book with an improvement for low solar elevations. The function I have created is based on this article.

[adriesse]

I would say the more models, the merrier. There may be some overhead in making them all available in the model chain(s), but that is neither obligatory nor urgent. Having a single muneer function with a selector argument (1990, 2004, 'pvgis') might make sense, but I would be equally happy with three separate functions. Or perhaps two of them call the third to reuse some code.

[AdamRJensen]

• Should I push it in this PR?
• Shall we have two muneer functions muneer1990 and muneer2004, or the Muneer from 1990 can be considered as deprecated and only implement the one from 2004?

Outside of PVGIS I have never heard of the Muneer transposition model, nor have I seen studies where it has an outstanding performance. Given that there are so many transposition models, I think we best only include the more recent Muneer model that is used by PVGIS.

How well does it match PVGIS btw? I have made my own implementation and can get it to match quite well for most times. I'd be happy to collaborate on investigating how well it matches pvgis.

[BernatNicolau]

• Should I push it in this PR?
• Shall we have two muneer functions muneer1990 and muneer2004, or the Muneer from 1990 can be considered as deprecated and only implement the one from 2004?

Outside of PVGIS I have never heard of the Muneer transposition model, nor have I seen studies where it has an outstanding performance. Given that there are so many transposition models, I think we best only include the more recent Muneer model that is used by PVGIS.

How well does it match PVGIS btw? I have made my own implementation and can get it to match quite well for most times. I'd be happy to collaborate on investigating how well it matches pvgis.

I have commited the muneer function from the book (muneer2004), note that the documentation is not updated. Can you replicate the results I obtain from the following code (comparison_allyear.csv) with your own version of the function so that we can compare differences? If you find a better way to collaborate please tell me :)

import pvlib
import pandas as pd

surface_tilt = 20
surface_azimuth = 180
latitude = 37
longitude = -4

# get PVGIS GHI, DHI for one month by requesting the TMY dataset:
pvgis_tmy, months, _, meta = pvlib.iotools.get_pvgis_tmy(
    latitude, longitude, usehorizon=False)
pvgis_all = pd.DataFrame()

for y in range(1, 13):
    print(y)
    pvgis_tmy_year = pvgis_tmy.loc[pvgis_tmy.index.month == y]

    # now get PVGIS POA for the same year/month:
    year = months[y-1]['year']  # the year for the TMY january
    pvgis_poa, _, meta = pvlib.iotools.get_pvgis_hourly(
        latitude, longitude, start=year, end=year, surface_tilt=surface_tilt, surface_azimuth=surface_azimuth, usehorizon=False)
    pvgis_poa = pvgis_poa.resample('H').mean()
    pvgis_poa = pvgis_poa.loc[pvgis_poa.index.month == y]  # keep only january

    # transpose GHI, DHI to POA:
    sp = pvlib.solarposition.get_solarposition(
        pvgis_poa.index, latitude, longitude)

    sp.index = pvgis_poa.index
    dni_extra = pvlib.irradiance.get_extra_radiation(pvgis_poa.index)

    poa_diffuse_muneer = pvlib.irradiance.muneer2004(
        surface_tilt, surface_azimuth, pvgis_tmy_year['dhi'], pvgis_tmy_year['ghi'], dni_extra, solar_azimuth=sp['azimuth'], solar_zenith=sp['zenith'])
    poa_muneer = pvlib.irradiance.get_total_irradiance(
        surface_tilt=surface_tilt,
        surface_azimuth=surface_azimuth,
        solar_zenith=sp['zenith'],
        solar_azimuth=sp['azimuth'],
        dni=pvgis_tmy_year['dni'],
        ghi=pvgis_tmy_year['ghi'],
        dhi=pvgis_tmy_year['dhi'],
        dni_extra=dni_extra,
        model="muneer",
    )

    aux = pd.DataFrame({
        'sky_diff-PVGIS': pvgis_poa['poa_sky_diffuse'],
        'sky_diff-muneer': poa_diffuse_muneer,
        'poa-PVGIS': pvgis_poa.iloc[:, 0:2].sum(axis=1),
        'poa-muneer': poa_muneer['poa_global']
    })
    out = pd.concat([aux, pvgis_tmy_year, sp, dni_extra], axis=1)
    pvgis_all = pd.concat([pvgis_all, out])


pvgis_all.to_csv('comparison_allyear.csv')

[adriesse]

Just wondering what might be the next step on this...?