Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Muneer transposition model #2184

Draft
wants to merge 55 commits into
base: main
Choose a base branch
from
Draft

Conversation

BernatNicolau
Copy link
Contributor

@BernatNicolau BernatNicolau commented Aug 26, 2024

  • Closes to Add Hofierka 2002 model #2117 (this comment)
  • I am familiar with the contributing guidelines
  • Tests added
  • Updates entries in docs/sphinx/source/reference for API changes.
  • 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`).
  • New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
  • Pull request is nearly complete and ready for detailed review.
  • Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

still not tested
Copy link
Member

@AdamRJensen AdamRJensen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

An initial review

pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
Comment on lines 1030 to 1034
.. [1] Muneer, T., 1990, Solar radiation model for Europe.
Building services engineering research and technology, 11: 153-163.

.. [2] Moon P and Spencer D E Illumination from a non-uniform sky
Trans. Illum. Eng. Soc. (London) 37 707-725 (1942)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These two sources follow different citation convention.

Also, could you add a DOI?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Resolved in 8b13595 . First time I do it, @AdamRJensen can you check if it is alright?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@BernatNicolau both look typed out correctly, but I don't think the second DOI value is correct (review comment)

BernatNicolau and others added 10 commits August 26, 2024 16:15
Co-authored-by: Adam R. Jensen <[email protected]>
Co-authored-by: Adam R. Jensen <[email protected]>
Adam initial review

Co-authored-by: Adam R. Jensen <[email protected]>
Adam initial review

Co-authored-by: Adam R. Jensen <[email protected]>
Adam initial review

Co-authored-by: Adam R. Jensen <[email protected]>
Adam initial review

Co-authored-by: Adam R. Jensen <[email protected]>
Copy link
Contributor

@echedey-ls echedey-ls left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks really good! Some minor documentation comments down below.

Since this is a new feature, we have to do some additional changes to the code:

  1. Docs. I suppose you have tried to find your function at the built readthedocs page, but had no luck. That happens because each function we want to publicity show in the API needs to be listed in some index pages. The appropriate list for this function is in docs/sphinx/source/reference/irradiance/transposition.rst
  2. Tests. We want to ensure this function works as expected in the future, so tests is a must-have. Look for test_irradiance.py and try to add a test. You can just mimic current code or read more about pytest. If you need some help with that let us know. By the way, the best way I have found to make test data is by using a spreadsheet, so you can double check you implementation.

At some point a whatsnew entry will need to be made. Best is to defer that for a bit so we can have more reviews till then.

pvlib/irradiance.py Outdated Show resolved Hide resolved
Comment on lines 1005 to 1007
Surface tilt angle in decimal degrees. Tilt must be >=0 and
<=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Surface tilt angle in decimal degrees. Tilt must be >=0 and
<=180. The tilt angle is defined as degrees from horizontal
(e.g. surface facing up = 0, surface facing horizon = 90)
Surface tilt angle from horizontal. E.g.,
facing up = 0°, facing horizon = 90°.
In degrees [°].

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Got the wording from perez. For consistency, should other functions also be updated or shall I keep the one that currently is (from perez)?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

echedey-ls edits are an improvement but I'm not in favor of changing documentation for other functions. You can change it here, but I think it's not a big deal if we use the old language.

pvlib/irradiance.py Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Show resolved Hide resolved
@BernatNicolau
Copy link
Contributor Author

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
Copy link
Contributor

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
Copy link
Member

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
Copy link
Contributor Author

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
Copy link
Member

adriesse commented Aug 27, 2024

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
Copy link
Contributor Author

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
Copy link
Member

cwhanse commented Aug 27, 2024

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

@BernatNicolau
Copy link
Contributor Author

@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.

muneer included in get_sky_diffuse code and documentation
muneer included in get_total_irradiance documentation
muneer function updated

TODO: tests, create documentation
added documentation, tests and modifications according to flake8
still quite new to the flake8. May not pass the test
@BernatNicolau
Copy link
Contributor Author

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 #2117.

This has solved the negative spikes:
Figure_1
Figure_2

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.

@BernatNicolau BernatNicolau marked this pull request as ready for review September 3, 2024 19:16
pvlib/irradiance.py Outdated Show resolved Hide resolved
Copy link
Contributor

@RDaxini RDaxini left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

a few additional suggestions

pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
pvlib/irradiance.py Outdated Show resolved Hide resolved
@kandersolar
Copy link
Member

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
Copy link
Contributor Author

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.
image

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

@adriesse
Copy link
Member

adriesse commented Sep 4, 2024

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
Copy link
Member

cwhanse commented Sep 4, 2024

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

@BernatNicolau
Copy link
Contributor Author

BernatNicolau commented Sep 5, 2024

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
Copy link
Member

adriesse commented Sep 5, 2024

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.

@kandersolar kandersolar modified the milestones: v0.11.1, v0.11.2 Sep 23, 2024
@AdamRJensen
Copy link
Member

AdamRJensen commented Sep 23, 2024

  • 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.

@RDaxini RDaxini mentioned this pull request Sep 26, 2024
5 tasks
@BernatNicolau BernatNicolau marked this pull request as draft September 27, 2024 07:06
@BernatNicolau
Copy link
Contributor Author

  • 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')

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

8 participants