How to reproject Natural Earth

Added by Dariusz 7 months ago

For example NE1_50M_SR_W.tif file from this package:
http://www.naturalearthdata.com/downloads/50m-raster-data/50m-natural-earth-1/ :
http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/raster/NE1_50M_SR_W.zip

GMT can produce ps map file this way:

gmt psimage NE1_50M_SR_W.tif -K -Dx0/0+w10i > NE1.ps
gmt psbasemap -O -R-180.0/-90.0/180.0/90.0r -JQ10i -Bafg >> NE1.ps
I don't know whether -R-180.0/-90.0/180.0/90.0 -JQ10i is correct, but final ps image seems to look well.
But how to use NE1_50M_SR_W.tif image as background for maps in other projections, like Albers conic equal-area projection or orthographic projection projection (such globe). This .tif background should be reprojected somehow. In second case wrapped around the globe.
Thank you.


Replies (9)

RE: How to reproject Natural Earth - Added by Andreas 7 months ago

Just use the tif file as you would use any other geo grid with e.g. grdimage and the -D option (http://gmt.soest.hawaii.edu/doc/latest/grdimage.html#d).

grdimage -D -P NE1_50M_SR_W.tif -R218/308/41/84 -JS270/90/15c  -Bafg -BWeSn > CA.ps

RE: How to reproject Natural Earth - Added by Dariusz 7 months ago

This doesn't work - produces only grid with coordinates on the white background. No Natural Earth.
That Natural Earth .tif file seems to be only graphics, like .bmp or .jpg .
Perhaps this .tif could be converted somehow into .grd file or similar before processing ? Is it possible to convert regular graphics into .grd file ?

RE: How to reproject Natural Earth - Added by Andreas 7 months ago

Oh, I see. Which version are you using?

My example was based on GMT6.0.0_r19808, but I see now that the same command under GMT 5.4.2 just gives a frame with white background, as you described. The PS file is huge though..

The .tif file is good.

$ grdinfo NE1_50M_SR_W.tif
NE1_50M_SR_W.tif: Title: Grid imported via GDAL
NE1_50M_SR_W.tif: Command:
NE1_50M_SR_W.tif: Remark:
NE1_50M_SR_W.tif: Pixel node registration used [Cartesian grid]
NE1_50M_SR_W.tif: Grid file format: gd = Import/export through GDAL
NE1_50M_SR_W.tif: x_min: -180 x_max: 180 x_inc: 0.0333333333333 name: x n_columns: 10800
NE1_50M_SR_W.tif: y_min: -90 y_max: 90 y_inc: 0.0333333333333 name: y n_rows: 5400
NE1_50M_SR_W.tif: z_min: 62 z_max: 255 name: z
NE1_50M_SR_W.tif: scale_factor: 1 add_offset: 0
+proj=longlat +datum=WGS84 +no_defs

RE: How to reproject Natural Earth - Added by Dariusz 7 months ago

Attempt to cut off a smaller area:

gmt grdcut NE1_50M_SR_W.tif -R-10/30/5/65 -V -GNE1_grdcut.grd
grdinfo NE1_grdcut.grd
grdinfo:
NE1_grdcut.grd: Title: Produced by grdcut
NE1_grdcut.grd: Command: grdcut NE1.tif -R-10/30/5/65 -V -GNE1_grdcut.grd
NE1_grdcut.grd: Remark:
NE1_grdcut.grd: Pixel node registration used [Cartesian grid]
NE1_grdcut.grd: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
NE1_grdcut.grd: x_min: -10 x_max: 30 x_inc: 0.0333333333333 name: x n_columns: 1200
NE1_grdcut.grd: y_min: 5 y_max: 65 y_inc: 0.0333333333333 name: y n_rows: 1800
NE1_grdcut.grd: z_min: 94 z_max: 255 name: z
NE1_grdcut.grd: scale_factor: 1 add_offset: 0
NE1_grdcut.grd: format: netCDF-4 chunk_size: 134,129 shuffle: on deflation_level: 3
GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]]
    UNIT["Degree",0.0174532925199433,
        AUTHORITY["EPSG","91122"]],
    AUTHORITY["EPSG","4326"]]
So .grd is produced, but grdimage doesn't work with this new file:
grdimage -D NE1_grdcut.grd -R-10/30/5/65 -JS10/35/10i -Bafg -BWeSn > NE1_grdimage.ps
grdimage: Using data type other than byte (unsigned char) is not implemented
grdimage (api_import_image): Map region exceeds 360 degrees [NE1_grdcut.grd]
Error returned from GMT API: GMT_IMAGE_READ_ERROR (22)

NE1_50M_SR_W.zip package contains i. e. such additional files:

NE1_50M_SR_W.prj

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]
NE1_50M_SR_W.tfw:
          0.03333333333333 
          0.00000000000000 
          0.00000000000000 
          -0.03333333333333 
          -179.98333333333333 
          89.98333333333333

RE: How to reproject Natural Earth - Added by Paul 7 months ago

You cannot pass a netcdf grid via the -D option. Remove the -D. But now it is just a z-grid and you will need to consider coloring - which is not what you started off to do. We will see if this is an easy fix in 5.4.x since it works in GMT6.

RE: How to reproject Natural Earth - Added by Joaquim 7 months ago

This option has worked on GMT for ages (e.g. my famous Jessica Rabbit example). Haven't had time to check with 5.4.3 but notice that correct syntax is -Dr, not -D

grdimage -Dr ...

RE: How to reproject Natural Earth - Added by Dariusz 7 months ago

Joaquim wrote:

grdimage -Dr ...

The original .tif file: Natural Earth appears on the map, but doesn't fit map generated by pscoast (-R -J), is much smaller.
The .grd file generated with grdcut: The same errors as previously.

RE: How to reproject Natural Earth - Added by Joaquim 7 months ago

Andreas, yes I can confirm that your command works fine on 6 but fails on 5. It probably never worked on 5 or at least I never remember to have tested it to plot only a sub-region of the image. It works fine though on 5 too if we plot the entire region (-Rg). So let's say this is an improvement behavior in GMT6.

RE: How to reproject Natural Earth - Added by Dariusz 7 months ago

GMT6 is on the way ...

(1-9/9)