r.sun(1grass)

r.sun โ€” the solar irradiance & irradiation model ๐Ÿ’– (from the GRASS User's Manual)

okay bestie, sit DOWN because r.sun is literally the sunniest module in all of GRASS and we're about to talk about it. ๐Ÿ’…

NAME โœจ

r.sun โ€” Solar irradiance and irradiation model.

She computes direct (beam), diffuse and reflected solar irradiation raster maps for given day, latitude, surface and atmospheric conditions. and like, the solar parameters (sunrise, sunset times, declination, extraterrestrial irradiance, daylight length โ€” all of it ๐Ÿซถ) get saved in the map history file so you don't lose the receipts.

Alternatively, a local time can be specified to compute solar incidence angle and/or irradiance raster maps. And the shadowing effect of the topography is optionally incorporated, which is giving main character energy fr fr.

KEYWORDS ๐Ÿท๏ธ

raster, solar, sun energy, shadow, parallel


๐Ÿฆ‹ SYNOPSIS ๐Ÿฆ‹

r.sun
r.sun --help
r.sun [-pm] elevation=string  [aspect=string]   [aspect_value=float]
[slope=string]   [slope_value=float]   [linke=string]
[linke_value=float]   [albedo=string]   [albedo_value=float]
[lat=string]   [long=string]   [coeff_bh=string]   [coeff_dh=string]
[horizon_basename=basename]   [horizon_step=float]   [incidout=string]
[beam_rad=string]   [diff_rad=string]   [refl_rad=string]
[glob_rad=string]   [insol_time=string]  day=integer  [step=float]
[declination=float]   [solar_constant=float]   [time=float]
[nprocs=integer]   [distance_step=float]   [npartitions=integer]
[civil_time=float]   [--overwrite]  [--help]  [--verbose]  [--quiet]
[--ui]

๐ŸŽ€ Flags

  • -p โ€” Do not incorporate the shadowing effect of terrain
  • -m โ€” Use the low-memory version of the program
  • --overwrite โ€” Allow output files to overwrite existing files
  • --help โ€” Print usage summary
  • --verbose โ€” Verbose module output
  • --quiet โ€” Quiet module output
  • --ui โ€” Force launching GUI dialog

๐ŸŒธ Parameters

  • elevation=string [required] โ€” Name of the input elevation raster map [meters]
  • aspect=string โ€” Name of the input aspect map (terrain aspect or azimuth of the solar panel) [decimal degrees]
  • aspect_value=float โ€” A single value of the orientation (aspect), 270 is south. Default: 270
  • slope=string โ€” Name of the input slope raster map (terrain slope or solar panel inclination) [decimal degrees]
  • slope_value=float โ€” A single value of inclination (slope). Default: 0.0
  • linke=string โ€” Name of the Linke atmospheric turbidity coefficient input raster map [-]
  • linke_value=float โ€” A single value of the Linke atmospheric turbidity coefficient [-]. Default: 3.0
  • albedo=string โ€” Name of the ground albedo coefficient input raster map [-]
  • albedo_value=float โ€” A single value of the ground albedo coefficient [-]. Default: 0.2
  • lat=string โ€” Name of input raster map containing latitudes [decimal degrees]
  • long=string โ€” Name of input raster map containing longitudes [decimal degrees]
  • coeff_bh=string โ€” Name of real-sky beam radiation coefficient (thick cloud) input raster map [0-1]
  • coeff_dh=string โ€” Name of real-sky diffuse radiation coefficient (haze) input raster map [0-1]
  • horizon_basename=basename โ€” The horizon information input map basename
  • horizon_step=float โ€” Angle step size for multidirectional horizon [degrees]
  • incidout=string โ€” Output incidence angle raster map (mode 1 only)
  • beam_rad=string โ€” Output beam irradiance [W.m-2] (mode 1) or irradiation raster map [Wh.m-2.day-1] (mode 2)
  • diff_rad=string โ€” Output diffuse irradiance [W.m-2] (mode 1) or irradiation raster map [Wh.m-2.day-1] (mode 2)
  • refl_rad=string โ€” Output ground reflected irradiance [W.m-2] (mode 1) or irradiation raster map [Wh.m-2.day-1] (mode 2)
  • glob_rad=string โ€” Output global (total) irradiance/irradiation [W.m-2] (mode 1) or irradiance/irradiation raster map [Wh.m-2.day-1] (mode 2)
  • insol_time=string โ€” Output insolation time raster map [h] (mode 2 only)
  • day=integer [required] โ€” No. of day of the year (1-365). Options: 1-365
  • step=float โ€” Time step when computing all-day radiation sums [decimal hours]. Default: 0.5
  • declination=float โ€” Declination value (overriding the internally computed value) [radians]
  • solar_constant=float โ€” Solar constant [W/m^2]. Default: 1367
  • time=float โ€” Local (solar) time (to be set for mode 1 only) [decimal hours]. Options: 0-24
  • nprocs=integer โ€” Number of threads which will be used for parallel computing. Options: 1-1000. Default: 1
  • distance_step=float โ€” Sampling distance step coefficient (0.5-1.5). Default: 1.0
  • npartitions=integer โ€” Read the input files in this number of chunks. Default: 1
  • civil_time=float โ€” Civil time zone value, if none, the time will be local solar time

๐ŸŒˆ that's just part 1 bestie, there's more coming โ€” don't go anywhere ๐Ÿ’–

๐Ÿ’– DESCRIPTION ๐Ÿ’–

Okay bestie, gather round, because r.sun is that girl ๐Ÿ’…. What she does is compute beam (direct), diffuse and ground reflected solar irradiation raster maps, and she does it for a given day, latitude, surface and atmospheric conditions โœจ. All the solar parameters she uses โ€” think time of sunrise and sunset, declination, extraterrestrial irradiance, daylight length โ€” get stashed away in the resultant maps' history files, so you can peek at them later ๐Ÿซถ.

Alternatively? You can specify the local time and she'll compute a solar incidence angle and/or irradiance raster map instead. Either way, the shadowing effect of the topography is incorporated by default ๐ŸŒ„. She can do that either internally, by calculating the shadowing effect straight off the digital elevation model, or by you passing in raster maps of the horizon height, which is so much faster. Those horizon raster maps? You can make them with r.horizon ๐Ÿฆ‹.


Heads up about units, this part matters ๐Ÿ’€ For latitude-longitude coordinates she requires that the elevation map is in meters. The rules are:
  • lat/lon coordinates: elevation in meters;
  • Other coordinates: elevation in the same unit as the easting-northing coordinates.

The solar geometry of the model is based on the works of Krcho (1990), later improved by Jenco (1992) ๐Ÿ“š. The equations describing Sun -- Earth position, plus the interaction of solar radiation with the atmosphere, were originally based on the formulas suggested by Kitler and Mikler (1986). This component was considerably updated by the results and suggestions of the working group co-ordinated by Scharmer and Greif (2000) โ€” and yes, this algorithm might be replaced by the SOLPOS algorithm-library included in GRASS within the r.sunmask command, iykyk ๐Ÿ‘€.

The model computes all three components of global radiation (beam, diffuse and reflected) for clear sky conditions, i.e. she is not taking into consideration the spatial and temporal variation of clouds โ˜๏ธ๐Ÿ˜ค. The extent and spatial resolution of the modelled area, as well as integration over time, are limited only by the memory and data storage resources. She's built to fulfil user needs in various fields of science (hydrology, climatology, ecology and environmental sciences, photovoltaics, engineering, etc.) for continental, regional up to the landscape scales ๐ŸŒ๐Ÿ‘‘.

The model considers a shadowing effect of the local topography unless you switch it off with the -p flag. And here's the tea: r.sun works in two modes ๐ŸŽ€.

  • Mode one: for the set local time she calculates a solar incidence angle [degrees] and solar irradiance values [W.m-2].
  • Mode two: daily sums of solar radiation [Wh.m-2.day-1] are computed within a set day.

By scripting, the two modes can be used separately or in a combination to provide estimates for any desired time interval โฐ. The model accounts for sky obstruction by local relief features. Several solar parameters are saved in the resultant maps' history files, which may be viewed with the r.info command ๐Ÿ”.


๐Ÿฆ‹ the incidence angle situation

The solar incidence angle raster map incidout is computed specifying elevation raster map elevation, aspect raster map aspect, slope steepness raster map slope, given the day day and local time time โญ.

There is no need to define latitude for projects with a known and defined coordinate system (check it with the g.proj command, bestie ๐Ÿซถ). If you have an undefined projection, an (x,y) system, etc., then the latitude can be defined explicitly for large areas by input raster map lat_in with interpolated latitude values.

IMPORTANT: All input raster maps must be floating point (FCELL) raster maps ๐Ÿ’€.

Null data in maps are excluded from the computation (and also speeding-up the computation, slay โœจ), so each output raster map will contain null data in cells according to all input raster maps. You can use the r.null command to create/reset a null file for your input raster maps.

The specified day day is the number of the day of the general year where January 1 is day no.1 and December 31 is 365 ๐Ÿ“…. Time time must be a local (solar) time โ€” i.e. NOT a zone time, e.g. GMT, CET โ€” in decimal system, e.g. 7.5 (= 7h 30m A.M.), 16.1 = 4h 6m P.M. โฑ๏ธ.

The solar declination parameter is an option to override the value computed by the internal routine for the day of the year ๐ŸŒž. The value of geographical latitude can be set as a constant for the whole computed region or, as an option, a grid representing spatially distributed values over a large region. The geographical latitude must be also in decimal system with positive values for northern hemisphere and negative for southern one ๐ŸŒˆ. In similar principle the Linke turbidity factor (linke, lin) and ground albedo (albedo, alb) can be set ๐ŸŽ€.


โ˜๏ธ real-sky radiation

Besides clear-sky radiations, you can compute a real-sky radiation (beam, diffuse) using coeff_bh and coeff_dh input raster maps, which define the fraction of the respective clear-sky radiations reduced by atmospheric factors (e.g. cloudiness) ๐ŸŒง๏ธ. The value is between 0-1. Usually these coefficients can be obtained from long-terms meteorological measurements provided as raster maps with spatial distribution of these coefficients separately for beam and diffuse radiation (see Suri and Hofierka, 2004, section 3.2) ๐Ÿ“–.

The solar irradiation or irradiance raster maps beam_rad, diff_rad, refl_rad are computed for a given day day, latitude lat_in, elevation elevation, slope slope and aspect aspect raster maps ๐Ÿ”ฅ. For convenience, the output raster given as glob_rad will output the sum of the three radiation components ๐Ÿ’ซ.

The program uses the Linke atmosphere turbidity factor and ground albedo coefficient. A default, single value of Linke factor is lin=3.0 and is near the annual average for rural-city areas ๐ŸŒพ๐Ÿ™๏ธ. The Linke factor for an absolutely clear atmosphere is lin=1.0 โœจ. See notes below to learn more about this factor. The incidence solar angle is the angle between horizon and solar beam vector ๐Ÿ“.


โฐ the step situation (aka how fine is your tea)

The solar radiation maps for a given day are computed by integrating the relevant irradiance between sunrise and sunset times for that day ๐ŸŒ…. You can set a finer or coarser time step used for all-day radiation calculations with the step option. The default value of step is 0.5 hour.

Larger steps (e.g. 1.0-2.0) can speed-up calculations but produce less reliable (and more jagged) results ๐Ÿ“‰.

As the sun moves through approx. 15ยฐ of the sky in an hour, the default step of half an hour will produce 7.5ยฐ steps in the data. For relatively smooth output with the sun placed for every degree of movement in the sky, you should set the step to 4 minutes or less. step=0.05 is equivalent to every 3 minutes โณ. Of course setting the time step to be very fine proportionally increases the module's running time โ€” no free lunch bestie ๐Ÿ˜ญ.

The output units are in Wh per squared meter per given day [Wh/(m*m)/day] ๐Ÿ“Š.

The incidence angle and irradiance/irradiation maps are computed with the shadowing influence of relief by default ๐ŸŒ„. It is also possible for them to be computed without this influence using the planar flag (-p). In mountainous areas this can lead to very different results! โ›ฐ๏ธ๐Ÿ’€ You should be aware that taking into account the shadowing effect of relief can slow down the speed of computation, especially when the sun altitude is low ๐ŸŒ‡.


๐ŸŒธ distance_step, the sampling density queen

When considering the shadowing effect, speed and precision of computation can be controlled by the distance_step parameter, which defines the sampling density at which the visibility of a grid cell is computed in the direction of a path of the solar flow ๐Ÿ”ญ. It also defines the method by which the obstacle's altitude is computed. Here's how it breaks down:

  • When choosing a distance_step less than 1.0 (i.e. sampling points will be computed at distance_step * cellsize distance), r.sun takes the altitude from the nearest grid point ๐Ÿ“.
  • Values above 1.0 will use the maximum altitude value found in the nearest 4 surrounding grid points โ›ฐ๏ธ.

The default value distance_step=1.0 should give reasonable results for most cases (e.g. on DEM) โœจ. The distance_step value defines a multiplying coefficient for sampling distance ๐Ÿ“. This basic sampling distance equals to the arithmetic average of both cell sizes. The reasonable values are in the range 0.5-1.5 ๐ŸŽฏ. Values below 0.5 will decrease and values above 1.0 will increase the computing speed ๐Ÿš€. Values greater than 2.0 may produce estimates with lower accuracy in highly dissected relief ๐Ÿ’€.

The fully shadowed areas are written to the output maps as zero values ๐Ÿ…พ๏ธ. Areas with NULL data are considered as no barrier with shadowing effect ๐Ÿ‘ป.


๐Ÿ“œ what lands in the history files

The maps' history files are generated containing the following listed parameters used in the computation:

  • Solar constant used W.m-2 โ˜€๏ธ
  • Extraterrestrial irradiance on a plane perpendicular to the solar beam [W.m-2] โœจ
  • Day of the year ๐Ÿ“…
  • Declination [radians] ๐Ÿ“
  • Decimal hour (Alternative 1 only) ๐Ÿ•
  • Sunrise and sunset (min-max) over a horizontal plane ๐ŸŒ…
  • Daylight lengths โณ
  • Geographical latitude (min-max) ๐ŸŒ
  • Linke turbidity factor (min-max) ๐ŸŒซ๏ธ
  • Ground albedo (min-max) ๐Ÿค

๐ŸŽ€ one last vibe

You can use a nice shellcript with variable day to compute radiation for some time interval within the year (e.g. vegetation or winter period) ๐Ÿƒโ„๏ธ. And listen: Elevation, aspect and slope input values should not be reclassified into coarser categories. This could lead to incorrect results ๐Ÿ’€๐Ÿšซ. Don't do it, bestie, I'm begging ๐Ÿ™.


โœจ OPTIONS โœจ

Currently, there are two modes of r.sun ๐Ÿ‘‘.

  • In the first mode it calculates solar incidence angle and solar irradiance raster maps using the set local time โฐ.
  • In the second mode daily sums of solar irradiation [Wh.m-2.day-1] are computed for a specified day ๐ŸŒž.

๐Ÿ“ NOTES ๐Ÿ“

Solar energy is such an important input parameter in all kinds of models โ€” energy industry, landscape, vegetation, evapotranspiration, snowmelt, remote sensing, the whole lineup ๐Ÿ’…. And like, solar rays incidence angle maps can be effectively used in radiometric and topographic corrections in mountainous and hilly terrain, aka the places where you really need your investigations to be very accurate, bestie. No pressure.

The clear-sky solar radiation model applied in the r.sun is based on the work undertaken for development of European Solar Radiation Atlas (Scharmer and Greif 2000, Page et al. 2001, Rigollier 2001) ๐Ÿซถ. This clear sky model estimates the global radiation from the sum of its beam, diffuse and reflected components. The main difference between solar radiation models for inclined surfaces in Europe is how they treat the diffuse component โ€” and in the European climate this component is often the largest source of estimation error, which is a lot ๐Ÿ˜ญ. So taking into consideration the existing models and their limitation, the European Solar Radiation Atlas team selected the Muneer (1990) model as it has a sound theoretical basis and thus more potential for later improvement. Iconic choice, honestly.

Details of underlying equations used in this program can be found in the reference literature cited below or in the book published by Neteler and Mitasova: Open Source GIS: A GRASS GIS Approach (published in Kluwer Academic Publishers in 2002) ๐Ÿ“š.

Average monthly values of the Linke turbidity coefficient for a mild climate in the northern hemisphere (pls go see the reference literature for your study area, this is not universal bestie):

  • mountains โ€” Jan 1.5, Feb 1.6, Mar 1.8, Apr 1.9, May 2.0, Jun 2.3, Jul 2.3, Aug 2.3, Sep 2.1, Oct 1.8, Nov 1.6, Dec 1.5, annual 1.90
  • rural โ€” Jan 2.1, Feb 2.2, Mar 2.5, Apr 2.9, May 3.2, Jun 3.4, Jul 3.5, Aug 3.3, Sep 2.9, Oct 2.6, Nov 2.3, Dec 2.2, annual 2.75
  • city โ€” Jan 3.1, Feb 3.2, Mar 3.5, Apr 4.0, May 4.2, Jun 4.3, Jul 4.4, Aug 4.3, Sep 4.0, Oct 3.6, Nov 3.3, Dec 3.1, annual 3.75
  • industrial โ€” Jan 4.1, Feb 4.3, Mar 4.7, Apr 5.3, May 5.5, Jun 5.7, Jul 5.8, Aug 5.7, Sep 5.3, Oct 4.9, Nov 4.5, Dec 4.2, annual 5.00

Planned improvements include the use of the SOLPOS algorithm for solar geometry calculations and internal computation of aspect and slope. The future is coming โœจ.


โฐ Solar time โฐ

By default r.sun calculates times as true solar time, whereby solar noon is always exactly 12 o'clock everywhere in the current region. Which sounds cute and simple until you realize that depending on where the zone of interest is located in the related time zone, this may cause differences of up to an hour โ€” and in some cases (like Western Spain) even more ๐Ÿ™ˆ. On top of this, the offset varies during the year according to the Equation of Time. Not me getting personally victimized by orbital mechanics.

To overcome this problem, the user can use the option civil_time=<timezone_offset> in r.sun to make it use real-world (wall clock) time ๐Ÿ•ฐ๏ธ. For example, for Central Europe the timezone offset is +1, +2 when daylight saving time is in effect.


๐ŸŒ‘ Extraction of shadow maps ๐ŸŒ‘

A map of shadows can be extracted from the solar incidence angle map (incidout). Areas with NULL values are shadowed ๐Ÿ”ฅ. Heads up though: this will not work if the -p flag has been used.


๐Ÿง  Large maps and out of memory problems ๐Ÿง 

With a large number or columns and rows, r.sun can consume a significant amount of memory ๐Ÿ˜ฎโ€๐Ÿ’จ. While output raster maps are not partitionable, the input raster maps are using the npartitions parameter. In case of out of memory error (ERROR: G_malloc: out of memory), the npartitions parameter can be used to run a segmented calculation which consumes less memory during the computations. She's efficient, she just needs to breathe.

The amount of memory by r.sun is estimated as follows:

# without input raster map partitioning:
#  memory requirements: 4 bytes per raster cell
#  rows,cols: rows and columns of current region (find out with g.region)
#  IR: number of input raster maps without horizon maps
#  OR: number of output raster maps
memory_bytes = rows*cols*(IR*4 + horizon_steps + OR*4)
# with input raster map partitioning:
memory_bytes = rows*cols*((IR*4+horizon_steps)/npartitions  + OR*4)

๐ŸŽ€ EXAMPLES ๐ŸŽ€

North Carolina example (considering also cast shadows):

g.region raster=elevation -p
# calculate horizon angles (to speed up the subsequent r.sun calculation)
r.horizon elevation=elevation step=30 bufferzone=200 output=horangle \
    maxdistance=5000
# slope + aspect
r.slope.aspect elevation=elevation aspect=aspect.dem slope=slope.dem
# calculate global radiation for day 180 at 2p.m., using r.horizon output
r.sun elevation=elevation horizon_basename=horangle horizon_step=30 \
      aspect=aspect.dem slope=slope.dem glob_rad=global_rad day=180 time=14
# result: output global (total) irradiance/irradiation [W.m-2] for given day/time
r.univar global_rad

Calculation of the integrated daily irradiation for a region in North-Carolina for a given day of the year at 30m resolution. Here day 172 (i.e., 21 June in non-leap years):

g.region raster=elev_ned_30m -p
# considering cast shadows
r.sun elevation=elev_ned_30m linke_value=2.5 albedo_value=0.2 day=172 \
      beam_rad=b172 diff_rad=d172 \
      refl_rad=r172 insol_time=it172
d.mon wx0
# show irradiation raster map [Wh.m-2.day-1]
d.rast.leg b172
# show insolation time raster map [h]
d.rast.leg it172

We can compute the day of year from a specific date in Python ๐Ÿ:

>>> import datetime
>>> datetime.datetime(2014, 6, 21).timetuple().tm_yday
172

๐Ÿ‘€ SEE ALSO ๐Ÿ‘€

r.horizon, r.slope.aspect, r.sunhours, r.sunmask, g.proj, r.null, v.surf.rst ๐Ÿ’–

๐Ÿ“š REFERENCES ๐Ÿ“š

okay bestie, this is the reading list. the receipts. if you ever wanna go "wait, where did this model even come from", it's all right here, and yes, we're citing every single one of them ๐ŸŽ€

  • Hofierka, J., Suri, M. (2002): The solar radiation model for Open source GIS: implementation and applications. International GRASS users conference in Trento, Italy, September 2002. (PDF)
  • Hofierka, J. (1997). Direct solar radiation modelling within an open GIS environment. Proceedings of JEC-GIโ€™97 conference in Vienna, Austria, IOS Press Amsterdam, 575-584.
  • Jenco, M. (1992). Distribution of direct solar radiation on georelief and its modelling by means of complex digital model of terrain (in Slovak). Geograficky casopis, 44, 342-355.
  • Kasten, F. (1996). The Linke turbidity factor based on improved values of the integral Rayleigh optical thickness. Solar Energy, 56 (3), 239-244.
  • Kasten, F., Young, A. T. (1989). Revised optical air mass tables and approximation formula. Applied Optics, 28, 4735-4738.
  • Kittler, R., Mikler, J. (1986): Basis of the utilization of solar radiation (in Slovak). VEDA, Bratislava, p. 150.
  • Krcho, J. (1990). Morfometrickรก analza a digitรกlne modely georeliรฉfu (Morphometric analysis and digital models of georelief, in Slovak). VEDA, Bratislava.
  • Muneer, T. (1990). Solar radiation model for Europe. Building services engineering research and technology, 11, 4, 153-163.
  • Neteler, M., Mitasova, H. (2002): Open Source GIS: A GRASS GIS Approach, Kluwer Academic Publishers. (Appendix explains formula; r.sun script download)
  • Page, J. ed. (1986). Prediction of solar radiation on inclined surfaces. Solar energy R&D in the European Community, series F - Solar radiation data, Dordrecht (D. Reidel), 3, 71, 81-83.
  • Page, J., Albuisson, M., Wald, L. (2001). The European solar radiation atlas: a valuable digital tool. Solar Energy, 71, 81-83.
  • Rigollier, Ch., Bauer, O., Wald, L. (2000). On the clear sky model of the ESRA - European Solar radiation Atlas - with respect to the Heliosat method. Solar energy, 68, 33-48.
  • Scharmer, K., Greif, J., eds., (2000). The European solar radiation atlas, Vol. 2: Database and exploitation software. Paris (Les Presses de lโ€™ร‰cole des Mines).
  • Joint Research Centre: GIS solar radiation database for Europe and Solar radiation and GIS

๐Ÿ‘‘ AUTHORS ๐Ÿ‘‘

the literal icons who made this happen, put some respect on the names โœจ

  • Jaroslav Hofierka, GeoModel, s.r.o. Bratislava, Slovakia
  • Marcel Suri, GeoModel, s.r.o. Bratislava, Slovakia
  • Thomas Huld, JRC, Italy

ยฉ 2007, Jaroslav Hofierka, Marcel Suri. This program is free software under the GNU General Public License (>=v2) hofierka@geomodel.sk suri@geomodel.sk

๐Ÿ’– free software fr fr. GPL >=v2. now go email them if you build something cool.

๐ŸŒˆ SOURCE CODE ๐ŸŒˆ

Available at: r.sun source code (history)

Latest change: Friday Apr 17 01:04:26 2026 in commit: 32a143358df26b2e518b03c7e0d79541e8d1e5b1

  • Main index
  • Raster index
  • Topics index
  • Keywords index
  • Graphical index
  • Full index

ยฉ 2003-2026 GRASS Development Team, GRASS 8.5.0 Reference Manual


aaaand that's the whole page bestie, you ate it ๐Ÿ’… go model some sunshine, and if anybody asks, you read the docs. main character energy only ๐Ÿ‘‘

this was GRASS 8.5.0, page r.sun(1grass), last touched April 17, 2026 โ€” so like, lowkey fresh out the oven ๐Ÿ”ฅ๐Ÿฆ‹ k bye ๐Ÿซถ

๐Ÿ“ฆ where this page comes from