A common question concerning the safety of photovoltaic (PV) power systems is the impact of reflected sunlight. PV modules have the potential to impact neighboring structures or activities, notably aviation. It is important to know where the reflected light will go and what the intensity of the light will be at any point in time.

Reflected sunlight is known as glare or glint. This post is a simple case study for assessing the glare impact from a PV installation in Qatar. It outlines the basic concepts behind glare assessments and provides R source code to facilitate glare studies.

**Sunlight Reflection and Best-Practice Regulations**

**Sunlight Reflection and Best-Practice Regulations**

In general, there are two types of glare created by any kind of reflective surface:

*Specular reflections*are the concentrated glare created by a direct sunbeam. Specular reflections have a glare characteristic similar to that of a mirror and has high sunlight intensity;*Diffuse reflections*represent reflected light that is scattered in all directions and is typically caused by uneven surfaces, such as rough water, gravel or dust. In general, diffuse reflections are low in intensity.

Guidelines for managing glare risks focus exclusively on specular reflections.^{1}

**Project Definition and Glare Parameters**

**Project Definition and Glare Parameters**

A 2 MW PV solar pilot project is currently planned to be built by Qatar Petroleum in Ras lafan Industrial City. The project has the following attributes:

- A maximum of 28,000 square meters (sqm) will be consumed by the installation of 1 MW of crystalline modules and 1 MW of thin film modules. Actual land use will depend on the solar conversion efficiency of the modules actually purchased.
- An additional 4,000 sqm will be consumed by open spaces between the panel rows and other balance of plant needs.
- In total, the project will have a maximum of 3,380 crystalline modules and 10,080 thin film modules, though actual count could due less.
- Using the QND95 coordinate system of Qatar, the front row or string of the PV modules will have a center point located at 222576 (Easting) and 464285 (Northing) . This point is located 192 meters directly North of the main office tower of the Multi-Purpose Administration Complex (MPAC) in Ras Lafan Industrial City.
- The layout of the solar PV array defines an optimal module tilt of 20-degrees, facing due South, to maximize annual energy production. The modules face directly into the windowed front of the MPAC office tower.
- The building plans for the main office tower confirm that the North wall facing the PV power plant has a center point located at 222683 (Easting) and 464045 (Northing). The planned height of the building 48.5 meters.

The project outline is shown below in blue and the MPAC outline is shown in red. Simplified balance of plant cabling is also shown with one cable (green) feeding the MPAC and a second (black) connecting with the MPAC grid substation.

**Summary of Findings**

**Summary of Findings**

*Reflected Light
*There is zero potential that a direct sunlight beam will be reflected by the solar PV modules and contact the windowed front of the MPAC office tower at any time during the year. Secular reflections directed south of the project total 401 years per year, between early April and mid-Setember. However, the angle of the reflected sunbeam is so steep the MPAC office tower is unaffected. Meanwhile, the intensity of the reflected light at its peak will never exceed 12.2 W/m

^{2}, assuming a reflection coefficient of 2% and the solar panels are clean.

*Deflected Light
*For the balance of the daylight production year, after the 401 reflected hours are accounted for, there are 4,406 production hours where specular reflections will be deflected away from the building and will be projected to the North in the direction of open sand and water spaces. The intensity of these specular reflections have a range starting at zero and a maximum of 5 W/m

^{2}, assuming a 2% reflection coefficient and the panels remain clean.

*Sunlight Intensity and Dust
*Solar panels in the GCC countries accumulate dust immediately after cleaning, which will degrade the amount of radiation striking the reflective glass of the PV modules. Dust accumulations will also scatter the light received and reduce the potential for direct sunbeam reflections. Heavy soiling will eliminate reflections completely. Finally, the glass used on the solar modules is designed to minimize sunlight reflections to less than 2% of the irradiance received under standard test conditions. Certification to this effect is necessary to maximize energy capture, power production and human safety.

**Sun Elevation and Reflected Light Angles**

**Sun Elevation and Reflected Light Angles**

The Sun elevation angle is the angle formed by the flat ground and a sunlight beam emitted by the sun. The Sun elevation angle starts at zero degrees in the East, it rises to a near vertical angle of 90-degrees at noon, and declines again to zero degrees in the West. The maximum angle at the project site is 87.47-degrees on June 20^{th} of each year.

