Dear @cjcoats,
I have been able to download and install wgrib.
Thanks and regards,
Catalyst
Hello Robert,
I am requesting a better solution than what I concocted since the script is taking a long time, due to the conditional if statement that it has to process for each station. As I mentioned earlier, the values of XB, YB, XC, YC before the do-loop are null for limited area domain while they are not for global domain runs, which doesn’t make sense since they are the values of unique station lat/lon values. I had given you one-day of MPAS limited area domain output. I was wondering if you or your colleagues had a chance to process those output files to understand and confirm the problem that I am stating. The processing time is very very long in the current condition. I appreciate any feedback.
With regards,
Kemal.
Dear @kgurer , @wyat.appel @rgilliam @cjcoats,
I want your guide on amet-config.R script , I don’t know if I am correct like this or how should the directory paths be:
amet_base <- Sys.getenv(’/home/catalyst/Desktop/Build_WRF/LIBRARIES/AMET-master’)
mysql_server <- “localhost” ## Name of MYSQL server
amet_login <- “ametsecure” ## AMET Root login for MYSQL serve
amet_pass <- “1234526” ## AMET Root password for MYSQL server
maxrec <- -1 ## Set MySQL maximum records for queries (-1 for no maximum)
Bldoverlay_exe_config <- paste(amet_base,"/home/catalyst/Desktop/Build_WRF/LIBRARIES/AMET-master/tools_src/bldoverlay/src/bldoverlay_v53.exe",sep="")
EXEC_sitex_daily_config <- paste(amet_base,"/home/catalyst/Desktop/Build_WRF/LIBRARIES/AMET-master/tools_src/sitecmp_dailyo3/sitecmp_dailyo3_v53.exe",sep="")
EXEC_sitex_config <- paste(amet_base,"/home/catalyst/Desktop/Build_WRF/LIBRARIES/AMET-master/tools_src/sitecmp/src/sitecmp_v53.exe",sep="")
Thanks and regards,
Catalyst
Thanks
Catalyst
Hi Catalyst,
It looks like your AMET base directory is /home/catalyst/Desktop/Build_WRF/LIBRARIES/AMET-master. Is that correct?
If so, for amet_base, you can either hard-code the directory explicitly, like this:
amet_base <- “/home/catalyst/Desktop/Build_WRF/LIBRARIES/AMET-master”
or you can keep it as the default setting of
amet_base <- Sys.getenv(‘AMETBASE’)
which retrieves amet_base from the environment setting of AMETBASE in the other AMET scripts. Now, one user reported an issue retrieving amet_base from the environment variable AMETBASE, but I was not able to reproduce the error on my system. He solved the issue by hard-coding the base directory, which is the first option I gave above.
For the executable directories, right now it appears you will be putting the AMET base directory in the path twice, since you’re both pasting in amet_base and specifying amet_base explicitly. So, instead you can just do:
Bldoverlay_exe_config <- paste(amet_base,"/tools_src/bldoverlay/src/bldoverlay_v53.exe",sep="")
EXEC_sitex_daily_config <-paste(amet_base,"/tools_src/sitecmp_dailyo3/sitecmp_dailyo3_v53.exe",sep="")
EXEC_sitex_config <- paste(amet_base,"/tools_src/sitecmp/src/sitecmp_v53.exe",sep="")
That will specify the paths to your executables correctly. Let me know if any of that is confusing or if you run into any other issues.
Wyat
@rgilliam Were you able to get back to Kemel on this question on the one-day MPAS limited area domain output? It got sandwiched in between two other unrelated posts somehow.
Hello Robert,
I had originally submitted my solution to read and process limited area MPAS model outputs and had also mentioned that that process was very very slow. That solution was not great at all and I worked out a much better solution that works for limited area MPAS model results. Here is the changes that I made:
- I removed the changes that I had made to “MET_site.mapping.R” file
- I put the following lines into “MET_observation.read.R” file starting on line 504 right after nsr assignment:
----------------------------- Beginning of changes ----------------------------
#…
new_site = NULL
new_slat = NULL
new_slon = NULL
new_selev = NULL
new_site_locname = NULL
new_report_type = NULL
new_ihour = NULL
new_imin = NULL
new_isec = NULL
new_stime = NULL
new_stime2 = NULL
minlat = 31.
maxlat = 43.
minlon = -128.5
maxlon = -113.5
minlat = min(slat)
maxlat = max(slat)
minlon = min(slon)
maxlon = max(slon)
print(paste("minlat = ",minlat))
print(paste("maxlat = ",maxlat))
print(paste("minlon = ",minlon))
print(paste("maxlon = ",maxlon))
for(s in 1:nsr){
if(slat[s] >= minlat && slat[s] <= maxlat && slon[s] >= minlon && slon[s] <= maxlon ){
new_site[length(new_site)+1] <- site[s]
new_slat[length(new_slat)+1] <- slat[s]
new_slon[length(new_slon)+1] <- slon[s]
new_selev[length(new_selev)+1] <- selev[s]
new_site_locname[length(new_site_locname)+1] <- site_locname[s]
new_report_type[length(new_report_type)+1] <- report_type[s]
new_ihour[length(new_ihour)+1] <- ihour[s]
new_imin[length(new_imin)+1] <- imin[s]
new_isec[length(new_isec)+1] <- isec[s]
new_stime[length(new_stime)+1] <- stime[s]
new_stime2[length(new_stime2)+1] <- stime2[s]
print(paste("length of array ",length(new_site)));
writeLines(paste("Inside. Site # = ",s," ",site[s]," lat = ",slat[s]," lon = ",slon[s]))
}
else{
writeLines(paste("Outside. Site # = ",s," ",site[s]," lat = ",slat[s]," lon = ",slon[s]))
}
}
nnsr <-length(new_site)
#…
site <-list(ns=ns, nsr=nsr, site=site, sites_unique=sites_unique, slat=slat, slon=slon,
elev=selev, site_locname=site_locname, report_type=report_type, ihour=ihour,
imin=imin, isec=isec, stime=stime, stime2=stime2)
site <-list(ns=nnsr, nsr=nnsr, site=new_site, sites_unique=new_site, slat=new_slat, slon=new_slon,
elev=new_selev, site_locname=new_site_locname, report_type=new_report_type, ihour=new_ihour,
imin=new_imin, isec=new_isec, stime=new_stime, stime2=new_stime2)
------------------------------ End of changes ----------------------------
One issue that puzzles me is that slat/slon variables are read from MPAS output files at each model output time within MET_observation.read.R file. When min/max of lat/lon are calculated (as I tried to do in the above changes), it should have showed the min/max of limited area domain, but it doesn’t. I didn’t invest too much time on why it doesn’t show that and hard coded the min/max lat/lon values of my limited area domain fine. I have confirmed with print/write statements that only the observations that fall within the modeling domain are used. Finally, the “site” variable assignment located at line 505 of MET_observation.read.R file is replaced with the new one given above.
Kemal.
Hello Robert,
I had cut and pasted the changes that I made MET_observation.read.R file into this box. That process apparently separated the comment sign “#” in front of “minlat, maxlat, minlon and maxlon” variables as well as “site” variable that are seen as bold face characters now. Those bold faced lines/characters are all comments. Please make a note of it.
Kemal.
Hello Robert and all,
I need to let you all know that I had to make another modification to the solution that I had just introduced to read limited area domain MPAS model results. The latest modification to MET_observation.read.R file is as follows starting from line 504:
--------------------- Beginning of modification -------------------------
new_sites = NULL
minlat = 31.
maxlat = 43.
minlon = -128.5
maxlon = -113.5
for(s in 1:nsr){
if(slat[s] >= minlat && slat[s] <= maxlat && slon[s] >= minlon && slon[s] <= maxlon ){
new_sites[length(new_sites)+1] <- site[s]
print(paste("length of array ",length(new_sites)));
writeLines(paste("Inside. Site # = “,s,” “,site[s],” lat = “,slat[s],” lon = ",slon[s]))
}
else{
writeLines(paste("Outside. Site # = “,s,” “,site[s],” lat = “,slat[s],” lon = ",slon[s]))
}
}
nnsr <-length(new_sites)
site <-list(ns=ns, nsr=nnsr, site=site, sites_unique=new_sites, slat=slat,
slon=slon, elev=selev, site_locname=site_locname,
report_type=report_type, ihour=ihour, imin=imin, isec=isec,
stime=stime, stime2=stime2)
------------------- End of modification -----------------------------
For those of you who need to process limited area domain MPAS model results, please use your own judgement. The solution that I present may lnot be perfect, and I am not an expert in AMET (at all) and in R, but the latest solution that I just presented seems to work and fast.
Good luck in your own work.
Kemal.
I have to apologies. Looking back at this thread I had no idea of the activity here in the last few months.
Talked with Wyat before about this, but seems notifications are diverted to the Clutter box in our email system. Will look into this.
Hey, to close this request/comment out on the CMAS forum I have an even better fix that is quite simple. Will follow with an email with specific codes though, because I also fixed the surface met matching for cases where SW radiation or surface pressure are not in the output. But for the limited area MPAS, this is what I did.
- It is better to handle this in the site mapping function MET_site.mapping.R with a simple “skip” of a site that does not fall within a domain. You can see this is in place for the wrf site mapping function where site lat lon is compared against the min and max lat lon and if outside the domain, the next statement just skips that site.
Using the closest MPAS vertex to each MADIS station, store the cell indices for that vertex
cind[sitenum,1]<-cells_on_vertex[1,vmap]
cind[sitenum,2]<-cells_on_vertex[2,vmap]
cind[sitenum,3]<-cells_on_vertex[3,vmap]
# Skip site if outside of a limited area MPAS domain.
if( min(cind[sitenum,]) <= 0 || max(cind[sitenum,]) > length(latv) ) {
nsites_outside_domain<-nsites_outside_domain+1
writeLines(paste(nsites_outside_domain,"Sites ****",sites_unique[s],"**** excluded because out of MPAS domain"))
next
}Note the first three lines exist already. I added the part after comment about limited area domain.
Hi Robert,
Thank you for the clean solution that you introduced for reading MPAS limited domain output files. I will make the addition to “MET_site.mapping.R” file on my
machine.
With regards,
Kemal.
Finding the closest MPAS-vertex or the enclosing MPAS-cell is quite a non-trivial problem, and thus calls for a more-sophisticated algorithm. Perhaps you should adapt the algorithms from the I/O API’s MODULE MPASFIO (file ioapi/modmpasfio.f90) routines FINDCELL and FINDVRTX (see https://cjcoats.github.io/ioapi/MODMPASFIO.html#find)
Note that if the stations are close together, these routines take advantage of that fact to speed up their “steepest-descent” algorithm even further for the second and subsequent calls…
Good information. Not sure this was available back when we started looking at the issue.
Will take a look at the algorithm. With that said, there is a variable cellsOnVertex that gives the Vertex cell indices each cell, so no calculation is needed as long as you know the cell the observation site lies. That is found by finding the cell (center) that has the minimum distance from the obs site.
But “brute force and ignorance” leads to a very inefficient algorithm-- O(N) instead of O(log(N)) in terms of the grid’s cell-count N. And the actual on-microprocessor performance is even more affected, since the O(N) algorithm does not have good memory-system locality.