mapImage() with filledContour producing weird patterns

[mdupilka]

Starting about a months or so ago, using filledContour = TRUE in the mapImage() has been producing weird patterns. Here is an example

Here is the sample code to produce it

require(oce)
data("coastlineWorldMedium", package = "ocedata")
int <- .1
lon <- seq(-179.95, 179.95, int)
lat <- seq(-39.95, 89.95, int) 
m <- matrix(0, nrow = length(lon), ncol = length(lat))
mapPlot(coastlineWorldMedium,
        type = "polygon",
        col = "grey95",
        #projection="+proj=laea +lat0=40 +lat1=60 +lon_0=-110",
        projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110",
        #longitudelim = c(minlon, maxlon), latitudelim = c(minlat, maxlat),
        longitudelim = c(-140, -80), latitudelim = c(45, 70),
        polarCircle = 0,
        clip = TRUE, drawBox = TRUE,
        cex.axis = .7  )

mapImage(longitude = lon, latitude = lat, z = m,
         filledContour = TRUE,
         breaks = 10, 
         col = colorRampPalette(c("blue", "red"))(10) )

With filledContour = FALSE this does not happen

[dankelley]

I'll take a look. By the way, I edited your comment to put

```R

at the start of the code and

```

at the bottom. That makes it more readable. Also, I changed the title to be more specific. (It's good when users can search past issues by title.)

I just tried it, and I get the same thing. I've not looked at it in detail.

My first step will be to remove anything specific that likely does not matter, e.g. the polar circle, the cex.axis and so forth. It's helpful to have issues boiled down to the very smallest limit. So I'll likely also chop the domain down (since drawing that filled map-image is slow).

I'll likely look at this tomorrow.

Thanks!!

PS. the reprex package is very handy for reporting issues. It guarantees the developers that the code does as stated, and it gets around the fact that sometimes people make issues that require private data for us to check, which slows things down a lot.

[dankelley]

I think it's a moire effect. The image is actually drawn by mapping long-lat space into page space, finding the four-sided elements, and filling each one as a polygon. Polygons have straight sides in page space and thus curved lines in (many) projected spaces. I think the patterns result from a non-overlapping. Possibly before lines were being drawn at the borders of polygons, and those filled in the space. Anyway, I'll look into this tomorrow.

PS. notice that I boiled things down a bit in the reprex, and used a matrix that varies and that has a pattern that the eye can notice.

library(oce)
#> Loading required package: gsw
data("coastlineWorldMedium", package = "ocedata")
int <- 0.1
lon <- seq(-179.95, 179.95, int)
lat <- seq(-39.95, 89.95, int)
m <- matrix(seq_along(lat), byrow=TRUE, nrow = length(lon), ncol = length(lat))
mapPlot(coastlineWorldMedium,
    projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110",
    longitudelim = c(-140, -80), latitudelim = c(45, 70)
)

mapImage(
    longitude = lon, latitude = lat, z = m,
    filledContour = TRUE,
    breaks = 10,
    col = oceColorsJet
)

^{Created on 2024-03-02 with reprex v2.1.0}

[mdupilka]

Thanks for such fast responses. This only started happening about a month or so ago that I noticed.

[mdupilka]

My first thought is that is did look like some kind of interference pattern

[dankelley]

Hm, I don't even understand the polygons I see, because they do not correspond to lon-lat boxes or x-y-page boxes.

At some point in the code, we interpolate to find data. (I think for the default case this is done near https://github.com/dankelley/oce/blob/c8d72c1a8f380761015bec384802117de8ac3db8/R/map.R#L3630) and so what I'll do next (tomorrow, though) is to zoom the region/matrix some more so I can put some browser() statements in to see what is going on in those spots that are coming up white.

This may boil down to the business of connecting a cartesian space to a non-cartesian one. As you likely know, @mdupilka (since I saw you wrote C code), this comes up in computer graphics: the rule is to map from viewer space back to simulation space, not vice-versa, because otherwise you may get holes. But I find the shapes of the holes curious.

Thanks for the issue, and for posting it on a weekend when I have more time.

library(oce)
#> Loading required package: gsw
data("coastlineWorldMedium", package = "ocedata")
int <- 0.5
lon <- seq(-179.95, 179.95, int)
lat <- seq(-39.95, 89.95, int)
m <- matrix(seq_along(lat), byrow = TRUE, nrow = length(lon), ncol = length(lat))
mapPlot(coastlineWorldMedium,
    projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110",
    longitudelim = c(-140, -80), latitudelim = c(65, 70)
)

mapImage(
    longitude = lon, latitude = lat, z = m,
    filledContour = TRUE,
    breaks = 10,
    col = oceColorsJet
)
mapGrid(0.5, 0.5, col = "magenta")

^{Created on 2024-03-02 with reprex v2.1.0}

[richardsc]

Just spitballing, because I'm on my phone, but could this be related to any changes when we migrated from sp/rgdal etc to sf?

[dankelley]

It could be that transition, I suppose. But we don't have a hashcode for a previous version that worked, so it will be hard to know that. (Anyway, we cannot go back because sp and rgdal are either removed or will be very soon.)

I don't understand the shapes of the blank regions. They are not in the shape of long-lat patches. Nor are they in the shape of xy on page space. So that's surprising to me. Anyway, I'll know more because I have used mapLocator() to find a blank spot, and have also exported the results of the binning. So now I can do which.min(abs(DIFFERENCE)) tricks to isolate the spots of the gridded matrix that are NA (which is what I think white must be). I tried altering the filled and fillGap parameters of the gridder but the spaces remain. Kind of a mystery.

[dankelley]

The gridder is finding holes. The weird patch shapes are just happening when there is a hole in one corner of a square "pixel" on the grid, etc.

[dankelley]

Hm. Maybe the below (R/bin.R) explains why things changed in the past few months.

This also explains why my setting fill and fillGap had no effect in my tests -- those arguments are no longer used in bin_mean_2D(). I might look up my old C code and see if there is any reason not to re-insert it (beyond the fact that very few oce users know C, so it's hard for them to help to find bugs and submit patches). I think they changed something about the C compiler, and maybe some old-fashioned actions I was undertaken got flagged as bad.

I'll look into this tomorrow.

[dankelley]

Here is a commit where I chopped the C file https://github.com/dankelley/oce/commit/811ba902b6ddb1bd00e8ab22274a0ab8ff8bf3dd; prior to that, I had commented a lot of stuff out.

I think I had forgotten about fill and fillgap when I chopped the C code. (I can see why I chopped it: 350 lines of 1980s style C code. I'm not even using Rcpp (which links R objects to C++ objects) so it's pretty primitive and code like that is hard to debug and maintain.

[dankelley]

The snapshot is my final comment for the day. This is how I am isolating cells in which to perform the average. Some of the cells are NA when computing the mean value, because of a lack of data therein.

[dankelley]

I've decided to go back to my C way of doing the 2D gridding, because it has a fillGap parameter that I must have written problems of this sort. Since the change is only a few months old, and since I should be able to do this without causing problems, I will go ahead and do this without asking for a co-developer vote ... although any watching this ought to feel free to object.

For me, the arguments in favour are:

1. It will solve this particular bug. (I think this bug might come up a fair bit.)
2. It will return to previous (very old) behaviour for binMean2D(), perhaps avoiding problems that had not yet been reported.
3. I don't think the move is too risky, because that C code is very old and likely was not causing problems. (I switched to using split() and cut() because I thought the code would be easier for others to maintain. But we lost the fillGap possibility with that change.)
4. I do not intend to do this for the 1D case, just the 2D case. That means I will be exposing only half the old C code. (The less C code, the better. Oh, and I imagine that cut() and split() are are calling C for speed.)

Note that I will do this in a new branch called issue2199.

[dankelley]

I've re-inserted the old C code, and the bad gaps are gone now. I am not 100% convinced that the gaps are being filled correctly, though. I'm attaching two images. The first is a snapshot of a detail of the second.

NOTE: this code is in the branch named 2199.

Image 1: detail

I see 4 slightly darker smudges in this view. I will look at the gap-filling code to see why that happens. (Am I using local values to fill the gap, or means? If the latter, that might explain the colour mismatch. I think it ought to be using the former.)

Image 2: full view

[dankelley]

@mdupilka and @richardsc -- I think I have this working now in the 2199 branch (commit b9536390d852f2e536bc942e0088e5915ffbe18c).

I hope others (especially @mdupilka, who is trying to use this code) can build from this branch (this or any later commit) and do some tests. If things are okay, I'll merge the 2199 branch into develop. I may do that anyway, if @richardsc and I discuss this in-person in the next few days and decide that this is the right path to follow. (I don't like leaving branches active for too long, because merging back and forth is a bit of a pain.)

My main test is as shown in the reprex shown below. Note that I'm using a semi-transparent colourmap, to make it easier to refer to problems in terms of coastline features.

Changes.

1. I went back to the old system (using C code) for the gridding mechanism that is used if mapImage() is called with filledContour=TRUE. This is preferable to the scheme that I set up in June of 2023 because it can fill gaps (as mentioned previously in this issue).
2. I cleaned up the code a bit to make it simpler for developers to work with.
3. I added an R interface to debugging information, so that if mapImage() is called with debug > 1, the C code will output information about gaps that are being filled, and how. (Previously, the only way to get information on the C code operations was to alter that code and rebuild.)

I ask that testers include the full debugging output, along with any datasets involved. Also, please use output to PNG, because if the problem only comes up in certain window geometries, I would need to somehow guess the exact window size you're using on your screen, and we would need to have the same machine types, etc.

If you click Details below, you'll see a reprex. Notice that I am using debug=3 in the mapImage() call, and that therefore some output is printed about how gaps are filled. The 2 direction(s) part means that the gap was filled across both x and y. I think the code can also fill gaps that are open on one side of the region, and then it would say 1 direction(s). I've not tested that condition, though, and the code is relatively complex because I'm working in C and not C++, so I use preprocessor macros for array lookup, etc.

``` r if (FALSE) { # DEK uses to make plots he can examine very closely png("issue2199_test_for_speed.png", unit = "in", width = 7, height = 7, res = 400 ) } library(oce) #> Loading required package: gsw data("coastlineWorldMedium", package = "ocedata") int <- 0.1 lon <- seq(-130, -90, int) lat <- seq(50, 89.95, int) m <- matrix(lat, byrow = TRUE, nrow = length(lon), ncol = length(lat)) par(mar = c(3, 3, 1, 1)) mapPlot(coastlineWorldMedium, projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110", longitudelim = c(-140, -80), latitudelim = c(65, 75) ) # use semi-transparent colours so we have coastline for reference mapImage( longitude = lon, latitude = lat, z = m, filledContour = TRUE, breaks = 256, col = function(n) paste0(oceColorsJet(n), "aa"), debug = 2 ) ```  #> mapImage(..., missingColor=NA, filledContour=TRUE, gridder=binMean2D, , ...) { #> only 1 break given, so taking that as number of breaks #> col is a function #> zclip:FALSE #> breaks[1:201]: 50.0, 50.2, ..., 89.8, 90.0 #> col[1:200]: #00007Faa, #000084aa, ..., #840000aa, #7F0000aa #> using 'breaks' to pin image colours #> pinned 401 low values and 0 high z values, out of a total of 160400 values #> using filled contours #> before trimming, length(xx): 160400 #> after trimming, length(xx): 67386 #> using binMean2D() #> binMean2D() ... #> calling C code bin_mean_2d #> in bin_mean_2d() with nx=67386 nxbreaks=165 nybreaks=165 fill=1 fillgap=-1 #> mean[47,149]=NA but mean[148,c(150,47)]=(76.8,76.9) yielding 76.85 (sum=76.85 N=1) #> set mean[47,149] to 76.85 #> mean[61,157]=NA but mean[c(60,62),157]=(77.8,77.9) yielding 77.85 (sum=77.85, N=1) #> mean[61,157]=NA but mean[156,c(158,61)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[61,157] to 77.85 #> mean[67,156]=NA but mean[c(66,68),156]=(77.8,77.9) yielding 77.85 (sum=77.85, N=1) #> mean[67,156]=NA but mean[155,c(157,67)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[67,156] to 77.85 #> mean[67,163]=NA but mean[c(66,68),163]=(78.4,78.5) yielding 78.45 (sum=78.45, N=1) #> set mean[67,163] to 78.45 #> mean[68,157]=NA but mean[c(67,69),157]=(77.9,78) yielding 77.95 (sum=77.95, N=1) #> mean[68,157]=NA but mean[156,c(158,68)]=(77.9,78) yielding 77.95 (sum=155.9 N=2) #> set mean[68,157] to 77.95 #> mean[69,150]=NA but mean[c(68,70),150]=(77.3,77.4) yielding 77.35 (sum=77.35, N=1) #> mean[69,150]=NA but mean[149,c(151,69)]=(77.3,77.4) yielding 77.35 (sum=154.7 N=2) #> set mean[69,150] to 77.35 #> mean[70,151]=NA but mean[c(69,71),151]=(77.4,77.5) yielding 77.45 (sum=77.45, N=1) #> mean[70,151]=NA but mean[150,c(152,70)]=(77.4,77.5) yielding 77.45 (sum=154.9 N=2) #> set mean[70,151] to 77.45 #> mean[71,143]=NA but mean[c(70,72),143]=(76.7,76.8) yielding 76.75 (sum=76.75, N=1) #> mean[71,143]=NA but mean[142,c(144,71)]=(76.7,76.8) yielding 76.75 (sum=153.5 N=2) #> set mean[71,143] to 76.75 #> mean[71,159]=NA but mean[c(70,72),159]=(78.1,78.2) yielding 78.15 (sum=78.15, N=1) #> mean[71,159]=NA but mean[158,c(160,71)]=(78.1,78.2) yielding 78.15 (sum=156.3 N=2) #> set mean[71,159] to 78.15 #> mean[72,144]=NA but mean[c(71,73),144]=(76.8,76.9) yielding 76.85 (sum=76.85, N=1) #> mean[72,144]=NA but mean[143,c(145,72)]=(76.8,76.9) yielding 76.85 (sum=153.7 N=2) #> set mean[72,144] to 76.85 #> mean[73,153]=NA but mean[c(72,74),153]=(77.6,77.7) yielding 77.65 (sum=77.65, N=1) #> mean[73,153]=NA but mean[152,c(154,73)]=(77.6,77.7) yielding 77.65 (sum=155.3 N=2) #> set mean[73,153] to 77.65 #> mean[73,160]=NA but mean[c(72,74),160]=(78.2,78.3) yielding 78.25 (sum=78.25, N=1) #> mean[73,160]=NA but mean[159,c(161,73)]=(78.2,78.3) yielding 78.25 (sum=156.5 N=2) #> set mean[73,160] to 78.25 #> mean[74,135]=NA but mean[c(73,75),135]=(76,76.1) yielding 76.05 (sum=76.05, N=1) #> mean[74,135]=NA but mean[134,c(136,74)]=(76,76.1) yielding 76.05 (sum=152.1 N=2) #> set mean[74,135] to 76.05 #> mean[74,145]=NA but mean[c(73,75),145]=(76.9,77) yielding 76.95 (sum=76.95, N=1) #> mean[74,145]=NA but mean[144,c(146,74)]=(76.9,77) yielding 76.95 (sum=153.9 N=2) #> set mean[74,145] to 76.95 #> mean[75,154]=NA but mean[c(74,77),154]=(77.7,77.8) yielding 77.73 (sum=77.73, N=1) #> mean[75,154]=NA but mean[153,c(155,75)]=(77.7,77.8) yielding 77.75 (sum=155.5 N=2) #> set mean[75,154] to 77.74 #> mean[75,161]=NA but mean[c(74,77),161]=(78.3,78.4) yielding 78.33 (sum=78.33, N=1) #> mean[75,161]=NA but mean[160,c(162,75)]=(78.3,78.4) yielding 78.35 (sum=156.7 N=2) #> set mean[75,161] to 78.34 #> mean[76,136]=NA but mean[c(75,77),136]=(76.1,76.2) yielding 76.15 (sum=76.15, N=1) #> mean[76,136]=NA but mean[135,c(137,76)]=(76.1,76.2) yielding 76.15 (sum=152.3 N=2) #> set mean[76,136] to 76.15 #> mean[76,146]=NA but mean[c(75,77),146]=(77,77.1) yielding 77.05 (sum=77.05, N=1) #> mean[76,146]=NA but mean[145,c(147,76)]=(77,77.1) yielding 77.05 (sum=154.1 N=2) #> set mean[76,146] to 77.05 #> mean[76,154]=NA but mean[c(75,77),154]=(77.74,77.8) yielding 77.77 (sum=77.77, N=1) #> mean[76,154]=NA but mean[153,c(155,76)]=(77.7,77.8) yielding 77.75 (sum=155.5 N=2) #> set mean[76,154] to 77.76 #> mean[76,161]=NA but mean[c(75,77),161]=(78.34,78.4) yielding 78.37 (sum=78.37, N=1) #> mean[76,161]=NA but mean[160,c(162,76)]=(78.3,78.4) yielding 78.35 (sum=156.7 N=2) #> set mean[76,161] to 78.36 #> mean[77,124]=NA but mean[c(76,78),124]=(75,75.1) yielding 75.05 (sum=75.05, N=1) #> mean[77,124]=NA but mean[123,c(125,77)]=(75,75.1) yielding 75.05 (sum=150.1 N=2) #> set mean[77,124] to 75.05 #> mean[77,137]=NA but mean[c(76,80),137]=(76.2,76.3) yielding 76.22 (sum=76.22, N=1) #> mean[77,137]=NA but mean[136,c(138,77)]=(76.2,76.3) yielding 76.25 (sum=152.5 N=2) #> set mean[77,137] to 76.24 #> mean[77,162]=NA but mean[c(76,87),162]=(78.4,78.4) yielding 78.4 (sum=78.4, N=1) #> mean[77,162]=NA but mean[161,c(163,77)]=(78.4,78.5) yielding 78.45 (sum=156.8 N=2) #> set mean[77,162] to 78.42 #> mean[78,137]=NA but mean[c(77,80),137]=(76.24,76.3) yielding 76.26 (sum=76.26, N=1) #> mean[78,137]=NA but mean[136,c(138,78)]=(76.2,76.3) yielding 76.25 (sum=152.5 N=2) #> set mean[78,137] to 76.25 #> mean[78,147]=NA but mean[c(77,86),147]=(77.1,77.1) yielding 77.1 (sum=77.1, N=1) #> mean[78,147]=NA but mean[146,c(148,78)]=(77.1,77.2) yielding 77.15 (sum=154.2 N=2) #> set mean[78,147] to 77.12 #> mean[78,155]=NA but mean[c(77,86),155]=(77.8,77.8) yielding 77.8 (sum=77.8, N=1) #> mean[78,155]=NA but mean[154,c(156,78)]=(77.8,77.9) yielding 77.85 (sum=155.6 N=2) #> set mean[78,155] to 77.82 #> mean[78,162]=NA but mean[c(77,87),162]=(78.42,78.4) yielding 78.42 (sum=78.42, N=1) #> mean[78,162]=NA but mean[161,c(163,78)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[78,162] to 78.44 #> mean[79,125]=NA but mean[c(78,85),125]=(75.1,75.1) yielding 75.1 (sum=75.1, N=1) #> mean[79,125]=NA but mean[124,c(126,79)]=(75.1,75.2) yielding 75.15 (sum=150.2 N=2) #> set mean[79,125] to 75.12 #> mean[79,137]=NA but mean[c(78,80),137]=(76.25,76.3) yielding 76.28 (sum=76.28, N=1) #> mean[79,137]=NA but mean[136,c(138,79)]=(76.2,76.3) yielding 76.25 (sum=152.5 N=2) #> set mean[79,137] to 76.26 #> mean[79,147]=NA but mean[c(78,86),147]=(77.12,77.1) yielding 77.12 (sum=77.12, N=1) #> mean[79,147]=NA but mean[146,c(148,79)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[79,147] to 77.14 #> mean[79,155]=NA but mean[c(78,86),155]=(77.82,77.8) yielding 77.82 (sum=77.82, N=1) #> mean[79,155]=NA but mean[154,c(156,79)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[79,155] to 77.84 #> mean[79,162]=NA but mean[c(78,87),162]=(78.44,78.4) yielding 78.43 (sum=78.43, N=1) #> mean[79,162]=NA but mean[161,c(163,79)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[79,162] to 78.44 #> mean[80,102]=NA but mean[c(79,84),102]=(72.9,72.9) yielding 72.9 (sum=72.9, N=1) #> mean[80,102]=NA but mean[101,c(103,80)]=(72.9,73) yielding 72.95 (sum=145.9 N=2) #> set mean[80,102] to 72.93 #> mean[80,125]=NA but mean[c(79,85),125]=(75.12,75.1) yielding 75.12 (sum=75.12, N=1) #> mean[80,125]=NA but mean[124,c(126,80)]=(75.1,75.2) yielding 75.15 (sum=150.3 N=2) #> set mean[80,125] to 75.14 #> mean[80,147]=NA but mean[c(79,86),147]=(77.14,77.1) yielding 77.13 (sum=77.13, N=1) #> mean[80,147]=NA but mean[146,c(148,80)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[80,147] to 77.14 #> mean[80,155]=NA but mean[c(79,86),155]=(77.84,77.8) yielding 77.83 (sum=77.83, N=1) #> mean[80,155]=NA but mean[154,c(156,80)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[80,155] to 77.84 #> mean[80,162]=NA but mean[c(79,87),162]=(78.44,78.4) yielding 78.44 (sum=78.44, N=1) #> mean[80,162]=NA but mean[161,c(163,80)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[80,162] to 78.44 #> mean[81,102]=NA but mean[c(80,84),102]=(72.93,72.9) yielding 72.92 (sum=72.92, N=1) #> mean[81,102]=NA but mean[101,c(103,81)]=(72.9,73) yielding 72.95 (sum=145.9 N=2) #> set mean[81,102] to 72.93 #> mean[81,125]=NA but mean[c(80,85),125]=(75.14,75.1) yielding 75.13 (sum=75.13, N=1) #> mean[81,125]=NA but mean[124,c(126,81)]=(75.1,75.2) yielding 75.15 (sum=150.3 N=2) #> set mean[81,125] to 75.14 #> mean[81,138]=NA but mean[c(80,83),138]=(76.3,76.3) yielding 76.3 (sum=76.3, N=1) #> mean[81,138]=NA but mean[137,c(139,81)]=(76.3,76.4) yielding 76.35 (sum=152.6 N=2) #> set mean[81,138] to 76.32 #> mean[81,147]=NA but mean[c(80,86),147]=(77.14,77.1) yielding 77.13 (sum=77.13, N=1) #> mean[81,147]=NA but mean[146,c(148,81)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[81,147] to 77.14 #> mean[81,155]=NA but mean[c(80,86),155]=(77.84,77.8) yielding 77.83 (sum=77.83, N=1) #> mean[81,155]=NA but mean[154,c(156,81)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[81,155] to 77.84 #> mean[81,162]=NA but mean[c(80,87),162]=(78.44,78.4) yielding 78.44 (sum=78.44, N=1) #> mean[81,162]=NA but mean[161,c(163,81)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[81,162] to 78.44 #> mean[82,102]=NA but mean[c(81,84),102]=(72.93,72.9) yielding 72.92 (sum=72.92, N=1) #> mean[82,102]=NA but mean[101,c(103,82)]=(72.9,73) yielding 72.95 (sum=145.9 N=2) #> set mean[82,102] to 72.94 #> mean[82,125]=NA but mean[c(81,85),125]=(75.14,75.1) yielding 75.13 (sum=75.13, N=1) #> mean[82,125]=NA but mean[124,c(126,82)]=(75.1,75.2) yielding 75.15 (sum=150.3 N=2) #> set mean[82,125] to 75.14 #> mean[82,138]=NA but mean[c(81,83),138]=(76.32,76.3) yielding 76.31 (sum=76.31, N=1) #> mean[82,138]=NA but mean[137,c(139,82)]=(76.3,76.4) yielding 76.35 (sum=152.7 N=2) #> set mean[82,138] to 76.33 #> mean[82,147]=NA but mean[c(81,86),147]=(77.14,77.1) yielding 77.13 (sum=77.13, N=1) #> mean[82,147]=NA but mean[146,c(148,82)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[82,147] to 77.14 #> mean[82,155]=NA but mean[c(81,86),155]=(77.84,77.8) yielding 77.83 (sum=77.83, N=1) #> mean[82,155]=NA but mean[154,c(156,82)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[82,155] to 77.84 #> mean[82,162]=NA but mean[c(81,87),162]=(78.44,78.4) yielding 78.44 (sum=78.44, N=1) #> mean[82,162]=NA but mean[161,c(163,82)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[82,162] to 78.44 #> mean[83,102]=NA but mean[c(82,84),102]=(72.94,72.9) yielding 72.92 (sum=72.92, N=1) #> mean[83,102]=NA but mean[101,c(103,83)]=(72.9,73) yielding 72.95 (sum=145.9 N=2) #> set mean[83,102] to 72.93 #> mean[83,125]=NA but mean[c(82,85),125]=(75.14,75.1) yielding 75.13 (sum=75.13, N=1) #> mean[83,125]=NA but mean[124,c(126,83)]=(75.1,75.2) yielding 75.15 (sum=150.3 N=2) #> set mean[83,125] to 75.14 #> mean[83,147]=NA but mean[c(82,86),147]=(77.14,77.1) yielding 77.13 (sum=77.13, N=1) #> mean[83,147]=NA but mean[146,c(148,83)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[83,147] to 77.14 #> mean[83,155]=NA but mean[c(82,86),155]=(77.84,77.8) yielding 77.83 (sum=77.83, N=1) #> mean[83,155]=NA but mean[154,c(156,83)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[83,155] to 77.84 #> mean[83,162]=NA but mean[c(82,87),162]=(78.44,78.4) yielding 78.43 (sum=78.43, N=1) #> mean[83,162]=NA but mean[161,c(163,83)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[83,162] to 78.44 #> mean[84,125]=NA but mean[c(83,85),125]=(75.14,75.1) yielding 75.12 (sum=75.12, N=1) #> mean[84,125]=NA but mean[124,c(126,84)]=(75.1,75.2) yielding 75.15 (sum=150.3 N=2) #> set mean[84,125] to 75.13 #> mean[84,137]=NA but mean[c(83,87),137]=(76.3,76.2) yielding 76.28 (sum=76.28, N=1) #> mean[84,137]=NA but mean[136,c(138,84)]=(76.2,76.3) yielding 76.25 (sum=152.5 N=2) #> set mean[84,137] to 76.26 #> mean[84,147]=NA but mean[c(83,86),147]=(77.14,77.1) yielding 77.13 (sum=77.13, N=1) #> mean[84,147]=NA but mean[146,c(148,84)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[84,147] to 77.14 #> mean[84,155]=NA but mean[c(83,86),155]=(77.84,77.8) yielding 77.83 (sum=77.83, N=1) #> mean[84,155]=NA but mean[154,c(156,84)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[84,155] to 77.84 #> mean[84,162]=NA but mean[c(83,87),162]=(78.44,78.4) yielding 78.43 (sum=78.43, N=1) #> mean[84,162]=NA but mean[161,c(163,84)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[84,162] to 78.44 #> mean[85,137]=NA but mean[c(84,87),137]=(76.26,76.2) yielding 76.24 (sum=76.24, N=1) #> mean[85,137]=NA but mean[136,c(138,85)]=(76.2,76.3) yielding 76.25 (sum=152.5 N=2) #> set mean[85,137] to 76.25 #> mean[85,147]=NA but mean[c(84,86),147]=(77.14,77.1) yielding 77.12 (sum=77.12, N=1) #> mean[85,147]=NA but mean[146,c(148,85)]=(77.1,77.2) yielding 77.15 (sum=154.3 N=2) #> set mean[85,147] to 77.13 #> mean[85,155]=NA but mean[c(84,86),155]=(77.84,77.8) yielding 77.82 (sum=77.82, N=1) #> mean[85,155]=NA but mean[154,c(156,85)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[85,155] to 77.83 #> mean[85,162]=NA but mean[c(84,87),162]=(78.44,78.4) yielding 78.43 (sum=78.43, N=1) #> mean[85,162]=NA but mean[161,c(163,85)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[85,162] to 78.44 #> mean[86,124]=NA but mean[c(85,87),124]=(75.1,75) yielding 75.05 (sum=75.05, N=1) #> mean[86,124]=NA but mean[123,c(125,86)]=(75,75.1) yielding 75.05 (sum=150.1 N=2) #> set mean[86,124] to 75.05 #> mean[86,137]=NA but mean[c(85,87),137]=(76.25,76.2) yielding 76.22 (sum=76.22, N=1) #> mean[86,137]=NA but mean[136,c(138,86)]=(76.2,76.3) yielding 76.25 (sum=152.5 N=2) #> set mean[86,137] to 76.24 #> mean[86,162]=NA but mean[c(85,87),162]=(78.44,78.4) yielding 78.42 (sum=78.42, N=1) #> mean[86,162]=NA but mean[161,c(163,86)]=(78.4,78.5) yielding 78.45 (sum=156.9 N=2) #> set mean[86,162] to 78.43 #> mean[87,136]=NA but mean[c(86,88),136]=(76.2,76.1) yielding 76.15 (sum=76.15, N=1) #> mean[87,136]=NA but mean[135,c(137,87)]=(76.1,76.2) yielding 76.15 (sum=152.3 N=2) #> set mean[87,136] to 76.15 #> mean[87,146]=NA but mean[c(86,88),146]=(77.1,77) yielding 77.05 (sum=77.05, N=1) #> mean[87,146]=NA but mean[145,c(147,87)]=(77,77.1) yielding 77.05 (sum=154.1 N=2) #> set mean[87,146] to 77.05 #> mean[87,154]=NA but mean[c(86,89),154]=(77.8,77.7) yielding 77.77 (sum=77.77, N=1) #> mean[87,154]=NA but mean[153,c(155,87)]=(77.7,77.8) yielding 77.75 (sum=155.5 N=2) #> set mean[87,154] to 77.76 #> mean[87,161]=NA but mean[c(86,89),161]=(78.4,78.3) yielding 78.37 (sum=78.37, N=1) #> mean[87,161]=NA but mean[160,c(162,87)]=(78.3,78.4) yielding 78.35 (sum=156.7 N=2) #> set mean[87,161] to 78.36 #> mean[88,154]=NA but mean[c(87,89),154]=(77.76,77.7) yielding 77.73 (sum=77.73, N=1) #> mean[88,154]=NA but mean[153,c(155,88)]=(77.7,77.8) yielding 77.75 (sum=155.5 N=2) #> set mean[88,154] to 77.74 #> mean[88,161]=NA but mean[c(87,89),161]=(78.36,78.3) yielding 78.33 (sum=78.33, N=1) #> mean[88,161]=NA but mean[160,c(162,88)]=(78.3,78.4) yielding 78.35 (sum=156.7 N=2) #> set mean[88,161] to 78.34 #> mean[89,135]=NA but mean[c(88,90),135]=(76.1,76) yielding 76.05 (sum=76.05, N=1) #> mean[89,135]=NA but mean[134,c(136,89)]=(76,76.1) yielding 76.05 (sum=152.1 N=2) #> set mean[89,135] to 76.05 #> mean[89,145]=NA but mean[c(88,90),145]=(77,76.9) yielding 76.95 (sum=76.95, N=1) #> mean[89,145]=NA but mean[144,c(146,89)]=(76.9,77) yielding 76.95 (sum=153.9 N=2) #> set mean[89,145] to 76.95 #> mean[90,153]=NA but mean[c(89,91),153]=(77.7,77.6) yielding 77.65 (sum=77.65, N=1) #> mean[90,153]=NA but mean[152,c(154,90)]=(77.6,77.7) yielding 77.65 (sum=155.3 N=2) #> set mean[90,153] to 77.65 #> mean[90,160]=NA but mean[c(89,91),160]=(78.3,78.2) yielding 78.25 (sum=78.25, N=1) #> mean[90,160]=NA but mean[159,c(161,90)]=(78.2,78.3) yielding 78.25 (sum=156.5 N=2) #> set mean[90,160] to 78.25 #> mean[91,144]=NA but mean[c(90,92),144]=(76.9,76.8) yielding 76.85 (sum=76.85, N=1) #> mean[91,144]=NA but mean[143,c(145,91)]=(76.8,76.9) yielding 76.85 (sum=153.7 N=2) #> set mean[91,144] to 76.85 #> mean[92,143]=NA but mean[c(91,93),143]=(76.8,76.7) yielding 76.75 (sum=76.75, N=1) #> mean[92,143]=NA but mean[142,c(144,92)]=(76.7,76.8) yielding 76.75 (sum=153.5 N=2) #> set mean[92,143] to 76.75 #> mean[92,159]=NA but mean[c(91,93),159]=(78.2,78.1) yielding 78.15 (sum=78.15, N=1) #> mean[92,159]=NA but mean[158,c(160,92)]=(78.1,78.2) yielding 78.15 (sum=156.3 N=2) #> set mean[92,159] to 78.15 #> mean[93,151]=NA but mean[c(92,94),151]=(77.5,77.4) yielding 77.45 (sum=77.45, N=1) #> mean[93,151]=NA but mean[150,c(152,93)]=(77.4,77.5) yielding 77.45 (sum=154.9 N=2) #> set mean[93,151] to 77.45 #> mean[94,150]=NA but mean[c(93,95),150]=(77.4,77.3) yielding 77.35 (sum=77.35, N=1) #> mean[94,150]=NA but mean[149,c(151,94)]=(77.3,77.4) yielding 77.35 (sum=154.7 N=2) #> set mean[94,150] to 77.35 #> mean[95,157]=NA but mean[c(94,96),157]=(78,77.9) yielding 77.95 (sum=77.95, N=1) #> mean[95,157]=NA but mean[156,c(158,95)]=(77.9,78) yielding 77.95 (sum=155.9 N=2) #> set mean[95,157] to 77.95 #> mean[96,156]=NA but mean[c(95,97),156]=(77.9,77.8) yielding 77.85 (sum=77.85, N=1) #> mean[96,156]=NA but mean[155,c(157,96)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[96,156] to 77.85 #> mean[96,163]=NA but mean[c(95,97),163]=(78.5,78.4) yielding 78.45 (sum=78.45, N=1) #> set mean[96,163] to 78.45 #> mean[102,157]=NA but mean[c(101,103),157]=(77.9,77.8) yielding 77.85 (sum=77.85, N=1) #> mean[102,157]=NA but mean[156,c(158,102)]=(77.8,77.9) yielding 77.85 (sum=155.7 N=2) #> set mean[102,157] to 77.85 #> mean[116,149]=NA but mean[148,c(150,116)]=(76.8,76.9) yielding 76.85 (sum=76.85 N=1) #> set mean[116,149] to 76.85 #> bin_mean_2d() encountered (and filled) 86 internal gaps #> setting up return value #> ... binMean2D() completed without problems #> } # mapImage() ^{Created on 2024-03-03 with [reprex v2.1.0](https://reprex.tidyverse.org)}

[mdupilka]

I will look into building from the branch and testing. I have never done that.

[dankelley]

@mdupilka I thought I had desmudged it, but I was wrong! I'd say don't bother building from the branch until I post here that all is okay. That might be tomorrow, though -- somehow I am not seeing the bug. If this continues, I may recode it using Rcpp, which lets me use matrices etc, which will simplify all the C pointers and preprocessor macros and such.

[mdupilka]

I am curious, Is this just an issue within OCE? Or does that bin_mean_2D function (or other) that you use affect any other packages that use it as well?

[dankelley]

This would affect any function in any library (or not in any library) that uses binMean2D() directly, or that uses mapImage() with filledContour set to TRUE.

It's all about filling gaps, which is something that binMean2D() can do (perhaps incorrectly -- that's what I'm investigating).

I switched from C code to standard R code involving split() and cut() because that is easier for users to understand. There are some extra advantages, e.g. a user can decide whether intervals are to be open on the left side or the right side, etc. My C code just makes a decision on that, basically by deciding whether to use < or <= on one side or the other.

The other good thing about using R for such work is that for testing a person can just do e.g. source("R/bin.R"); source("myfile.R") and the test is very fast, as compared to building the package which takes maybe 5 minutes.

I need to set this aside for a day to percolate. Maybe I should revert to the version as of last week, which will have the weird pattern but otherwise will perform as expected for other tasks. Then it's just a matter of filling gaps. That would likely have to be in C but it would be new so I'd use Rcpp (letting me use matrices etc, not all a whack of pointers all over the place) which might be easier to get right. That approach would keep other behaviour the same, which is an advantage for the main use case (which almost certainly is binMean2D(), not mapImage() with filledContour=TRUE).

I am at the point where I don't even quite see why I am gridding, since polygon() doesn't need that. But there is likely a reason, I guess. This is why I need to leave it aside. Often I see Clark on Mondays or Wednesdays, and maybe we could chat about it over my blackboard.

[mdupilka]

Thanks for the info. I have been using the filledContour because it adds some nice smoothing, and you say it is faster. But for my operational programs, I use a grid spacing of 0.1 deg which is quite fine. Not using the filledContour does not really make too much of a difference. However, I do have other non-operational plotting with grids of 2.5 to 5 deg. It makes a big difference there not using the filledContour.

[dankelley]

I'm looking at the code right now.

I think (as of today) that the best plan is:

1. Revert the code to what it was last week, so that binMean2D() etc. work as advertised in the present CRAN version. This lets folks decide whether the left-side or right-side of bins ought to be open, etc. Having such choices will be familiar to R users, and I think it is a good thing anyway.

2. This leaves the problem of gaps, which I think is a general thing that might come up a lot in other contexts, so the right scheme might be to abstract that into another function, which could be called by binMean2D() etc after doing the cut() & split() processing in R.

This leaves 2 mysteries and a decision.

• Mystery 1: why are the gaps showing up in the northern region, where there will be more points (given a grid in lon-lat) to define average-within-patch in x-y space? I just don't understand that.

• Mystery 2: why is my present gap-filling scheme making a different colour in the filled regions? The code can now spit out the interpolation, and it seems right. Maybe this is a plotting thing with the filled-contour code.

• Decision: how to fill gaps? Right now I do something simple. If the i,j point is NA, then I look across i to find two non-NA points, and interpolate linearly for one estimate. Then I do the same across j to get a second estimate. Then I average the two estimates. Neither one will give an estimate if the gap is too wide, though, so a slice of NA values will use adjacent points, not points far away along the slice. Now, this scheme is crude, throwing away other nearby data. But it is easy to code, and gets around problems of how deep to go into surrounding values in all directions, and how to use those values to get a patch. Ultimately, it's a matter of making up data, and there is no "right" way so I just did something simple, realizing that in most cases gaps will consist of small regions.

[dankelley]

I wrote a new R function called fillGapMatrix() to do gap-filling of the gridded matrix that is needed for this issue (that is, for mapImage() called with filledContour equal to TRUE). At https://youtu.be/cY504JBZeDQ I've put a little explainer for how it works. Of course, it is also documented.

This function is in the 2199 branch, but NOTE that this branch is not ready for building yet, and I don't recommend that anybody try. The point of this comment is just to say that I am starting to execute my plan of gaining some success in the present issue, without breaking any old code.

PS. the function works by calling C++ code. This is a lot easier to write than C code, because the latter doesn't not have a matrix type, whereas C++ using Rcpp has some nice types for matrices and other things.

PPS. The code is peppered with debugging output, but that won't show up unless fillGap() is called with debug=1 set.

[dankelley]

I think I have this working now, in the 2199 branch. If you click on Details below, you'll see R code that produces a graph with similar geometry as that in the original comment by @mdupilka. (That comment had no png() specifications so I had to guess on the geometry. That's why the graph is a bit different.) In the code, note the use of a user-defined "gridder" function, which fills in NA-valued regions (AKA gaps).

NOTES:

1. The gaps get filled in if the user creates a function to do the gridding. That is illustrated in the R code.
2. The docs for mapImage() now explain how to work with a user-supplied gridding function.
3. I got rid of the test C code that I had revived in working on this over the weekend. Now, the function behaviours (apart from this fix to mapImage(), etc) are as on CRAN. This is important to me, and explains the delay in getting to a solution.
4. I don't understand the gaps. Maybe @richardsc and I can talk about that a bit tomorrow, if he has time before his class.
5. I will make a comment immediately after posting this, which is for another example.

**Code** ```R PNG <- function(name) { if (!interactive()) { # For the png settings, see # https://codedocean.wordpress.com/2014/02/03/anti-aliasing-and-image-plots/ png(name, type = "cairo", antialias = "none", family = "Arial", unit = "in", width = 7, height = 3.5, res = 400, pointsize = 8 ) par(mar = c(2, 2, 1, 1), mgp = c(2, 0.7, 0)) } } GNP <- function() { if (!interactive()) { dev.off() } } library(oce) data("coastlineWorldMedium", package = "ocedata") int <- 0.1 lon <- seq(-130, -90, int) lat <- seq(50, 89.95, int) m <- matrix(lat, byrow = TRUE, nrow = length(lon), ncol = length(lat)) PNG("2199a_default.png") mapPlot(coastlineWorldMedium, projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110", longitudelim = c(-140, -80), latitudelim = c(65, 75) ) g <- function(...) binMean2D(..., fill = TRUE) mapImage( longitude = lon, latitude = lat, z = m, filledContour = TRUE, breaks = 256, # col = function(n) paste0(oceColorsJet(n), "aa") col = oceColorsJet ) GNP() PNG("2199a_custom_gridder.png") mapPlot(coastlineWorldMedium, projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110", longitudelim = c(-140, -80), latitudelim = c(65, 75) ) g <- function(...) binMean2D(..., fill = TRUE) mapImage( longitude = lon, latitude = lat, z = m, filledContour = TRUE, breaks = 256, gridder = g, # col = function(n) paste0(oceColorsJet(n), "aa") col = oceColorsJet ) GNP() ``` **Resultant plot 1: default** Notice the ugly moire-type pattern.  **Resultant plot 2: with custom gridding function** Note that the pattern is no longer present, because the gaps got filled in. 

[dankelley]

This is similar to the previous comment so I won't explain much. I will say, though, that the "smudge" problem that had consumed me was not actually caused by missing values. It was because I was fiddling with the colour palette, making semi-transparent colours. My hypothesis is that some values just happened to be difficult to render. (My machine is running a beta OS, and some rendering elements are changing from day to day.)

**Code** ```R PNG <- function(name) { if (!interactive()) { # For the png settings, see # https://codedocean.wordpress.com/2014/02/03/anti-aliasing-and-image-plots/ png(name, type = "cairo", antialias = "none", family = "Arial", unit = "in", width = 7, height = 3.5, res = 400, pointsize = 8 ) par(mar = c(2, 2, 1, 1), mgp = c(2, 0.7, 0)) } } GNP <- function() { if (!interactive()) { dev.off() } } library(oce) data("coastlineWorldMedium", package = "ocedata") int <- 0.1 lon <- seq(-130, -90, int) lat <- seq(50, 89.95, int) m <- matrix(lat, byrow = TRUE, nrow = length(lon), ncol = length(lat)) PNG("2199a_default.png") mapPlot(coastlineWorldMedium, projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110", longitudelim = c(-140, -80), latitudelim = c(65, 75) ) g <- function(...) binMean2D(..., fill = TRUE) mapImage( longitude = lon, latitude = lat, z = m, filledContour = TRUE, breaks = 256, # col = function(n) paste0(oceColorsJet(n), "aa") col = oceColorsJet ) GNP() PNG("2199a_custom_gridder.png") mapPlot(coastlineWorldMedium, projection = "+proj=lcc +lat_1=30 +lat_2=45 +lon_0=-110", longitudelim = c(-140, -80), latitudelim = c(65, 75) ) g <- function(...) binMean2D(..., fill = TRUE) mapImage( longitude = lon, latitude = lat, z = m, filledContour = TRUE, breaks = 256, gridder = g, # col = function(n) paste0(oceColorsJet(n), "aa") col = oceColorsJet ) GNP() ``` **First graph, without gap filling** Notice the white regions.  **Second graph, with gap-filling** 

[dankelley]

PS to @mdupilka and @richardsc the R code is in the oce-test repo, which is worth cloning, to keep up with things. It is at https://github.com/dankelley/oce-issues in the 21xx directory, subdirectory 2199. On a unix/mac machine, just type make there to run the tests. Build oce from the 2199 branch first, of course.

That's it for me for today. Thanks for the patience.

[mdupilka]

I cloned the oce from the 2199 branch, built and installed it. I ran it on my example that I started the issue with and it did the job of filling in the weird pattern. I will keep testing it on my other projects.

[mdupilka]

If I now set filledContour = FALSE, the plot is completely screwed up. This is what I get whether using the gridder function or "binMean2D". I get this type of pattern in any of my other projects. Hopefully there is an easy fix to this, as it is unusable right now.

[dankelley]

Thanks for the test. My guess is that it will be an easy fix, but please supply the code that made the plot. (For any issue, please supply code.)

This looks like something funny is going on with the projection. But it cannot be totally broken because contouring still works. Please see the reprex below. By the way, I made that simply by copying the sample code from ?mapContour and then typed reprex::reprex(). It's easy to make a self-contained test and it can be helpful since people work in different time zones. I know you won't see this comment for perhaps 6 hours, so I'll try to make a reprex myself.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
if (requireNamespace("ocedata", quietly=TRUE)) {
    data(levitus, package="ocedata")
    par(mar=rep(1, 4))
    mapPlot(coastlineWorld, projection="+proj=robin", col="lightgray")
    mapContour(levitus[["longitude"]], levitus[["latitude"]], levitus[["SST"]])
}

^{Created on 2024-03-06 with reprex v2.1.0}

[dankelley]

Here is a reprex I made from the code in ?mapImage. I'll use this as a base. I suspect it worked before. My first step will be to try this example with the CRAN version.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
data(topoWorld)

# Northern polar region, with color-coded bathymetry
par(mfrow=c(1,1), mar=c(2,2,1,1))
cm <- colormap(zlim=c(-5000, 0), col=oceColorsGebco)
drawPalette(colormap=cm)
mapPlot(coastlineWorld, projection="+proj=stere +lat_0=90",
        longitudelim=c(-180,180), latitudelim=c(70,110))
mapImage(topoWorld, colormap=cm)
mapGrid(15, 15, polarCircle=1, col=gray(0.2))
mapPolygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col="tan")

^{Created on 2024-03-06 with reprex v2.1.0}

[dankelley]

Here is what I get from the CRAN version. Next, I'll try building from the "develop" branch.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
data(topoWorld)

# Northern polar region, with color-coded bathymetry
par(mfrow=c(1,1), mar=c(2,2,1,1))
cm <- colormap(zlim=c(-5000, 0), col=oceColorsGebco)
drawPalette(colormap=cm)
mapPlot(coastlineWorld, projection="+proj=stere +lat_0=90",
        longitudelim=c(-180,180), latitudelim=c(70,110))
mapImage(topoWorld, colormap=cm)
mapGrid(15, 15, polarCircle=1, col=gray(0.2))
mapPolygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col="tan")