Elevation angle defines the angle of specular reflections. Specifically, the sunlight beam striking the surface of a module is either deflected away from the Sun (to the North) or reflected back toward the Sun (to the South). The reflection/deflection boundary is defined by a line perpendicular to the module surface, as shown in the two images below.

A module set to a 20-degree tilt has a deflection/reflection boundary at 70-degree, as confirmed above using the angles highlighted in red and as defined by the horizontal plane. The angle of the incoming sunbeam is measures measured relative to the reflection/deflection boundary. The same angle is then applied on the opposite side of the boundary to determine the direction of the outgoing sun beam. The images make clear that sunlight will generate specular reflections directed toward the Northern building tower only when the Sun is between 70-degrees and 90-degrees.

**RLIC Sun Elevation and MPAC Reflection Angles**

**RLIC Sun Elevation and MPAC Reflection Angles**

The following chart defines the sun elevation angle (y-axis) at the project site as a function of time of year (x-axis) and time of day (colored lines). This data was generated using the solaR package for R, which contains several algorithms and open source code for Sun position tracking. Alternatively, C-code of for solar tracking can be obtained here.

Direct specular reflections toward the MPAC building tower only occur when the elevation angle is greater than 70-degrees. The highlighted box in the image at left (click to enlarge) confirms this is possible only between 11:00 and 13:00 hours during the months of April to September. In total, there are 401 hours during the year when the project has potential to reflect light toward the MPAC building.

Reflected sunlight directed toward the MPAC office tower will leave the module array at angles ranging from 53-degrees (when the sun is almost directly overhead at 87.47-degrees) to 70-degrees (when the Sun is on the deflection/reflection boundary). These are very steep angles given the distance between the solar project and the office tower. By the time the specular reflection arrives at the coordinates of the North wall of the tower, they have an elevation ranging between 250 and 500 meters above the ground. The estimated height is calculated using basic trigonometry:

(1)

where ref.height is the reference height of the reflected sunlight beam at the front wall of the building, DRB is the angle of the deflection/reflection boundary in radians, Elev is the Sun elevation angle in radians, and setback is the distance between the first string of modules closets to the building and the front wall of the building.

The estimated reference height as a function of time is shown in the chart above. The estimated height of the reflections is significantly higher than the 6-story building at 48.5 meters. Direct contact of the sunlight beam would only be possible if the building ranged between 30 and 60 floors.

**Glare Intensity**

**Glare Intensity**

Estimated glare intensity as a function of time is shown in the chart below. The estimation relies on 3 inputs:

- Extra-terrrestrial radiation by hour and adjusted for time of year
- A “quick-n-dirty” scalar to convert total solar irradiance to ground irradiance. The scalar is an average of the ratio of ground irradiance observed at the project site divided by total solar irradiance.
- The glare coefficient or percent of irradiance reflected by the PV modules and sourced from the manufacturer. It is assumed that the glare coefficient is constant as a function of elevation angle.

The maximum glare estimated is 12.2 W/m^{2}. This exceeds the US aviation threshold of 7 W/m^{2}. The estimate suggests that project glare could impact pilot safety and cause a vision “after-image.” In practice, the estimated glare is less than the glare expected from structural glass on office towers. Hence, well established risk mitigation would include the use of low cost glare-resistant sunglasses or visors. The estimated glare also suggests that site data collection is justified once the project is constructed to validate the accuracy of the desk-top estimates.

**Research on Glare Intensity**

**Research on Glare Intensity**

Glare intensity can vary depending on the PV module and module glass used. It can also vary given local atmospheric conditions and dust accumulation on the PV panels. As a result, glare intensity in practice can vary widely across projects and locations.

Numerous aviation studies have been conducted to assess the intensity of reflected sunlight on the human eye for purposes of reinforcing aviation safety. One report, *A Study of the Harzardous Glare Potential to Aviators from Utility-Scale Flat-Panel Photovoltaic Systems*, undertook research at the Las Vegas, Nevada airport. The study relied on ocular safety metrics, which quantify the potential for a post-flash glare after-image on the back of the human eye. The study presented the graph above.

The graph is a source of a potential concern since the percent of reflected light for elevation angles greater than 70 degrees is very high and requires on-site validation. However, glare intensity is shown to be low for the majority of the Sun’s path through the sky. The study concluded the following:

