grdproject - different -R, same filesize

Added by Andreas 7 months ago

The attached script creates 6 subgrids and projects them, based on the different -R input (country codes). The area varies greatly, but the filesize is (roughly) the same. The bigger area, the lower resolution, to fit into the ~fixed filesize. Why is this?

From a directory listing with filesizes and filenames:

719K DK_earth_relief_15m.nc
798K DK,RU_earth_relief_15m.nc
701K IS_earth_relief_15m.nc
766K NO_earth_relief_15m.nc
786K NO,RU_earth_relief_15m.nc
797K RU_earth_relief_15m.nc

country_relief.sh Magnifier (165 Bytes)


Replies (10)

RE: grdproject - different -R, same filesize - Added by Andreas 7 months ago

Another thing I ran into;

The W/E/S/N from the DCW isnt always spot on. E.g.

$ pscoast -ENO,RU -M | gmtinfo -I1
-R-180/180/41/82
$ pscoast -ECA,RU -M | gmtinfo -I1
-R-180/180/41/84
$ pscoast -ECA -M | gmtinfo -I1
-R218/308/41/84
$ pscoast -ERU -M | gmtinfo -I1
-R19/191/41/82
$ pscoast -EDK,CA,RU -M | gmtinfo -I1
-R-180/180/41/84
$ pscoast -EGL,RU -M | gmtinfo -I1
-R-180/180/41/84
$ pscoast -EGL -M | gmtinfo -I1
-R286/349/59/84
$ pscoast -EGL,CA -M | gmtinfo -I1
-R218/349/41/84
$ pscoast -EGL,US -M | gmtinfo -I1
-R172/349/18/84
$ pscoast -EGL,US,RU -M | gmtinfo -I1
-R-180/180/18/84

The results giving -180/180 as W/E doesnt seem right. -ERU seems to give a good -R, but when RU is in combination with other country codes, trouble arises?!

And last, when giving -RNO,RU to a grdproject command, the grid is cut and does not include the area implied by the given -R (see script and generated figure).

project_NO,RU.sh Magnifier (383 Bytes)

project_NO,RU.pdf (290 KB)

RE: grdproject - different -R, same filesize - Added by Joaquim 7 months ago

I don't find it so strange that the combined NO,RU comes out like that. After all to represent [-10 191] you'll have to choose which longitude limits to use. What is not correct is that it doesn't respect the FORMAT_GEO_OUT

pscoast -ENO,RU -M --FORMAT_GEO_OUT=+D | gmtinfo -I1
-R-180/180/41/82

RE: grdproject - different -R, same filesize - Added by Paul 7 months ago

Well, the FORMAT_GEO_OUT is misplaced (needs to go with gmtinfo) but it still fails...

RE: grdproject - different -R, same filesize - Added by Andreas 7 months ago

1.What about the constant filesize - different resolution thing?

2.The script project.sh takes a country code argument and projects the earth_relief_30m (polar stereographic). The projected grid is cut in many cases, removing a part which should be included according to -R. Maybe the problem is all country code combinations which gives an -R with long -180/180? But then again RU,US plots just fine.

$ pscoast -EDK,RU -M | gmtinfo -I1
-R8/191/41/82 <-- PLOTS FINE
$ pscoast -EGL,RU -M | gmtinfo -I1
-R-180/180/41/84 <-- CUT!
$ pscoast -EIS,RU -M | gmtinfo -I1
-R-180/180/41/82 <-- CUT!
$ pscoast -ENO,RU -M | gmtinfo -I1
-R-180/180/41/82 <-- CUT!
$ pscoast -ERU -M | gmtinfo -I1
-R19/191/41/82 <-- PLOTS FINE
$ pscoast -EUS,RU -M | gmtinfo -I1
-R-180/180/18/82 <-- PLOTS (kind of) FINE

Using 6.0.0_r19803.

project.sh Magnifier (362 Bytes)

project_DK,RU.pdf (382 KB)

project_IS,RU.pdf (415 KB)