A
#> Error in eval(expr, envir, enclos): object 'A' not found

^{Created on 2024-03-06 with reprex v2.1.0}

[dankelley]

This is with the "develop" branch. This confirms that the problem is only in the "2199" branch.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
data(topoWorld)

# Northern polar region, with color-coded bathymetry
par(mfrow=c(1,1), mar=c(2,2,1,1))
cm <- colormap(zlim=c(-5000, 0), col=oceColorsGebco)
drawPalette(colormap=cm)
mapPlot(coastlineWorld, projection="+proj=stere +lat_0=90",
        longitudelim=c(-180,180), latitudelim=c(70,110))
mapImage(topoWorld, colormap=cm)
mapGrid(15, 15, polarCircle=1, col=gray(0.2))
mapPolygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col="tan")

^{Created on 2024-03-06 with reprex v2.1.0}

[dankelley]

This has been fixed in the "2199" branch, commit f5abf04d446680c50da422688b399d230e6d6974 (if you click on that link, you'll see that I had tried to simplify some code without really looking closely ... silly me).

@mdupilka please git-pull the "2199" branch and rebuild, then test on your own example. If it is still broken, please supply a "reprex" example. If it's OK, just a note stating that would be fine. I hope to close this issue this week, because branches tend to diverge and it can be a lot of work cleaning them up.

library(oce)
#> Loading required package: gsw
data(coastlineWorld)
data(topoWorld)

# Northern polar region, with color-coded bathymetry
par(mfrow=c(1,1), mar=c(2,2,1,1))
cm <- colormap(zlim=c(-5000, 0), col=oceColorsGebco)
drawPalette(colormap=cm)
mapPlot(coastlineWorld, projection="+proj=stere +lat_0=90",
    longitudelim=c(-180,180), latitudelim=c(70,110))
mapImage(topoWorld, colormap=cm)
mapGrid(15, 15, polarCircle=1, col=gray(0.2))
mapPolygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]],
    col="tan")