- Specular reflections from solar panels varied with elevation angle, with maximum glare at angles close to 90 degrees;
- The potential for hazardous glare is low and similar to that of smooth water;
- Other common materials presented greater risk to aviators, including white Portland cement, snow and the structured glass of building towers.

The U.S. Federal Aviation Administration (FAA) has also conducted comparative studies of the specular reflection from PV modules.^{2} The following graphic on the relative reflectivity of different surfaces was sourced from the FAA. It is interesting to note that reflections from vegetation are deemed a greater hazard than solar PV panels by the federal aviation regulator.

In practice, aviation safety risk cannot be ignored. *However,* *glare risk to aviators is easily mitigated at low cost through the use of protective eyewear, such as glare resistant sunglasses and visors.* Secondary precautions, if required, include limiting flight path activity for very short and discrete periods during the year.

Finally, another study presented the following relative reflectivity and presented the following graph.^{3}. Again, independent research confirms that reflections from solar PV modules are possible, but less than the risk posed by steel structures, water or office tower glass.

All 3 studies also make clear that the reflective glare does not have one simple scalar that defines glare intensity. Variation in study results are a function of the different PV modules tested and differences in measurement methods. As a result, projects should rely on EPC specifications for equipment purchasing that have no glare potential or the lowest glare potential possible.

**Conclusion**

**Conclusion**

Relocation of the PV project is not justified based on the potential of a sunlight glare impact or safety risk to the MPAC office tower. Published aviation studies suggest that the risk to pilots is also low and easily mitigated at low cost. The long-term health and safety focus of the project should include activities to measure the specular reflections generated by the solar modules to ensure minimal aviation risk and proper precautions. Data collection is required for Qatar specific conditions and results should be incorporated into front-end design engineering as a core standard for PV modules.

**R Source Code**

