wrapping data creates long cells

[danlipsa]

If we use the following script (provided by @aashish24 ) that adds two wrapped longitudinal values 180 (same as -180) and 185 (same as -175) we get long cells around the wrapping as seen in the attached picture. See the attached script:

test_robinson_wrap.txt

Here is the image:

[danlipsa]

@aashish24 @doutriaux1 This should work only for cylindrical projections such as mercator isn't it? But it does not work for that either - see the following image:

[danlipsa]

@aashish24 @doutriaux1 I think mercator should work, but I am not sure I have generated the graphics method correctly. Here is the image for lambert projection which seems to work correctly:

[danlipsa]

@aashish24 @doutriaux1 can you paste a script that shows the bug. I think the projected data should be a rectangle. Sample projections are lambert (which works) or mercantor. It is not clear to me what should happen to the wrapped data if the projected data is not a rectangle such as for the robinson projection.

[doutriaux1]

@danlipsa i think it's an issue with proj4 converting 365 -> 5 hence the cel lthat used to go

355 -> 365 now looks like 5 -> 355 hence the big blocks so i think depending on projection it sometimes happens sometimes not.

[danlipsa]

@doutriaux1 It is not clear to me what the correct result should be? What if we have data from -180 to 540 (=180 + 360)? Should we have two globes side by side?

[aashish24]

@doutriaux1 It is not clear to me what the correct result should be? What if we have data from -180 to 540 (=180 + 360)? Should we have two globes side by side?

Data and vis should be duplicated in that case is my understanding. Just like google maps. If you zoom out and can see more than -180 to 180 you will see US on the right side of Asia.

[doutriaux1]

yes for regular porj it is duplicated, but some projection bring everything back between -180/180 (polar for exampe)

[aashish24]

Data and vis should be duplicated in that case is my understanding. Just like google maps. If you zoom out and can see more than -180 to 180 you will see US on the right side of Asia.

I just want to clarify. If we are using projects that wrap data, may be we just need to discard the data? Or constain user to chose only values that makes sense for a projection?

[danlipsa]

@aashish24 In google maps the projected globe is a rectangle. In my mind appending data and showing that in the map makes sense only in that case - when the projected globe is a rectangle. If the projected globe is a ellipsoid, what do you do? I never seen one ellipsoid and a half side by side - this is what you would get for -180 to 360.

[danlipsa]

@doutriaux1 @aashish24 Can I assume that longitude is between -180, 180 eventually extended at the high end, low end or both ends?

[doutriaux1]

I kind of want to say no... Users are going to send you data between 0/360 as well, now I'm not sure what you real mean by longitude is between -180, 180 I guess it depends at which stage you're talking about, I think what proj4 does is that yes, try to put everything between -180/180 (maybe in different units like m for example)

[doutriaux1]

@danlipsa if i remember correctly robinson with data 0/360 should look different from robinson with data between -180/180 (centered around 180 vs 0)

[doutriaux1]

@danlipsa using the following script with 1.5.1

import vcs
import cdms2

f=cdms2.open("/usr/local/uvcdat/1.5.1/sample_data/clt.nc")

s1=f("clt",longitude=(0,360))
s2=f("clt",longitude=(-180,180))

x=vcs.init()
iso = x.createisofill()

iso.projection="robinson"

x.plot(s1,iso)
x.png("0-360")
x.clear()
x.plot(s2,iso)
x.png("180-180")

we used to get:

0-360

-180-180

[danlipsa]

Thanks @doutriaux1 So it seems that the middle of the range is used as the central meridian. What I am trying to do is to restrict the longitude to 360 degrees for projections that cannot show duplicate data, such as robinson. What I can do then remove the extra degrees from the right side and print a warning. The users might end-up with a different direction of view than what they wanted.

[doutriaux1]

@danlipsa, I don't think that's the way to go at it. The proper way is to look at the cell's vertice direction (clockwise vs counterclockwise) if the direction changes it means the cell is being wrapped and we should not spread across the page.