project_GL,RU.pdf (419 KB)

project_NO,RU.pdf (435 KB)

project_US,RU.pdf (516 KB)

project_RU.pdf (381 KB)

RE: grdproject - different -R, same filesize - Added by Paul 7 months ago

The key problem we have is that our algorithm for determining overall -R for a set of segments does not work if any segment exceed 180. This only applies to Russia in DCW but could apply to any plate tectonic polygon so not limited to coastlines. The current algorithm keeps track of which quadrant has data and in the GL,RU case we find all 4 quadrants have data and that yields -180/180. It fails to find -74/191 which is the shortest longitude range; it also fails to find 19-348 which is longer but still contains the data. I think a new algorithm needs to do one segment at the time, compare the next segment with what we have so far, and adjust.

I have not looked at your grdproject stuff yet but since you are not specifying any output grid increment I think it will inherit the number of points in the input as nx, ny.

RE: grdproject - different -R, same filesize - Added by Paul 7 months ago

I have implemented a new scheme in r19804 which seems to work fine. All tests pass and give the correct result for your cases. Note:

  1. To just get the DCW region, use pscoast directly, e.g., gmt pscoast -EGL,RU, and append +r<inc> if rounding is needed.
  2. To run via gmt info you will have to add -fg, otherwise it is all Cartesian (we cannot guess that xxx.yyy is a longitude).
  3. You should redo your grdproject cases and see if anything odd happens there. The regions will be correct but you may want to use -D.

RE: grdproject - different -R, same filesize - Added by Andreas 7 months ago

Thank you very much Paul. The area of interest (-R) issues looks really good now.

Another thing I noticed with grids created by grdproject is that the rightmost column of raster cells have funky values (I dont think this is a new issue caused by the recent revisions). See included figures. I used project.sh posted above to create the grids.

Regarding the grdproject stuff; the default behaviour of making the projected grid have the same number of points as the original (as opposed to have the ~correspoding resolution in projected units?), appears to be a bit counter intuitive to me, but can be changed with -D as you pointed out. Thanks!

RE: grdproject - different -R, same filesize - Added by Paul 7 months ago

Thanks Andreas, I found what the problem is: During writing of the projected grid the writer code still thought the grid was geographic and global (360) lon/lat so it applied a rule that said if the eastern col different from the western then we duplicate it (this is fine if your grid is gridline registered and spans 360). It is not a good rule for projected grids in meter. See if r19808 works better for you.

RE: grdproject - different -R, same filesize - Added by Andreas 7 months ago

Great - looks perfect, Paul! Thank you.

I'll try to be quiet for a while now.

RE: grdproject - different -R, same filesize - Added by Andreas 7 months ago

Regarding FORMAT_GEO_OUT not being respected (as stated above); Using -R<country code> (instead of -RW/E/S/N) seems to cause some inconsistency.

$ echo -900000 -2800000 | mapproject -Js0/90/75/1:1 -F -C -I -RIS --FORMAT_GEO_OUT=D
342.181111085   63.666861986
$ echo -900000 -2800000 | mapproject -Js0/90/75/1:1 -F -C -I -RIS --FORMAT_GEO_OUT=-D
342.181111085   63.666861986
$ echo -900000 -2800000 | mapproject -Js0/90/75/1:1 -F -C -I -RIS --FORMAT_GEO_OUT=+D
342.181111085   63.666861986
$ echo -900000 -2800000 | mapproject -Js0/90/75/1:1 -F -C -I -R300/360/60/65 --FORMAT_GEO_OUT=D
-17.8188889145  63.666861986
$ echo -900000 -2800000 | mapproject -Js0/90/75/1:1 -F -C -I -R300/360/60/65 --FORMAT_GEO_OUT=-D
-17.8188889145  63.666861986
$ echo -900000 -2800000 | mapproject -Js0/90/75/1:1 -F -C -I -R300/360/60/65 --FORMAT_GEO_OUT=+D
342.181111085   63.666861986

(1-10/10)