**R Source Code**

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 |
# Header ---- # # Name: GlareStudy.R # Title: Estimation of sunlight glare reflected by solar PV panels # Version: 1.0 # Date: 2014-Sep-24 # Author: Brad Horn # License: GPL (>=2) # # Description: The script is used to estimate the glare impact from solar PV modules on an # building. The program is used to determine (1) if the height of the reflected # sunbeam at the building location is higher or lower than the height of the # building; (2) the irrqdiance intensity of the reflected sunbeam (W/m^2); and # (3) if the azimuth angle of the Sun confirms of eliminates the glare risk. # # The program requires the following inputs: # - the geapgraphic coordinates of the 4 corners for the solar array (QND95) # - the geogrphic coordnates of the 4 corners of the building (QND95) # - the max building elevation or height above the ground # - the number of reflection points to assess across the front of the array and # - the number of reflection receipt points across the fron of the building # - the start and end dates of the assessment period for solar postion # # if the number of reflection points is set to 1, then the glare assessments is # done using a module in the middle of the front row of the solar array closest to # the building. If the number of points is set to 2, then the glare assessment is # done for the first and last module in the front row of the solar array. Finally, # if the number of reflection points is greater than 2, then the first and last # module in the front row are used and the remaining points are equally spaced. # # The Sun position calculation relies on zenith angle only to estimate the height # of the reflected sunbeam on the front of the building. Yes/no flags test if # the reflected height at the location of the building is less than the height of # the building. # # if a glare risk exists, then the reflected irradiance intensity on the building # is estimated using a glare coefficient. The coefficient defines the percent of # light recieved by the solar module that is then reflected into space. Irradiance # hitting the plane of array is sourced from an external file, which contains a # TMY profile similiar to that used by PVsyst. # # Finally, if the glare risk striking the building is confirmed to be high, then # a third estimation is done using the solar azimuth angle to make sure the # reflected sunbeam is actually hitting the building or not. # # Details: The glare coefficient is assumed to be fixed, regardless of zenith angle. See # OEM module specs to define the amount of irradiance relfected by any particular # module. The code also calculates extra-terrestrial irradiance for any date/time. # The variable "ground.SolI" then converts total solar irradiance to an estimate of # ground-based irradiance. The scalar used for the conversion is an average for # the day and based on actual ground measurement. The resulting irradiance data # is "quick-n-dirty." When glare risk is found to be a real event for potential # safety, then actual irradiance profile should be used. # Dev.Notes: NA # Depends solaR package for R (http://cran.r-project.org/web/packages/solaR/index.html) # Internal functions for unit conversion are defined at the end of the script # References: NA ##-------------------------------------------------------------------------------------------## # 1.Input Data ---- # define the calendar start and end dates for the solar position and angle estimates start.date <- as.POSIXct("2015-01-01 01:00:00", tz = "Asia/Qatar") end.date <- as.POSIXct("2015-12-31 23:00:00", tz = "Asia/Qatar") # define the height and tilt angle of the array modules (meters; degrees) mod.height <- 2 mod.tilt <- 20 # define the 4 corners of the solar array and the building proj.4corners <- data.frame(Points = 1:4, Easting = c(222419.036, 222588.892, 222686.450, 222464.591), Northing = c(464383.242, 464481.466, 464281.824, 464288.849)) bldg.4corners <- data.frame(Points = 1:4, Easting = c(222622.222, 222744.443, 222738.868, 222622.22), Northing = c(464057.895, 464031.579, 464015.790, 464015.79)) # define the two that define the front row of the modules and the building surface # NOTE: The order of points needs to align the closest corners of the array and building # if the vectors are (3, 4) and (2, 1) the points 3/2 and 4/1 are being matched up proj.points <- c(3, 4) bldg.points <- c(2, 1) # define the building height and the number of reflection source/receipt points to assess bldg.height <- 48.25 n.points <- 1 # define glare intensity using a glare coefficient or the % of sunlight expected to be reflected. # Assumes fixed relfectivity, regardless of zenith angle. Define the scalar used to convert # extraterrestrial irradiance to ground irradiance glare.coef <- 0.02 ground.SolI <- 0.46 # 2.Coordinate Data Transforms and Definitions ---- # convert QND95 4 corner coordinates to WGS84 decimal degree coordinates and return a dataframe proj.4corners <- as.data.frame(QND2Dec(proj.4corners)) bldg.4corners <- as.data.frame(QND2Dec(bldg.4corners)) # define all assessment points along the front of the array and building if(n.points == 1){ proj.xy <- data.frame(Points = 1, Longitude = proj.4corners[proj.points[2], "Longitude"] + 0.5 * (proj.4corners[proj.points[1], "Longitude"] - proj.4corners[proj.points[2], "Longitude"]), Latitude = proj.4corners[proj.points[1], "Latitude"] + 0.5 * (proj.4corners[proj.points[2], "Latitude"] - proj.4corners[proj.points[1], "Latitude"])) bldg.xy <- data.frame(Points = 1, Longitude = bldg.4corners[bldg.points[1], "Longitude"] + 0.5 * (bldg.4corners[bldg.points[2], "Longitude"] - bldg.4corners[bldg.points[1], "Longitude"]), Latitude = bldg.4corners[bldg.points[2], "Latitude"] + 0.5 * (bldg.4corners[bldg.points[1], "Latitude"] - bldg.4corners[bldg.points[2], "Latitude"])) } if(n.points == 2){ proj.xy <- data.frame(Points = 1:n.points, Longitude = proj.4corners[proj.points, "Longitude"], Latitude = proj.4corners[proj.points, "Latitude"]) bldg.xy <- data.frame(Points = 1:n.points, Longitude = bldg.4corners[bldg.points, "Longitude"], Latitude = bldg.4corners[bldg.points, "Latitude"]) } if(n.points > 2){ proj.xy <- data.frame(Points = 1:n.points, Longitude = seq(proj.4corners[proj.points[2], "Longitude"], proj.4corners[proj.points[1], "Longitude"], length.out = n.points), Latitude = seq(proj.4corners[proj.points[1], "Latitude"], proj.4corners[proj.points[2], "Latitude"], length.out = n.points)) bldg.xy <- data.frame(Points = 1:n.points, Longitude = seq(bldg.4corners[bldg.points[1], "Longitude"], bldg.4corners[bldg.points[2], "Longitude"], length.out = n.points), Latitude = seq(bldg.4corners[bldg.points[2], "Latitude"], bldg.4corners[bldg.points[1], "Latitude"], length.out = n.points)) } # define the geographic setback distance between the array and the building for all assessment points setback <- 1000 * geodist(proj.xy[, "Longitude"], proj.xy[, "Latitude"], bldg.xy[, "Longitude"], bldg.xy[, "Latitude"]) # 2. Sun Position Calculations ---- # fBTd: sequence of days for the assessment year temp.year <- year(start.date) temp.day.seq <- fBTd(mode = "serie", year = temp.year) - (3 * 60 * 60) attr(temp.day.seq, "tzone") <- "Asia/Qatar" # get PV module latitude temp.lat <- proj.xy[, "Latitude"] # fSolID - daily movements of the Sun and Earth as seen for target latitude temp.solD <- fSolD(temp.lat, temp.day.seq) # fSolI - solar angles by hour (radians) temp.sort <- c("aman", "AzS", "AlS", "Bo0") temp.solI <- fSolI(temp.solD, sample = "hour")[, temp.sort] # convert solar angles to degrees temp.solI[, c("AzS", "AlS")] <- R2D(temp.solI[, c("AzS", "AlS")]) # clean the solar angle data temp.solI[, "AzS"] <- replace(temp.solI[, "AzS"], temp.solI[, "aman"] == 0, NA) temp.solI[, "AlS"] <- replace(temp.solI[, "AlS"], temp.solI[, "aman"] == 0, NA) temp.solI[, "Bo0"] <- replace(temp.solI[, "Bo0"], temp.solI[, "aman"] == 0, NA) solI.angles <- data.frame(TimeStamp = time(temp.solI), Month = as.factor(month(time(temp.solI))), Hour = as.factor(hour(time(temp.solI))), Viz.TF = coredata(temp.solI[, "aman"]), Azimuth = coredata(temp.solI[, "AzS"]), Elevation = coredata(temp.solI[, "AlS"]), Zenith = 90 - coredata(temp.solI[, "AlS"]), SolI = ground.SolI * coredata(temp.solI[, "Bo0"])) # determine the day of year with the min/max elevation angle at 12 noon noon.angles <- solI.angles[which(hour(solI.angles[, 1]) == 12), ] noon.max <- max(noon.angles[,"Elevation"]) noon.min <- min(noon.angles[,"Elevation"]) noon.max.date <- noon.angles[which(noon.angles[,"Elevation"] == noon.max),] noon.min.date <- noon.angles[which(noon.angles[,"Elevation"] == noon.min),] # 3.Sun Reflection Height and Glare Intensity at Building ---- # mod.perp defines the relection/deflection boundary and is a line perpendicualr to the module. # the height calc relies on zenith angle only, ignoring azimuth angle impacts. The code is case # specific with a reflection sink to the North and utilizing only the elevation angles that will # reflect North. Udate angle.sort if deflection angles to the South are relevant. mod.perp <- 90 - mod.tilt angle.sort <- which(solI.angles[, "Elevation"] > 70 & solI.angles[, "Elevation"] <= 90) solI.angles <- transform(solI.angles, Bldg.Refl.Height = NA) solI.angles <- transform(solI.angles, Glare = NA) solI.angles[angle.sort, "Bldg.Refl.Height"] <- tan(D2R(mod.perp - (solI.angles[angle.sort, "Elevation"]-mod.perp))) * setback + mod.height solI.angles[angle.sort, "Glare"] <- solI.angles[angle.sort, "SolI"] * glare.coef # Plot Data ---- # color definition for sunlight hours ramp1 <- c("darkred", "red", "coral3", "orange", "yellow3", "gold", "yellow", "greenyellow", "green3", "cadetblue4", "navy", "darkorchid", "magenta") # subset data plot.sort <- which(hour(solI.angles[, 1]) >= 6 & hour(solI.angles[, 1]) <= 18) # plot elevation angle by hour ggplot(solI.angles[plot.sort, ], aes(TimeStamp, Elevation, group = Hour, color = Hour)) + geom_line(size = 1.2) + scale_color_manual(values = ramp1) + scale_y_continuous(limits = c(0, 90), breaks = seq(0, 90, 10)) + scale_x_datetime(breaks = date_breaks("1 month"), labels = date_format("%b")) + labs(title = expression(atop(bold("Solar Elevation Angles By Hour"), atop("Ras Laffan Industrial City", "Easting = 222580, Northing = 464288"))), x = "Month", y = "Solar Elevation Angle") # plot azimuth angle by hour ggplot(solI.angles[plot.sort, ], aes(TimeStamp, Azimuth, group = Hour, color = Hour)) + geom_line(size = 1.2) + scale_color_manual(values = ramp1) + scale_y_continuous(limits = c(-120, 120), breaks = seq(-120, 120, 20)) + scale_x_datetime(breaks = date_breaks("1 month"), labels = date_format("%b")) + labs(title = expression(atop(bold("Solar Azimuth Angles By Hour"), atop("Ras Laffan Industrial City", "Easting = 222580, Northing = 464288"))), x = "Month", y = "Solar Azimuth Angle") # plot reflected height at building by hour refl.sort <- which(solI.angles[, "Bldg.Refl.Height"] > 0) ggplot(solI.angles[refl.sort, ], aes(TimeStamp, Bldg.Refl.Height, group = Hour, color = Hour)) + geom_line(size = 1.2) + scale_color_manual(values = c("coral3", "orange", "yellow3")) + scale_y_continuous(limits = c(250, 500), breaks = seq(250, 500, 25)) + scale_x_datetime(breaks = date_breaks("1 month"), labels = date_format("%b")) + labs(title = expression(atop(bold("Solar Reflection Height By Hour"), atop("Ras Laffan Industrial City", "Estimated Height at Front of MPAC Building"))), x = "Month", y = "Reflection Height (meters above ground)") # plot glare intensity by hour glare.sort <- which(solI.angles[, "Glare"] > 0) ggplot(solI.angles[glare.sort, ], aes(TimeStamp, Glare, group = Hour, color = Hour)) + geom_line(size = 1.2) + scale_color_manual(values = c("coral3", "orange", "yellow3")) + scale_y_continuous(limits = c(11.5, 12.5), breaks = seq(11.5, 12.5, 0.1)) + scale_x_datetime(breaks = date_breaks("1 month"), labels = date_format("%b")) + labs(title = expression(atop(bold("Solar Glare By Hour"), atop("Ras Laffan Industrial City", "Estimated Irradiance Intensity"))), x = "Month", y = "Irradiance Intensity (W/m^2)") #. Internal Functions ---- # convert degrees to radians D2R <- function(x){ rad <- x * pi / 180 return(rad) } # convert radians to degrees R2D <- function(x){ deg <- x * 180 / pi return(deg) } # convert Easting and Northings coordinates (QND95) to decimal degrees (WGS84) # requires a dataframe with Easting and Northing columns QND2Dec <- function(coord){ coordinates(coord) <- c("Easting", "Northing") proj4string(coord) <- CRS("+proj=tmerc +lat_0=24.45 +lon_0=51.21666666666667 +k=0.99999 +x_0=200000 +y_0=300000 +ellps=intl +towgs84=-119.425,-303.659, -11.0006,1.1643,0.174458,1.09626,3.65706 +units=m +no_defs") xy.lonlat <- spTransform(coord, CRS("+proj=longlat +datum=WGS84")) dimnames(xy.lonlat@coords)[[2]] <- c("Longitude", "Latitude") dimnames(xy.lonlat@bbox)[[1]] <- c("Longitude", "Latitude") return(xy.lonlat) } # Convert Decimal Degree coordinates (WGS84) to Eastings and Nothings (QND95) # requires a dataframe with Longitude and Latitude columns Dec2QND <- function(coord){ coordinates(coord) <- c("Longitude", "Latitude") proj4string(coord) <- CRS("+proj=longlat +datum=WGS84") xy.eastnorth <- spTransform(coord, CRS("+proj=tmerc +lat_0=24.45 +lon_0=51.21666666666667 +k=0.99999 +x_0=200000 +y_0=300000 +ellps=intl +towgs84=-119.425,-303.659, -11.0006,1.1643,0.174458,1.09626,3.65706 +units=m +no_defs")) dimnames(xy.eastnorth@coords)[[2]] <- c("Easting", "Northing") dimnames(xy.eastnorth@bbox)[[1]] <- c("Easting", "Northing") return(xy.eastnorth) } # Calc distance between two geographic points in decimal degrees # Uses the method "great circle distances" and assumes a perfect sphere # Relies on personal degree-to-radian (D2R) and radian-to-degree (R2D) conversion functions # South and west coordinates are assumed to be negative geodist <- function(lat1, lon1, lat2, lon2) { gd <- (sin(D2R(lat1)) * sin(D2R(lat2))) + (cos(D2R(lat1)) * cos(D2R(lat2)) * cos(D2R(abs(lon2-lon1)))) gd <- R2D(acos(gd)) return(111.23 * gd) } |