355 --- 365                                             355 --- 5          5 -- 355
  |            |                                                 |         |   ==    |       |
355 --- 365                                             355 -- 5           5 --  355

left is rotating clockwise right counterclockwise, so right should be split in two cells :

5--0         and 355 -- 360

Does this make sense?

[danlipsa]

@doutriaux1 Thanks for the suggestion. Maybe that is a better way to approach this, as I don't need to know which projections handle duplicate data and which projections don't. I will take a look at the wrapping.

[danlipsa]

@doutriaux1 Not sure if this is going to work for every projection. For instance for aeqd as longitude increases the projected points go around a circle.

[doutriaux1]

We can start to implement it for the projection we know of (robinson, molleweide, etc..) when it works make it do this for all and see if it breaks anything.

[danlipsa]

@doutriaux1 Is there any other way in uvcdat to specify the central meridian besides the middle of the range of longitude?

[danlipsa]

@doutriaux1 Another question. If the longitude range is bigger than 360 where is the central meridian? Is is the middle of the range of longitude as well?

[doutriaux1]

@danlipsa yes there is a way

import vcs
x=vcs.init()
p=x.createprojection()
p.type = "robinson"
p.list()

output:

In [4]: p.list()
 ---------- Projection (Proj) member (attribute) listings ---------
secondary method = Proj
name = __projection_458848618473684
type = robinson
sphere = 1e+20
centralmeridian = 1e+20
falseeasting = 1e+20
falsenorthing = 1e+20

[doutriaux1]

1.e20 means the user did not define it, in the gctp packages that meant let gctp decide. I'm not sure in proj4. I think we do middle of longitude range.

[aashish24]

@doutriaux1 @danlipsa and I found a solution. @danlipsa will be pushing a fix shortly :smile:

[doutriaux1]

That's awesome! Thx.

[danlipsa]

@doutriaux1 @aashish24 Do we need a test for this one? Otherwise I'll close it.

[doutriaux1]

I think we should add one yes, otherwise the bug might come back and we won't notice it.

[danlipsa]

Sounds good. @aashish24 Can we sneak this in for the things to do for the next two weeks? Should not take long.

[aashish24]

Sounds good. @aashish24 Can we sneak this in for the things to do for the next two weeks? Should not take long.

Lets finish what we have and then we can add this one in if we ran out of tasks.

[aashish24]

I added 2.6 and Backlog.

[durack1]

@danlipsa I've encountered some similar issues with the ocean grids - so for the models: CMCC-CESM, CMCC-CM, CMCC-CMS, CNRM-CM5, MPI-ESM-MR, MPI-ESM-P, MRI-CGCM3 and I've attached a couple of examples of these below:

[danlipsa]

@durack1 What is the data and the vcs script you use to create the plot?

[durack1]

Would a tarball with these data and a demo script be helpful? If yes I'll get onto that tomorrow.

It would be great to more systematically test across all the available data grids that we commonly use.. These could be incorporated into a more comprehensive test suite data pool

[danlipsa]

@durack1 Yes, that would be great. You could attach the file to this issue. If it turns out to be a different bug, which is possible, we can create a different issue that better describe the problem.

[durack1]

@danlipsa so here are some files, and you can simply use:

import glob,vcs
import cdms2 as cdm
# Ocean grid - basin
files = glob.glob('*.nc')
for i in files:
    c = vcs.init()    
    fH = cdm.open(i)
    tmp = fH('basin')
    c.plot(tmp)
    outFileName = '.'.join([i.split('.')[:-2]),'png'])
    print 'Saving: ',outFileName
    c.png(outFileName)
    c.close()
    fH.close()

basin.tar.gz

[durack1]

@danlipsa let me know if you need any more input here - I hope you can reproduce my issue.

[danlipsa]

@durack1 Thanks Paul. I will take a look at this let you know.

[danlipsa]

I'll close this as I opened https://github.com/UV-CDAT/uvcdat/issues/1848 instead.