How to reproject Natural Earth

Added by Dariusz 10 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 10 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 10 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 10 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 10 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 10 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 10 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 10 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 10 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 10 months ago

GMT6 is on the way ...

(1-9/9)