dankelley / oce

R package for oceanographic processing
http://dankelley.github.io/oce/
GNU General Public License v3.0
145 stars 42 forks source link

Orientation not being assigned correctly from raw ADCP files #1588

Closed hhourston closed 5 years ago

hhourston commented 5 years ago

I'm working on some upward and downward facing moored ADCP data, collected using RDI Workhorse model instruments. For some raw ADCP files, the orientation is being read in wrong ("upward" instead of "downward" or vice-versa) by the function read.adp().

I am by no means an R expert. I am wondering if this issue is arising because the orientation is being assigned as the recorded orientation for the first ensemble, at which point the ADCP may be in the process of being deployed and is oriented with the wrong end up at that time. Is there a way to assign the orientation as the average orientation over all the ensembles?

dankelley commented 5 years ago

You can alter anything within an oce-created object, but things are pretty tricky when it comes to orientation, coordinate systems, and the transformation matrix. The firmware of your instrument may also matter, in terms of how oce decodes the data.

@richardsc is more familiar with these issues than I am. I think he may be out of contact for a while, but perhaps he'll see this issue and add something helpful.

In a lot of cases, the first thing we do is to ask the reporter to make a sample of the data available to us privately. I mention privacy because Clark and I both know what it is like to need to hold data private for a while, e.g. before publishing or completing a thesis. With ADCP files, the "ensembles" are often just 1 to 10kb long, so using the unix command split is a good way to chop out say 200kb that can be sent in an email to us. (You can get my email at my profile page, and I can forward a file to Clark.)

If I were you, I'd split out a sample file and send it, so Clark and I can have a look. But it might take a week or more before Clark will be in contact, so you might also want to just figure out things by yourself, if you are good at geometry. For example, the vertical velocity might just be of the wrong sign ... but the other two components are a lot trickier because of subtle issues like just what the instrument reports for a pitch when the machine is upside down, etc.

hhourston commented 5 years ago

Thank you for your suggestions! I've sent a sample file your way. I suspected it might be tricky... I'll keep poking away at it though.

richardsc commented 5 years ago

I have definitely seen this issue before and will take a look at the data file if @dankelley can forward to me. It likely won’t be until at least next week sometime though.

richardsc commented 5 years ago

Maybe related to #256 and #275

hhourston commented 5 years ago

Thanks for directing me to those!

richardsc commented 5 years ago

I didn’t look close enough to see if the solution is in there, but at least posted the links so I could find them later :smile:

dankelley commented 5 years ago

I'm going to use this spot to take some notes on the data format. I'll do it by snapshots of the RDI manual, with one comment per snapshot. The snapshots will show the document name at the top, so I don't have to try to figure out how to refer to the doc.

dankelley commented 5 years ago
Screen Shot 2019-08-10 at 8 12 30 am
dankelley commented 5 years ago
Screen Shot 2019-08-10 at 8 09 16 am
dankelley commented 5 years ago
Screen Shot 2019-08-10 at 8 06 11 am
dankelley commented 5 years ago

