I have been using AMET v1.4 to compare MPAS “global” model results against observations without any problem. However, when I try to compare MPAS “limited area” model results, then I get an error:
“Mapping MADIS obs sites to MPAS grid
Error in cwgt[sitenum, 1] <- (XP * (YB - YC) + YP * (XC - XB) + XB * YC - :
replacement has length zero”
I tried to track down the error and found out that the eror occurs in function “mpas_site_map”. However, I cannot locate this function (I assume it is somewhere in one of the R scripts, but cannot find it) and hence cannot continue debugging. First I thought that limited area domain must have to do with this, but when I compare WRF “limited area” model results against observations, I don’t have any problem, either. So, I am not sure what is happening. Output seems to suggests that, based on some location finding algorithm, the analysis script is trying to co-locate the site on the model domain. But without locating this function, I cannot move forward with this error. Can you help me to solve this error?
I cut and pasted the full output below for your review. Thank you for your help.
With regards,
Kemal.
-------------- beginning of output ------------------
Date/Time START
Mon Nov 25 11:14:12 PST 2019
Loading required package: RMySQL
Loading required package: DBI
Loading required package: date
Loading required package: ncdf4
Read Model Variables
passed 1
Matching MPAS output file with observations: /share/mmb01/ST/MPAS/2016/limited.area/output.files/history.2016-07-18_00.00.00.nc
passed 2
Creating new database if it does not exist: CREATE DATABASE IF NOT EXISTS amet_test ;
Database Action: An AMET ** stations ** table exists in database amet_test . Will use this existing table.
Database Action: An AMET ** project_log ** table exists in database amet_test . Will not re-create this existing table.
Database Action: mpas_2016_test_surface AMET table exists. Will add new model-obs pairs to this existing table.
Reading MPAS Model Output File: /share/mmb01/ST/MPAS/2016/limited.area/output.files/history.2016-07-18_00.00.00.nc
Read latCell
Read lonCell
Read latVertex
Read lonVertex
Read time
Read t2m
Read q2
Read u10
Read v10
Read swdnb
Read surface_pressure
passed 3
MPAS time period skipped because initial model time
passed 4
Opening MADIS metar for time: 20160718 00:00:00
/aa/kgurer/AMET-1.4/obs/MET/point/metar/netcdf/20160718_0000
passed 5
Mapping MADIS obs sites to MPAS grid
Error in cwgt[sitenum, 1] <- (XP * (YB - YC) + YP * (XC - XB) + XB * YC - :
replacement has length zero
Calls: mpas_site_map
Execution halted
Date/Time END
Mon Nov 25 11:14:14 PST 2019
--------------- end of output ----------------------
Thank you for your kind reply. I appreciate that. It is great to hear that there is a work-around for this problem, which I will certainly try. Having said that, we still need to find out if there is a bug in AMET code. I had not tried with AMET 1.3 version before when I was writing my question, but then I tried yesterday with the previous version and sure enough version 1.3 also gives the same error. Since I don’t have problem with WRF limited area model comparisons with observations and I can plot 3 km resolution MPAS limited area model results without any problem, I am guessing there is a small error in AMET code in finding the location of observation point over the limited area MPAS domain. SInce I cannot find the “mpas_site_map” function, I got stuck. Just to clarify, I change “only” the path of model output from global MPAS model outputs to the path for the limited area domain MPAS model output (nothing else changes in AMET script) and then the error comes up.
That function is here: $AMETBASE/R_db_code/MET_site.mapping.R
With that said, I suspect the error is related to the limited area MPAS. At the time of development there was no such option in MPAS. The function is based on the full global MPAS. We have no limited area domains, so would be greatful for a one day MPAS output for testing. Without even looking at the code, I suspect this is occurring when and observation is outside of the limited area domain and a calculation results in zero or NA. It’s likely a 10 min fix that involves catching the NA’s and skipping to the next side like the WRF side of things when an ob is outside of the domain. As it stands now, the assumptions is that all sides reside in the global domain.
Going to ping you via email for an FTP to get a sample output.
I’m trying to install AMET1.4. I’ve installed all other required software. I installed R (3.6.1) but couldn’t find the package stats in it. I only found stats19, statsr, statsguru, statsExpressions, STAT and others. Which one should I install for the AMET to work properly? Also, the wgrib software couldn’t be downloaded successfully on the link you provided. Is there an alternative to it?
Are you trying to find the package stat or stats? Unless I’m mistaken, AMET requires the package stats and not stat.
Also, did you attempt to install the stat(s) R library separately from your base R install? Some of the libraries needed by AMET will need to be downloaded separately from the base R install.
Yes, that indicates that the stats package is installed and I think was actually already loaded and available in your R session (as R did not respond that it was loading the required package as it would with other packages). Since it’s in the base R installation, R likely pre-loads that package when you start an R session.
Hopefully you won’t run into any further issues with R.
I must be blind. We have a label for it: “looker-blind”. One looks at it but cannot see it, like you open the fridge, the food that you are looking for is in front of you but you cannot see it. I searched the function in all R scripts and F files. I either mistyped the search script, so I couldn’t find; or it found it but I didn’t see it. In any case, thank you for pointing to the script where the function is.
I solved the problem. As I suspected, when a station lat/lon is outside of the model domain, closest MPAS vertex to each station becomes null (cind[sitenum,2] and cind[sitenum,3]) and hence XB, XC, YB, YC points become null and the calculation stops with zero length in cwgt. I calculated minlat, maxlat, minlon and maxlon before station do-loop and put a multiple conditional if statement within the station do-loop and allowed the script to continue running only if station lat/lon is within the domain, then it worked. Here are the lines that I added:
Within the do-loop, right after site index “inds”:
if(slat[inds]>=minlat && slat[inds]<=maxlat && slon[inds]>=minlon && slon[inds]<=maxlon){ next }
AMET is working now. However, since the script is determining whether that station is within the domain and number of lat/lon vertices are quite large, script is taking a long time (still short enough) to move forward. If we can create a list of acceptable stations from the stations falling into the model domain before the do-loop, we can totally avoid this determination and script will run faster.
Thank you very much for the solution and pointing me the right script. I still cannot believe that I missed it. I apologize taking your time.
Well done sir. I will get a fix in our official archive asap.
Need to apologize about the latency in my response. Wyat and I just found that the CMAS forum emails are being put in my clutter directory for some strange reason.
And I messed up and did not send you the password to the FTP. If you could still upload a sample output, it would be appreciated. Will send you the password in a private email.
No problem. Just send me an email to my work email with the login credentials and I will upload the files. In the meantime, I sent you an email directly with
the files that I provided using our ftp server. If you cannot download the files from our ftp site, I will upload the files to scienceftp2.epa.gov. Please let me know.
The wgrib software couldn’t be downloaded successfully on the link you provided. Is there an alternative to it?
Also, 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)
