rxtracto_3D() axial extraction error

[makratofil]

Apologies for bringing up yet another issue! I'm not sure if you'll be able to reproduce this given the certificate issues you've been having with the server I'm extracting data from, but wanted to check and see if you've came across this before:

I'm essentially doing the same extraction in the other issue except for a bounding box for several tags off the coast of California. I'm using the apdrc server as it has the global HYCOM dataset at the spatial and temporal resolution that I need, and have been having a hard time finding a similar dataset on other ERDDAP servers.

sshInfo <- rerddap::info('hawaii_soest_f75b_adc6_12ab', url = 'https://apdrc.soest.hawaii.edu/erddap/') #SSH
sshInfo

# set variable information
xpos <- c(-119.911, -117.6807)
ypos <- c(32.81567, 34.2595)
tpos <- c("2012-01-20", "2012-01-20")
zpos <- c(0., 0.) # set altitude to zero 

ssh3D <- rxtracto_3D(sshInfo, parameter = "ssh", xcoord = xpos, ycoord = ypos, tcoord = tpos)
plotBBox(ssh3D, plotColor = 'algae')

The initial warning I got:

Warning message:
  In CPL_crs_from_input(x) :
  GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.

And the resulting plot looks like this:

I suppose the strangest thing is that I didn't encounter this kind of issue with other tags up until this point, and the previous bounding boxes/tags I extracted data for were in the same region and came from the same parent dataset.

[rmendels]

@makratofil the entire r-spatial group of packages have undergone rapid changes that also depend on newer versions of gdal, geos and proj4. It has major caused headaches for everyone, breaking a lot of packages. I tested that plotdap was still working on the vignette a week or two ago, but it may not be working in this case. Or it may be incompatibilities between versions you have of sf, gdal etc etc. Or I may have different versions of the underlying libraries. Will let you know.

[rmendels]

@makratofil and yes the certificate problem hasn't been fixed yet on that server, so I can not access it.

[rmendels]

@makratofil any chance you could save the extract ssh3D to an .RData file and either post it or send it to me? For plotdap, that is all I need.

[makratofil]

@rmendels Yep- I'll email the file to you directly!

[rmendels]

@makratofil First off, the warning message is now something that rgdal throws because of the incompatibilities introduced in new versions of proj4 wanting projection information as WKT2 strings, rather than the older crs strings. It is just a warning however. Also the 'pplotBBox()' plot is correct, from your code (rather than what you sent me) try the following:

ssh <- drop(ssh3D$ssh)
image(ssh)

and you should get

It is possible the extract is working incorrectly, but I can't test that at the moment.

[makratofil]

@rmendels Great, thank you! Please let me know if you do find some issue with the extract- but good to know that the plotting function works!

[rmendels]

@makratofil Okay I went to the ERDDAP web page and downloaded a netcdf file of the data for a slightly larger area, and below is the image. Note you get the same kind of patchiness, and also that for the more restricted area you extracted it looks correct.

So I may close this issue.

[makratofil]

@rmendels Thanks for doing that! I haven't had much time to revisit this recently but regarding the abrupt edge/line in my plot it may be due to the rgdal/proj4 versions I currently have installed- I'll try updating them and re-running this and let you know if that works. Thanks!

[makratofil]

@rmendels I've just been able to revisit this issue and after updating my r-spatial packages the error (edge/line through the plot) still occurred.

Looking further into the issue I believe it is associated the xcoord specification. The SSH dataset longitude range is from around 74-434 deg east and I previously specified xcoord as c(-120, -118). Making xcoord c(242, 240) generated a correct extraction (no abrupt edge/line through plot). The documentation for the package says that either form may be specified for extraction (-180-180 or 0-360). Are the make360 and make180 internal functions only applied for datasets that cross the dateline, or for all datasets? I.e., does the extraction function internally convert the longitude coords in accordance to the erddap dataset coords or is this something to be done manually, as I did? Just curious- I've been able to successfully extract using the 0-360 xcoord specification thus far so can close this issue.

