Python API Reproject Raster Coverages

The reprojecting of a RasterCoverage from one coordinate system to another is a commonly needed operation. To make it really simple the Python API of IlwisObjects provides a function that does exactly that. But if you want to know what happens "behind the scenes" further down you can find a step-by-step guide that leads you through the whole process.

.reprojectRaster

The method IlwisObjects provides is called .reprojectRaster() and can be used on an existing RasterCoverage. It will produce a new coverage as output.

As input it needs:

• A name for the new coverage
• The EPSG number of the new coordinate system
• Name of interpolation method (nearestneighbour, bilinear or bicubic)

# load a raster
rc = RasterCoverage("n000302.mpr")
# call the reproject method
rcReproj = rc.reprojectRaster("newraster", 2050, "bilinear")
# set connection for output and store it
rcReproj.store(workingDir+rasterDir+"/aa_testreproj", "GTiff", "gdal")

Step by Step

What we actually need to do to resample a RasterCoverage is call the resample operation with the help of Engine.do(). Resample is a bit more flexible as reproject as it can change pixelsizes and sizes of rasters. But what do we need for the resample operation to work?

What we need	How we get it
Name of RasterCoverage	simply call .name() on your coverage
Kind of interpolation	Just choose from nearestneighbour, bilinear or bicubic and pass it as a string
New GeoReference	The tricky part, will be described in the next steps

Get the new GeoReference

What you need for a new GeoReference:

• EPSG number of new coordinate system
• Coordinates of minimum and maximum corners of the new envelope
• The new size of the grid

The EPSG is the easiest part, as you can just lookt it up on the Internet.

New Envelope

To get the new envelope we first need to get the CoordinateSystem and Envelope of the original raster. We can simply get them by using the methods .coordinateSystem() and .envelope() on our original coverage.

Additionally, we need the new coordinate system. We can create it by calling the constructor CoordinateSystem() and pass it the string with the new EPSG number, e.g. "code=epsg:2050".

Once we have it we can call the .convertEnvelope() function from the original coordinate system and pass it the new coordinate system as well as the original envelope. As output it will give us our new, converted envelope.

sourceCsy = rc.coordinateSystem()
sourceEnv = rc.envelope()

targetCsy = CoordinateSystem("code=epsg:2050")

targetEnv = sourceCsy.convertEnvelope(targetCsy, sourceEnv)

New Grid Size

Now we still have to get the new grid size. For this we will need our new envelope and the original GeoReference, which we get by simply calling .geoReference() on the original coverage.

WE can then use the .envelope2Box() method from the GeoReference to convert our new envelope into a box, from which we can get the size via .size(). From the size we can then get its x and y by using .xsize and .ysize.

sourceGeoref = rc.geoReference()

bo = georef.envelope2Box(targetEnv)
sz = bo.size()

Assembling the new GeoReference

Now we've got all that is necessary to create a new GeoReference, all that is left to do si assemble all the parts togehter in the way IlwisObjects can read them. For that we need to create a string in the form of: "code=georef:type=corners,csy=epsg:2050,envelope=-4.48539e+06 -4.57006e+06 5.11477e+06 4.39656e+06,gridsize=1201 1122,name=grf1" You can freely choose the name of the GeoReference. To get the corners of the new envelope we can simply call it with the native python function str(). As the x and y size of the new box are actually integers, they also have to get turned into strings with the str() function.

newGeoRefStr = ("code=georef:type=corners,csy=epsg:2050,envelope=" + str(targetEnv) + ",gridsize=" + str(sz.xsize) + " " + str(sz.ysize) + ",name=grf1")

targetGeoRef = GeoReference(newGeoRefStr)

Reprojecting the Raster Coverage

As mentioned in the beginning, all we have to do now, is use the Engine.do() function with the resample option and pass all arguments it needs.

# call Engine.do and give it the name of the operation, the coverage's name, the name of the new georef and the interpolation method
rcReproj = Engine.do("resample", rc.name(), targetGeoRef.name(), "bilinear")

This is how you reproject a RasterCoverage in IlwisObjects, now you've got a new coverage for further processing or saving to disk.