From the above, it looks as though I can insert something in the code to examine the FLD, to find the appropriate bit (the one at the bottom two lines of the table in the just-previous comment ... finding a way to state the location of that bit in the file is a challenge, with LSB vs MSB, endianness etc.

Given that a FLD starts with 0x00 followed by 0x00, the proper spot to interpret this might be at https://github.com/dankelley/oce/blob/41ab9a456ac61dbf8739c48f7c5cebac77ea820c/R/adp.rdi.R#L1216 in the code. However, I am not 100% sure on that because the code complexity is high (lots of if blocks), owing to the need to handle a wide variety of instrument types.

dankelley commented 5 years ago

NOTE: There is 1 question for @hhourston at the end and 2 questions @richardsc at the end. (Sorry this comment is so long.)

I have instrumented the code to check for changes in the two bytes that are used for this. You can get this version at git commit 32935158991b37b3b20c2e11a18f9bc3e13583dc of the "develop" branch.

What I'm finding is that the first "configuration", stored in header$systemConfiguration (as inferred from the first profile/ensemble in the data file) is

11001010-01000001

and these continue up until profile/ensemble number 27, at which point we get

01001010-01000001

In this notation -- at least as oce is working now -- the first digit in these character sequences means the orientation, with 0 meaning downward and 1 meaning upward.

At a quick glance, it looks as though all the profiles/ensembles after 27 have a 0 in the first character.

Q1 for @hhourston -- As you can see in https://github.com/dankelley/oce/issues/1588#issuecomment-520140680 above, the 8 bits examined by read.adp.rdi() also contain info on instrument frequency and beam convex/concave setup (plus some things I don't know about). When I read the file with d<-read.ETC and then do summary(d), I see a freq of 300kHz (inferred because the right-most 3 bits of the first byte are 010). Is that correct? Also, when I do d[["beamPattern"]] I get "convex". Is that right? Knowing whether these two things are right will give me a bit more confidence that I am not looking at the bits wrong.

Q1 for @richardsc -- Have you come across this sort of thing in his previous work with RDI, with the initial orientation as inferred by read.adp.rdi() not matching the known orientation in the water? Maybe e.g. in SLEIWEX, adps were set up always in "pods" so the on-deck orientation matched the in-water orientation.

UPDATE: Q2 is moot; see next comment

Q2 for @richardsc -- Sorry, Clark, but just thinking out loud here ... I have to admit that I am perplexed by the RDI docs. They say to take two bytes (numbers 5 and 6) and "convert them" to binary (how??) and then they give a table with 8 bits, not 16 bits. HUH??? Why does the table not list 16 bits? And, how exactly are we to convert ... to a 2-byte int, a 2-byte unsigned int, etc? If it's the latter, it means just stringing the bits together, as I am doing presently in the code near https://github.com/dankelley/oce/blob/32935158991b37b3b20c2e11a18f9bc3e13583dc/R/adp.rdi.R#L92 I know I've been confused about this before. Maybe RDI is thinking that it's impossible to look within the bits of a one-byte entity (is that something for C on msdos/windows machines?). I guess the table is referring just to the first byte in this sequence.

dankelley commented 5 years ago

Oh, hang on. The RDI docs actually show both bytes ... I had only captured one of them in a comment above. The full table is as below. Therefore, question 2 for @richardsc in the previous comment is now moot ... sorry!

Screen Shot 2019-08-10 at 10 07 59 am
dankelley commented 5 years ago

I've updated oce ("develop" branch commit 02b1048ea7539abb7d829196a4aa8a6f0a3e3210) to handle multiple values of @metadata$orientation.

Here's what works so far:

Here's what does not work so far:

Since the test data file is not public (well, I forgot to ask @hhourston whether the data are secret, so I am assuming so) I will not post test code here, nor sample graphs. However, I will email these things to @hhourston and @richardsc in a few minutes.

hhourston commented 5 years ago

@dankelley Both the frequency and beamPattern are correct for that particular file. Thanks for sharing the data format notes!

richardsc commented 5 years ago

Hi @hhourston! Sorry for the slow follow-up on this -- I was on vacation and then swamped with all the stuff waiting on my desk when I got back.

I have encountered this issue generally before, though I'm not sure if I've encountered it specifically for RDI ADCPs. There is definitely a caution here about how to handle such data, and whether it should be oce or the user (who knows their data and deployment) who should account for this.

I'm rebuilding oce now, and will take a look at the data file myself, but first, a few general points (basically, just my Monday-morning brain thinking aloud).

  1. The fact that the ADCP records a time series of the orientation is handy, but it presents some problems for how to determine the correct "deployed" orientation for doing things like beam-to-ENU transformations. Clearly simply reading the first one and using that to assign a @metadata value is a bit dangerous.

  2. As @dankelley mentioned, if your instrument is recording in ENU (i.e. the ADCP is doing the transformations on a ping-by-ping basis, then you shouldn't have to worry about the velocities -- the orientation is taken into account by the ADCP firmware. So, just trim your file to the "in-water" times, and you can safely ignore the orientation flag.

  3. When oce reads a raw file, it sets the orientation flag based on the first one in the file, and any subsequent transformations are based off of that. That will be problematic in the case of an instrument that recorded in "beam" coordinates (or potentially XYZ/instrument coordinates) because the wrong formula will get used (see the RDI coordinate transformation manual).

  4. We could change oce so that it includes the time series of orientation in the coordinate transformation (i.e. every ensemble uses a value for orientation to ensure the correct formula is applied). This will probably slow the calculation down, but I think it's done is C already so probably isn't a huge overhead.

  5. It might make sense for oce to detect and issue a warning() when it reads a file for which the orientation is not constant. That way the user has a heads-up that they should be careful.

  6. In such instances, my approach in the past has been to simply trim the data to the in-water times, manually set the orientation to what I know it's supposed to be (e.g. adp <- oceSetMetadata(adp, 'orientation', 'upward')) and then proceed with the toEnu() coordinate transformation.

  7. For some reason I can't open the emailed file with the newest version of oce ... I get:

    > d <- read.adp.rdi('a1_20050503_20050504_0221m.000')
    Error in decodeHeaderRDI(buf, debug = debug - 1) : 
      first two bytes in file must be 0x7f 0x7f, but they are 0x00 0x00

    EDIT: I figured it out. It was a PEBKAC situation

  8. However, I see that the example adp dataset has been updated, as I see:

    > str(adp@metadata)
     $ oceCoordinate       : chr "enu"
     $ orientation         : chr [1:25] "upward" "upward" "upward" "upward" ...
  9. The question I have is how does having a character vector or orientations affect any downstream processing such as beamToXyz(), toEnu(), xyzToEnu() etc ... ?

richardsc commented 5 years ago

It occurred to me that, at least for the simplest question that @hhourston asked:

Is there a way to assign the orientation as the average orientation over all the ensembles?

The answer is (buried in my previous post):

adp <- oceSetMetadata(adp, 'orientation', 'upward')

So, if you know the orientation, you can just set it directly.

Though, now that Dan has made some changes to how that is read in as a vector, we might have to think about downstream consequences.

dankelley commented 5 years ago

There are definitely going to be consequences but they do not fit the title of this issue, so I am going to make a new issue. (The point of that is to make it easier for users to find discussion based on searching through issue titles.)

I think this issue (as defined by its title and first comment) has been addressed, actually.

dankelley commented 5 years ago

This issue has been closed, since it seems to have been addressed, and since discussion has stalled. The original reporter should feel free to reopen it, though, or to open new issues that might be related. Thanks.

PS. This is a standardized reply.