Thank you for taking the time to look at this!

[rmendels]

@makratofil that suggests the problem lies in the fact that longitudes go beyond 360, which is non-standard. Don't know if the problem is what sent as the bounds for the extract, how those bounds are mapped back in output, or when the extract is mapped. I will need to get in and check.

[rmendels]

@makratofil I take it back, I still can't extract from that site because of certificates. Can you do me a favor. Run the call both ways, into different variables, but include the option "verbose = TRUE" in the extract. Please send me the verbose output from the two calls, as well as the the results of the two extracts. You may have sent some of this previously, but I misplace things easily.

[makratofil]

@rmendels Much apologies, have yet another clarification in the source of the issue! Played around with it a bit and the extraction output is fine if I specify xcoord = c(-120, -118) or c(240, 242) and the coordinates listed in the verbose outputs are the same between the two, suggesting that the issue isn't with the -180-180 or 0-360 formatting (although the non-standard scale may still be at play).

For the analysis I'm doing, I have a standard way of defining the bounding box or surrounding region to extract data across all tags, so I'm not specifying a ball park estimate of bounding coordinates. The extraction goes wrong when I use the specific bounding x coordinates: c(-119.911, -117.6807) or c(240.089, 242.3193). I think I pointed this out previously, but I've only had this issue with particular tags/bounding sets that are very similar regions/bounds of other sets that I've been able to successfully extract. I've confirmed the successful extractions by graphing the data on the server/or manually downloading and graphing myself, as you did above.

I've provided the code and verbose outputs below- the verbose outputs didn't come up on the console immediately following execution of the call, if that is normal? Let me know if this is the information you're looking for. And thanks much again!

For all extractions, these variables remain the same:

ypos = c(32.81567, 34.2595)
tpos = c("2012-01-20", "2012-01-20")
sshInfo <- rerddap::info('hawaii_soest_f75b_adc6_12ab', url = 'https://apdrc.soest.hawaii.edu/erddap/') #SSH
sshInfo

<ERDDAP info> hawaii_soest_f75b_adc6_12ab 
 Base URL: https://apdrc.soest.hawaii.edu/erddap/ 
 Dimensions (range):  
     time: (2008-12-28T00:00:00Z, 2013-08-20T00:00:00Z) 
     latitude: (-78.64, 46.8781) 
     longitude: (74.16003, 434.08003) 
 Variables:  
     emp: 
     mld: 
         Units: m 
     mlp: 
         Units: m 
     qtot: 
     s_trend: 
         Units: PSU 
     ssh: 
         Units: m 
     t_trend: 

The extraction piece:

ssh3D <- rxtracto_3D(sshInfo, parameter = "ssh", xcoord = xpos, ycoord = ypos, tcoord = tpos, verbose = T)

