Reprojecting Rasters

Originally posted: 30 April 2008
*Special thanks to Markus Neteler for guidance in getting this right*

Reprojecting vector maps today is easy, in free and commercial GIS software. Reprojecting raster maps is not.
It can be done for free using the GDAL set of tools. Below, I show how to do that in Linux. I am still trying to figure out how to invoke GDAL tools properly in Windows, and I will update this page once I work it out.

Here I am trying to reproject a clipped section of the MODIS map.
The original map is a beautiful, 173 megabyte TIF file from NASA’s Visible Earth site.

This image?is a radically-reduced jpeg of the original file. Note that GIMP reads this from the top left as 0,0, which equals 90N 180W in geographic coordinates. The Prime Meridian, in pixels, is at 10800,0.

It is 21,600 pixels wide by 10,800 pixels high; in other words, 1×1 degree = 60×60 pixels; each pixel = 1 minute of lat and long.

I can turn the .tif into a GeoTIFF file with this information, using <gdal_translate -a_ullr> at the command-line. Don’t worry if you are not sure it is installed. If you try to run the command and it is not installed, it will advise:

The program ‘gdal_translate’ is currently not installed.? You can install it by typing:
sudo apt-get install gdal-bin
bash: gdal_translate: command not found

The parameters for the -a_ullr option are the upper-left x and y, and then lower-right x and y coordinates. (This is abbreviated on the manpage as: -a_ullr ulx uly lrx lry.) On this file, that is?180w, 90n, 180e, 90s. The lat-long projection, with an ellipsoid and datum of wgs84, has an EPSG code of 4326. Here is the full command:

gdal_translate -a_ullr -180 90 180 -90 -a_srs “EPSG:4326” modis_21600_RGB.tif modis_21600_RGB_geo.tif

I am going to clip this to show the area of central Asia. Clipping boundaries:

North Limit:	50N
South Limit:	20N
West Limit:	25E
East Limit:	95E

Then I can crop the GeoTiff using <gdal_translate -projwin>, creating a new file called casia_modis_geo.tif:

gdal_translate -projwin?25 50 95 20 modis_21600_RGB_geo.tif?casia_modis_geo.tif

The MODIS image, trimmed to a manageable 21 MB geotiff file. (Again, this image is a reduced jpeg)

Now I want to reproject this lat-long image into an equidistant conic projection (eqdc). Equidistant conic? projections have evenly-spaced parallels and straight, radial meridians. North-south scale is consistent everywhere, and east-west scale is also consistent with north-south scale at the standard parallels. Between the two standard parallels, east-west scale is slightly compressed (about 1% in most designs) and outside them, east-west scale is stretched. For a more detailed discussion of this projection, please review the websites of Carlos Furuti , National Geographic, and the USGS. For this map I will specify a southern standard parallel (lat_1) of 32 degrees North; a northern standard parallel (lat_2) of 39 North, and a central meridian of 60 East.

I am going to use the gdalwarp command-line program to do this reprojection. Note that when you install gdal, you actually install a variety of utilities, including gdal_translate used above. You also install gdalwarp, which can reproject files using EPSG codes or PROJ4 definitions. Since equidistant-conics are unusual projections, ideal for a large east-west region, there is no EPSG code that I think will work. Here is the PROJ4 definition of this projection:

+proj=eqdc +lat_1=32n +lat_2=39n +lat_0=35n +lon_0=60e +ellps=WGS84 +datum=WGS84 +no_defs

Here is a brief explanation:

+proj=eqdc??? “eqdc” is the PROJ4 abbreviation for equidistant conic projection.
+lat_1=32n??? The first, or southern standard parallel, is 32 North.
+lat_2=39n??? The second, or northern standard parallel, is 39 North.
+lat_0=35n??? The midline latititude is 35 North (not sure if this is required).
+lon_0=60e??? The standard meridian is 60 East.
+ellps=WGS84??? The spheroid (mathematical model of the shape of the earth) is World Geodetic Survey 1984.
+datum=WGS84??? The altitude standard is also WGS84.

Here are the main options in the gdalwarp command which I will use:

[-s_srs srs_def] = source spatial reference set. Definition can be an EPSG code or a PROJ4 parameter-set. Default assumption is lat-long; if your geotiff is lat-long, you can omit this option.
[-t_srs srs_def] = target spatial reference set. Definition can be an EPSG code or a PROJ4 parameter-set. NOTE! when we place the PROJ4 definition within a gdalwarp command, enclose the definition within single-quotes because the definition contains spaces. Like any command-line program, gdalwarp would interpret the space as a separator between options such as the projection, the resampling style, and the input and output filenames.
[rn,rb,rc,rcs]??? Resampling methods: nearest neighbor, bilinear, cubic,
srcfile ??? This is the source or input file (required).
dstfile??? ?This is the destination or output file (required; cannot simply overwrite source file).

Here is the whole command:

gdalwarp -t_srs ‘+proj=eqdc +lat_1=32n +lat_2=39n +lat_0=35n +lon_0=60e +ellps=WGS84 +datum=WGS84’ -rc casia_modis_geo.tif modis_eqdc_geo.tif

And here is the output, shown in a reduced-size JPEG. The actual output file is 34 MB, and shows rich detail. But for the moment, grab a globe and compare this projection with the actual shapes of India, the Caspian, and the Black Seas. Notice how much more accurate this conic projection is than the lat-long projection above.