^{Created on 2024-03-06 with reprex v2.1.0}

[mdupilka]

I built the new 2199 code and ran it on my major project. All looks to be working with the filledContour using the function for the gridder and using both TRUE and FALSE. I plot a number of different fields, some using the filling and some not. It looks back to the normal way things were a month or so ago. Thanks for all the work. Hopefully, you can eventually find out why the gridding changed to produce the weird pattern.

[dankelley]

CLARK -- question for you at the end.

The reason has to do with whether there are any lon-lat points within a given x-y grid cell.

We represent the provided (lon,lat) data on an (x,y) space. But the resultant (x,y) values are not on a grid. Therefore, we set up a regular (x',y') grid, and use binAverage2D() to assign the (x,y) data to appropriate cells.

Depending on the project, the zoom level, etc., there can be (x',y') cells that do not contain any (x,y) data. So they are set to NA.

The code decides (see https://github.com/dankelley/oce/blob/a9d1c17653bde3315e7a6c7b68eaed8d89ef7dae/R/map.R#L3635) on how fine to make the (x',y') grid based on the data density in the view, etc. No user control is provided for that. A coarser grid would pick up more data, but then fine-scale variation would be lost. For large scale views, there can be (depending on the projection) quite a lot of variation in the range of (x,y) cell size.

QUESTION FOR CLARK, @richardsc --

I think maybe I should add another parameter to mapImage() to control that cell size. It would end up in the definition of NN, near the above-mentioned link. Do you object to that? (I can make it so that the default will give the present results, to avoid behaviour changing for users who do not have the gap problem.)

[dankelley]

I added a parameter gridCoarseness to mapImage() in the 2199 branch, commit d0a05d733a5eb8f6eb4901865407ea739ce7db0b.

If you pull the oce-issues repo and type make there, you'll get files 2199c_1.png, 2199c_2.png, etc. I'm pasting them below, under Details. The demonstration shows how altering the new gridCoarseness parameter can fill in blanks, but also coarsens the grid. The top margin shows the value of this parameter. If not supplied, 10 is used.

I ask that @mdupilka rebuild, check on the application at hand, and also take a look at the updated docs (type ?mapImage and scroll to gridCoarseness to see if my text makes any sense).

Perhaps @mdupilka can consider closing the issue. I don't like leaving two active oce branches for too long, because merge conflicts come up when that is done.

    

[richardsc]

This seems like a good idea. Though, I'm not sold on the argument name -- to me "coarseness" implies that a higher number will have a lower resolution (e.g. coarse gravel has fewer, larger, rocks). Maybe that's actually how it works, here, but regardless it's a little confusing.

Does the number have a specific meaning? (sorry, haven't pulled and rebuilt yet to look at the docs). What about gridSpacing? gridResolution?

[dankelley]

Higher numbers mean coarser grids. I think the best plan is to suggest alterations in the docs as they are. The user is not specifying a spacing or a resolution, just whether to increase or decrease the spacing. (The code counts all the grid cells that are in the current view. Then it divides that by coarseness and takes a square root. That becomes the number of cells to take in each of the x and y directions.

Think of the default, 10, as meaning that we hope to get about 3 lon-lat data within each cell. I picked that number as the default likely for an eastern Canada view, and it was close to OK in the present case, except at high latitudes. The right number will depend on the view and the projection. My guess is that 10 might be ok for quite a few places but now the user can control things. But remember, we are not letting them state the number of x or y values. We could do that, I suppose, but I do like having an automatic value that depends on the data.

I'd prefer to keep the existing behaviour the same by default, but we could do that: if no grid given then do as we do now. Otherwise, if grid is a single number, make a square grid like that. Or if it's two numbers, first is for x and second is for y.

Lots of possibilities could come up. Now's the time to make a decision.

Thanks.

[dankelley]

Come to think of it, if I stick with coarseness I can make the default be 1, not 10. That will be less confusing. I could also put the factor outside the square root. Then it means "factor by which to widen and heighten the grid cells".

[mdupilka]

Are all the changes in effect now? Should I try out the 2199?

[dankelley]

I think it's likely that Clark and I will agree on some alterations to the new gridCoarseness parameter. He doesn't seem to like the name, but I like it so we'll see if I've convinced him.

I think the default value of gridCoarseness ought to be 1 (as discussed above).

Maybe wait until tomorrow to try it. I don't think anything is broken anymore, but the code is still in flux.

[dankelley]

I went ahead and changed so that gridCoarseness has a default of 1, not 10. Inside the code, the value 10 is retained, just now multiplied by gridCoarseness. I thought that a default of 1 would make more sense to users.

The test code at https://github.com/dankelley/oce-issues/blob/main/21xx/2199/2199c.R shows the results of altering this parameter. The results are shown below. Arguably, I should change that factor 10 value, but I don't want to, because that was the old default.

I've also gone through the whole mapImage() documentation to try to make it clearer and less disjoint.

In the test plot, notice that the default (middle panel) has some of the weird patterns, but that increasing the parameter to 1.5 makes them go away. I imagine 1.2 might do the same, but this all depends on the projection and the size of the plotted domain (so it depends on page geometry, for example).

Sorry, you'll need to click on the plot to see it in a larger size. This was a plot that I made to view on my screen, not in the half-screen width that GH has chosen for it's web product.

I have a seminar in a while and a meeting this afternoon but I'll still see emails, so if either @richardsc or @mdupilka want to chime in, I'll notice it. I might not be able to do anything until later today or tomorrow though. My hope is that @mdupilka will test this on real data, and also take a moment to look at the docs to see if they halfway make sense. It's a bit hard to explain this matter of holes in grids without standing at a blackboard (as Clark and I have done) and plotting some highly detailed plots that need zooming on grids (as I have done).

[dankelley]

Oh, I just realized that the docs (but not the code) for mapImage() said that filledContour can be a number. So that means we don't need a new gridCoarseness parameter. I'm changing the code and docs to do that. Funny how doing some chores can clear the head. I will update the preceding demonstration of gridCoarseness in-place, once the code has been built and tested locally.

[dankelley]

@mdupilka -- the 2199 branch is ready for testing now. Please do ?mapImage to see the interface, which I think may be quite useful now. I propose you pull the oce-issues repo and do make clean ; make in the 2199 directory (or however you run R non-interactively on a windows machine, if that's what you use) and look at 2199c.R, which demonstrates the new way of working.

I apologize for the very many comments in this thread. I will be silent now, until I hear back from you.

[mdupilka]

Okey dokey. I will do testing today.

[mdupilka]

I just built the package and am running my project I will do various tests. I see you have added a gap = 2 parameters to the gridder function. I really have no idea of the ins and outs of how the gridder works, and I suspect most users will not either. So I will use the function you suggest in the mapImage documentation. Where might I find the 2199c.R?

[dankelley]

clone github.com/dankelley/oce-issues http://github.com/dankelley/oce-issues then visit 21xx dir, then 2199 dir

On Mar 7, 2024, at 2:18 PM, mdupilka @.***> wrote:

I just built the package and am running my project I will do various tests. I see you have added a gap = 2 parameters to the gridder function. I really have no idea of the ins and outs of how the gridder works, and I suspect most users will not either. So I will use the function you suggest in the mapImage documentation. Where might I find the 2199c.R?

— Reply to this email directly, view it on GitHub https://github.com/dankelley/oce/issues/2199#issuecomment-1984158920, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAAYJDPF2AKA5F3N7OLM5QTYXCVQFAVCNFSM6AAAAABEDJEVLSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSOBUGE2TQOJSGA. You are receiving this because you were assigned.

[mdupilka]

I just got an error message: Error in binMean2D(..., fill = TRUE, gap = 2) : unused argument (gap = 2)

[dankelley]

I'm in a meeting right now. I suggest you rebuild oce from the 2199 branch before testing. But I can't look right now.

[mdupilka]

I downloaded the 2199 branch and did a rebuild and install. But there is still no gap argument. It works without it. For my purposes I am finding a filledContour = 2 is cleaning everything up. I use a lon-lat grid spacing of 0.1 deg. Projection "+proj=lcc +lat_0=63 +lat_1=33 +lat_2=45 +lon_0=-110"

[dankelley]

Sorry, I guess you got that gap error from running the code suggested in the docs for binMean2D(). That should read fillgap = 2, not gap = 2. I'll update the source in a few minutes.

[dankelley]

The 2199 source is now updated.