The verbose output for xcoord = c(-120, -118)

 * Found bundle for host apdrc.soest.hawaii.edu: 0x26bd1520 [can pipeline]
        * Re-using existing connection! (#14) with host apdrc.soest.hawaii.edu
          * Connected to apdrc.soest.hawaii.edu (128.171.151.240) port 443 (#14)
            > GET /erddap/griddap/hawaii_soest_f75b_adc6_12ab.nc?ssh[(2012-01-20T00:00:00Z):1:(2012-01-20T00:00:00Z)][(32.8387):1:(34.239)][(118.00003):1:(120.00003)] HTTP/1.1
            Host: apdrc.soest.hawaii.edu
            User-Agent: libcurl/7.64.1 r-curl/4.3 crul/0.9.0
            Accept-Encoding: gzip, deflate
            Accept: application/json, text/xml, application/xml, */*

              < HTTP/1.1 200 
            < Date: Thu, 23 Jul 2020 19:52:36 GMT
            < Server: Apache/2.4.6 (CentOS) OpenSSL/1.0.2k-fips PHP/7.3.9
            < Last-Modified: Thu, 23 Jul 2020 19:52:36 GMT
            < xdods-server: dods/3.7
            < erddap-server: 1.80
            < Content-Disposition: attachment;filename=hawaii_soest_f75b_adc6_12ab_522a_2704_c136.nc
            < Content-Encoding: gzip
            < Content-Type: application/x-netcdf
            < Transfer-Encoding: chunked
            < 
              * Connection #14 to host apdrc.soest.hawaii.edu left intact

And the verbose output for xcoord = c(240,242)

* Found bundle for host apdrc.soest.hawaii.edu: 0x26bd56a0 [can pipeline]
    * Re-using existing connection! (#13) with host apdrc.soest.hawaii.edu
      * Connected to apdrc.soest.hawaii.edu (128.171.151.240) port 443 (#13)
        > GET /erddap/griddap/hawaii_soest_f75b_adc6_12ab.nc?ssh[(2012-01-20T00:00:00Z):1:(2012-01-20T00:00:00Z)][(32.8387):1:(34.239)][(240.00003):1:(242.00003)] HTTP/1.1
        Host: apdrc.soest.hawaii.edu
        User-Agent: libcurl/7.64.1 r-curl/4.3 crul/0.9.0
        Accept-Encoding: gzip, deflate
        Accept: application/json, text/xml, application/xml, */*

          < HTTP/1.1 200 
        < Date: Thu, 23 Jul 2020 17:59:17 GMT
        < Server: Apache/2.4.6 (CentOS) OpenSSL/1.0.2k-fips PHP/7.3.9
        < Last-Modified: Thu, 23 Jul 2020 17:59:17 GMT
        < xdods-server: dods/3.7
        < erddap-server: 1.80
        < Content-Disposition: attachment;filename=hawaii_soest_f75b_adc6_12ab_e31f_83d1_877c.nc
        < Content-Encoding: gzip
        < Content-Type: application/x-netcdf
        < Transfer-Encoding: chunked
        < 
          * Connection #13 to host apdrc.soest.hawaii.edu left intact

And finally, the verbose output for the extraction that generates the abrupt axial edge in the plot:

* Found bundle for host apdrc.soest.hawaii.edu: 0x26bd54e0 [can pipeline]
* Re-using existing connection! (#4) with host apdrc.soest.hawaii.edu
  * Connected to apdrc.soest.hawaii.edu (128.171.151.240) port 443 (#4)
    > GET /erddap/griddap/hawaii_soest_f75b_adc6_12ab.nc?ssh[(2012-01-20T00:00:00Z):1:(2012-01-20T00:00:00Z)][(32.8387):1:(34.239)][(240.08003):1:(242.32003)] HTTP/1.1
    Host: apdrc.soest.hawaii.edu
    User-Agent: libcurl/7.64.1 r-curl/4.3 crul/0.9.0
    Accept-Encoding: gzip, deflate
    Accept: application/json, text/xml, application/xml, */*

      < HTTP/1.1 200 
    < Date: Thu, 23 Jul 2020 16:54:39 GMT
    < Server: Apache/2.4.6 (CentOS) OpenSSL/1.0.2k-fips PHP/7.3.9
    < Last-Modified: Thu, 23 Jul 2020 16:54:39 GMT
    < xdods-server: dods/3.7
    < erddap-server: 1.80
    < Content-Disposition: attachment;filename=hawaii_soest_f75b_adc6_12ab_5688_16d1_071e.nc
    < Content-Encoding: gzip
    < Content-Type: application/x-netcdf
    < Transfer-Encoding: chunked
    < 
       * Connection #12 to host apdrc.soest.hawaii.edu left intact

[rmendels]

@makratofil Thanks. I wish I could access that server so I could get into debug mode. The first one ( xcoord = c(-120, -118)) has produced an incorrect URL. My best guess is still that it has to do with the underlying longitudes not being on either (-180, 180) nor (0, 360). The reason what you may see may not be exactly what you ask for, is that I find the closest coordinates actually in the dataset and ask for those. This is described at the end of the Vignette, under something like "What happens when you make a request". Until I can access the server, it is going to very difficult to see what is going on.