4

My goal is to produce a binary map indicating whether the moon is visible at a certain time of day with GRASS GIS. (i.e. A mountain could be blocking the view of the moon.) r.sunmask would be the perfect module. Pretend the moon is the sun and just provide the azimuth and altitude. However the module runs too slowly. r.sun is then the recommended substitute. However r.sun does not provide for direct azimuth/altitude input, just the time of the day. I could somehow do a funky backwards calculation of the time when the sun would match the azimuth/altitude of the moon, but I'm not sure that would work 100% of the time.

What is a good way to accomplish this task?

p.s. Is web Mercator (EPSG:3857) suitable for this task or do I need a different projection?

2 Answers 2

2

I have no idea how to solve your problem with GRASS, but I am quite confident that it shouldn't be too complicated to do what you want with R (provided you know a little about R ;) )

Just made up a little example. Note that I have no idea about the accuracy or other caveates hidden somewhere.

library(oce)
t <- ISOdatetime(2016, 10, 22, 3, 00, 00, tz = "UTC") - 1*3600 # one time zone difference to UTC
#the position coordinates
lat <- 49.5
lon <- 11
###
#moonAngle() from package oce gives you infos about the moon for a certain date and location.
###
ma <- moonAngle(t, lon, lat)

#extract the azimuth and altitude from the returned list
azimuth <- ma$azimuth #141
altitude <- ma$altitude #45.5

library(raster)
#download example DEM
alt <- getData('alt', country='DEU') #Example DEM
plot(alt)

library(insol)
#produce the shadow map with the DEM and the vector of the light beam
sv=normalvector(altitude,azimuth)
shadow <- doshade(alt, sv)

plot(shadow)
#save result as tiff
writeRaster(shadow, "Shadowmap", format = "GTiff", overwrite=TRUE)

library(sp)
#Check if your position is in the shadow
S <- SpatialPoints(list(lon,lat))
plot(S,add=T)
res <- extract(shadow,S)
#print the result: 0 = no shadow, 1 = shadow
res #0

https://cran.r-project.org/web/packages/oce/oce.pdf

https://cran.r-project.org/web/packages/insol/insol.pdf DEM Shadow map

3
  • This would work. The heart of the solution is the doshade() function from the insol package. This also brings up the larger question of what is the best tool for the job. The product is a web app that does several calculations: moonshade, viewshed, etc. I have no urge to reinvent the wheel, but the performance needs to be snappy. R, matlab, octave, etc. are possibilities but I could only find a viewshed extension for matlab (which I don't own.) I'll leave this question as unanswered, because I would rather not have to load data into another question. I also don't know how fast R would run.
    – Justin
    Oct 22, 2016 at 8:43
  • Viewshed (of whatever?) can be incorporated from GRASS into R via the rgrass7 package. loc <- initGRASS(gisBase="/usr/lib/grass70",home=tempdir(),override=TRUE ) # Import DEM to GRASS and set region execGRASS('r.in.gdal', flags=c('o',"overwrite"), input= "DEM_20k.tif", output='rastgrass') execGRASS("g.region", parameters=list(raster="rastgrass") )
    – Bernd V.
    Oct 22, 2016 at 12:21
  • execGRASS("r.viewshed", flags=c('c','r','b',"overwrite"), parameters=list(input="rastgrass", output=output, coordinates=coords, observer_elevation=elevation, max_distance=10000)) R itself is fast as ligthning, the GRASS part takes longer. But for a tool like you have in mind, you probabaly need other stuff. Am no coder, just have an idea how to solve everyday problems ;)
    – Bernd V.
    Oct 22, 2016 at 12:22
2

If you have your data on your webserver, you can have/keep it in a GRASS database in GRASS format. Then you can skip the import part (r.in.gdal), which takes some time depending on the data size.

R is fast because data is handled in memory. If you have enough memeory available and if you want to speedup GRASS you can write the temoprary maps only to memory (r.external.out format=mem) and not to disk (which saves you some more time).

Might be worth having a look at WPS and the ZOO project: http://zoo-project.org/trac/wiki/ZooWebSite/QGIS_WPS_Client http://zoo-project.org/site/

P.S.: You can use Web Mercator (EPSG:3857 and the like) for visualization, but never use it for analysis! See: http://earth-info.nga.mil/GandG/wgs84/web_mercator/

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.