\n",
"Find the longitude and latitude of some places of interest in the British isles (West of Greenwich!) and using the MODLAND MODIS tile calculator and the geotransform, repeat the above experiment. Note that the [MODIS calculator](http://landweb.nascom.nasa.gov/cgi-bin/developer/tilemap.cgi) calculates both the projected coordinates in the MODIS sinusoidal projection, as well as the pixel number, so it is a helpful way to check whether you got the right result.\n",
"
\n",
" | \n",
"
\n",
"\n",
"| Park name | Longitude [deg] | Latitude [deg] | \n",
"
\n",
"\n",
"| Dartmoor national park | \t-3.904 | \t50.58 | \n",
"
\n",
"\n",
"| New forest national park | \t-1.595 | \t50.86 | \n",
"
\n",
"\n",
"| Exmoor national park | \t-3.651 | \t51.14 | \n",
"
\n",
"\n",
"| Pembrokeshire coast national park | \t-4.694 | \t51.64 | \n",
"
\n",
"\n",
"| Brecon beacons national park\t | -3.432 | \t51.88 | \n",
"
\n",
"\n",
"| Pembrokeshire coast national park | \t-4.79 | \t51.99 | \n",
"
\n",
"\n",
"| Norfolk and suffolk broads\t | 1.569 | \t52.62 | \n",
"
\n",
"\n",
"| Snowdonia national park | \t-3.898\t | 52.9 | \n",
"
\n",
"\n",
"| Peak district national park\t | -1.802\t | 53.3 | \n",
"
\n",
"\n",
"| Yorkshire dales national park | \t-2.157\t | 54.23 | \n",
"
\n",
"\n",
"| North yorkshire moors national park\t | -0.8855\t | 54.37 | \n",
"
\n",
"\n",
"| Lake district national park\t | -3.084\t | 54.47 | \n",
"
\n",
"\n",
"| Galloway forest park\t | -4.171\t | 54.87 | \n",
"
\n",
"\n",
"| Northumberland national park\t | -2.228\t | 55.28 | \n",
"
\n",
"\n",
"| Loch lomond and the trossachs national park\t | -4.593\t | 56.24 | \n",
"
\n",
"\n",
"| Tay forest park\t | -4.025 | \t56.59 | \n",
"
\n",
"\n",
"| Cairngorms national park\t | -3.545\t | 57.08 | \n",
"
\n",
"
\n",
"
\n",
"Now try modifying that method so that you can\n",
"\n",
"1. Select the output data type diffrent to Float32\n",
"2. Provide a given projection and geotransform (e.g. if you don't have a GDAL filename)\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reprojecting things"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Previously, we have used the [MODLAND grid converter](http://landweb.nascom.nasa.gov/cgi-bin/developer/tilemap.cgi) site to go from latitude/longitude pairs to MODIS projection. However, in practice, we might want to use a range of different projections, and convert many points at the same time, so how do we go about that?\n",
"\n",
"In GDAL/OGR, most projection-related tools are in the `osr` package, which needs to be imported like e.g. `gdal` itself. The main tools are the `osr.SpatialReference` object, which defines a projection object (with no projection to start with), and the `osr.CoordinateTransformation` object. \n",
"\n",
"Once you instantiate `osr.SpatialReference`, it holds no projection information. You need to use methods to set it up, using EPSG codes, Proj4 strings, or whatever. These methods typically start by `ImportFrom` (e.g. `ImportFromEPSG`, `ImportFromProj4`, ...).\n",
"\n",
"The `CoordinateTransformation` requires a source and destination spatial references that have been configured. Once this is done, it expose the method `TransformPoint` to convert coordinates from the source to the destination projection.\n",
"\n",
"Let's see how this works by converting some latitude/longitude pairs to the Ordnance Survey's [National Grid](http://en.wikipedia.org/wiki/Ordnance_Survey_National_Grid) projection. The projection is also available in [spatialreference.org](http://spatialreference.org/ref/epsg/27700/), where we can gleam its EPSG code (27700). The EPSG code for longitude latitude is [4326](http://spatialreference.org/ref/epsg/4326/). Let's see this in practice:\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-30T18:03:33.081894Z",
"start_time": "2017-10-30T18:03:33.047183Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Site: Snowdonia national park Lon/Lat: -3.898,52.9, OS x,y,z: 272430,335305,-51.8129\n"
]
}
],
"source": [
"from osgeo import osr, ogr\n",
"\n",
"# Define the source projection, WGS84 lat/lon. \n",
"wgs84 = osr.SpatialReference( ) # Define a SpatialReference object\n",
"wgs84.ImportFromEPSG( 4326 ) # And set it to WGS84 using the EPSG code\n",
"\n",
"# Now for the target projection, Ordnance Survey's British National Grid\n",
"osng = osr.SpatialReference() # define the SpatialReference object\n",
"# In this case, we get the projection from a Proj4 string\n",
"osng.ImportFromEPSG( 27700)\n",
"# or, if using the proj4 representation\n",
"#osng.ImportFromProj4 ( \"+proj=tmerc +lat_0=49 +lon_0=-2 \" + \\\n",
"# \"+k=0.9996012717 +x_0=400000 +y_0=-100000 \" + \\\n",
"# \"+ellps=airy +datum=OSGB36 +units=m +no_defs\" )\n",
"\n",
"\n",
"# Now, we define a coordinate transformtion object, *from* wgs84 *to* OSNG\n",
"tx = osr.CoordinateTransformation( wgs84, osng)\n",
"# We loop over the lines of park_data, \n",
"# using the split method to split by newline characters\n",
"park_name, lon, lat = \"Snowdonia national park\", -3.898,52.9\n",
"\n",
"# create a geometry from coordinates\n",
"point = ogr.Geometry(ogr.wkbPoint)\n",
"point.AddPoint(lon, lat)\n",
"\n",
"# Actually do the transformation using the TransformPoint method\n",
"point.Transform ( tx )\n",
"\n",
"osng_x = point.GetX()\n",
"osng_y = point.GetY()\n",
"osng_z = point.GetZ()\n",
"\n",
"# Print out\n",
"print \"Site: {:s} Lon/Lat: {:g},{:g}, OS x,y,z: {:g},{:g},{:g}\".format(park_name, \n",
" lon, lat, osng_x, osng_y, osng_z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can test the result is reasonable by feeding the data for `osng_x` and `osng_y` in the OS own [coordinate conversion website](http://www.ordnancesurvey.co.uk/gps/transformation) and making sure that the calculated longitude latitude pair is the same as the one you started with."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reprojecting whole rasters"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using the Python bindings\n",
"\n",
"We can reproject one raster file from one projection to another easily by using the gdal.Warp method. Let's suppose we want to convert our original `lc_h17v03.img` dataset (which is MODIS sinusoidal projection) into Lat/Long projection (WGS84, EPSG code 4326).\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-30T18:03:34.646322Z",
"start_time": "2017-10-30T18:03:34.539674Z"
}
},
"outputs": [],
"source": [
"gg = gdal.Warp(\"lc_h17v03.img\", \"lc_h17v03_wgs84.tif\", \n",
" format=\"GTiff\", width=2400, height=2400, dstSRS=\"EPSG:4326\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are lots of options you can pass on to `gdal.Warp`. You can see all of them by reading the help associated with the `gdal.WarpOptions` method:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-31T12:05:14.586934Z",
"start_time": "2017-10-31T12:05:14.294597Z"
}
},
"outputs": [],
"source": [
"import gdal\n",
"gdal.WarpOptions?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Using command line tools"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The easiest way to reproject a raster file is to use GDAL's [`gdalwarp`](http://www.gdal.org/gdalwarp.html) tool. As an example, let's say we want to reproject the landcover file from earlier on into latitude/longitude (WGS84):"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-30T18:03:36.059884Z",
"start_time": "2017-10-30T18:03:35.968033Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR 1: Output dataset lc_h17v03_wgs84.tif exists,\n",
"but some command line options were provided indicating a new dataset\n",
"should be created. Please delete existing dataset and run again.\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"# in case you don't have library path set\n",
"# use 'locate libnetcdf` or similar if its not in here\n",
"LD_LIBRARY_PATH=/opt/anaconda/lib:$LD_LIBRARY_PATH\n",
"export LD_LIBRARY_PATH\n",
"\n",
"gdalwarp -of GTiff -t_srs \"EPSG:4326\" -ts 2400 2400 test_lc_h17v03.img lc_h17v03_wgs84.tif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see here that the command takes a number of arguments:\n",
"\n",
"1. `-of GTiff` is the outut format (in this case GeoTIFF)\n",
"2. `-t_srs \"EPSG:4326\"` is the **to** projection (the **from** projection is already specified in the source dataset), in this case, lat/long WGS84, known by its [EPSG code](http://spatialreference.org/ref/epsg/4326/)\n",
"3. `-ts 2400 2400` instructs `gdalwarp` to use an output of size `2400*2400`.\n",
"3. `test_lc_h17v03.img` is the **input dataset**\n",
"4. `lc_h17v03_wgs84.tif` is the **output dataset**\n",
"\n",
"Note that `gdalwarp` will reproject the data, and decide on the pixel size based on some considerations. This can result in the size of the raster changing. If you wanted to still keep the same raster size, we use the `-ts 2400 2400` option, or select an appropriate pixel size using `-tr xres yres` (note it has to be in the target projection, so degrees in this case). We can use `gdalinfo` to see what we've done."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-30T18:03:36.532357Z",
"start_time": "2017-10-30T18:03:36.469958Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"/home/ucfajlg/python/geogg122-1/Chapter4_GDAL\n",
"Driver: HFA/Erdas Imagine Images (.img)\n",
"Files: test_lc_h17v03.img\n",
"Size is 2400, 2400\n",
"Coordinate System is:\n",
"PROJCS[\"Sinusoidal\",\n",
" GEOGCS[\"GCS_Unknown datum based upon the custom spheroid\",\n",
" DATUM[\"Not_specified_based_on_custom_spheroid\",\n",
" SPHEROID[\"Custom_spheroid\",6371007.181,0],\n",
" TOWGS84[0,0,0,-0,-0,-0,0]],\n",
" PRIMEM[\"Greenwich\",0],\n",
" UNIT[\"Degree\",0.017453292519943295]],\n",
" PROJECTION[\"Sinusoidal\"],\n",
" PARAMETER[\"longitude_of_center\",0],\n",
" PARAMETER[\"false_easting\",0],\n",
" PARAMETER[\"false_northing\",0],\n",
" UNIT[\"Meter\",1]]\n",
"Origin = (-1111950.519667000044137,6671703.117999999783933)\n",
"Pixel Size = (463.312716527916677,-463.312716527916507)\n",
"Corner Coordinates:\n",
"Upper Left (-1111950.520, 6671703.118) ( 20d 0' 0.00\"W, 60d 0' 0.00\"N)\n",
"Lower Left (-1111950.520, 5559752.598) ( 15d33'26.06\"W, 50d 0' 0.00\"N)\n",
"Upper Right ( 0.000, 6671703.118) ( 0d 0' 0.01\"E, 60d 0' 0.00\"N)\n",
"Lower Right ( 0.000, 5559752.598) ( 0d 0' 0.01\"E, 50d 0' 0.00\"N)\n",
"Center ( -555975.260, 6115727.858) ( 8d43' 2.04\"W, 55d 0' 0.00\"N)\n",
"Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined\n",
" Description = Layer_1\n",
" Metadata:\n",
" LAYER_TYPE=athematic\n",
"\n"
]
}
],
"source": [
"%%bash\n",
"# in case you don't have library path set\n",
"# use 'locate libnetcdf` or similar if its not in here\n",
"LD_LIBRARY_PATH=/opt/anaconda/lib:$LD_LIBRARY_PATH\n",
"export LD_LIBRARY_PATH\n",
"\n",
"pwd\n",
"gdalinfo test_lc_h17v03.img"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-30T18:03:36.594672Z",
"start_time": "2017-10-30T18:03:36.534446Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Driver: GTiff/GeoTIFF\n",
"Files: lc_h17v03_wgs84.tif\n",
"Size is 2400, 2400\n",
"Coordinate System is:\n",
"GEOGCS[\"WGS 84\",\n",
" DATUM[\"WGS_1984\",\n",
" SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
" AUTHORITY[\"EPSG\",\"7030\"]],\n",
" AUTHORITY[\"EPSG\",\"6326\"]],\n",
" PRIMEM[\"Greenwich\",0],\n",
" UNIT[\"degree\",0.0174532925199433],\n",
" AUTHORITY[\"EPSG\",\"4326\"]]\n",
"Origin = (-19.999999994952233,59.999999994611805)\n",
"Pixel Size = (0.008333919248404,-0.004166959624202)\n",
"Metadata:\n",
" AREA_OR_POINT=Area\n",
"Image Structure Metadata:\n",
" INTERLEAVE=BAND\n",
"Corner Coordinates:\n",
"Upper Left ( -20.0000000, 60.0000000) ( 20d 0' 0.00\"W, 60d 0' 0.00\"N)\n",
"Lower Left ( -20.0000000, 49.9992969) ( 20d 0' 0.00\"W, 49d59'57.47\"N)\n",
"Upper Right ( 0.0014062, 60.0000000) ( 0d 0' 5.06\"E, 60d 0' 0.00\"N)\n",
"Lower Right ( 0.0014062, 49.9992969) ( 0d 0' 5.06\"E, 49d59'57.47\"N)\n",
"Center ( -9.9992969, 54.9996484) ( 9d59'57.47\"W, 54d59'58.73\"N)\n",
"Band 1 Block=2400x3 Type=Byte, ColorInterp=Gray\n",
" Description = Layer_1\n",
" Metadata:\n",
" LAYER_TYPE=athematic\n"
]
}
],
"source": [
"%%bash\n",
"# in case you don't have library path set\n",
"# use 'locate libnetcdf` or similar if its not in here\n",
"LD_LIBRARY_PATH=/opt/anaconda/lib:$LD_LIBRARY_PATH\n",
"export LD_LIBRARY_PATH\n",
"\n",
"gdalinfo lc_h17v03_wgs84.tif"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's see how different these two datasets are:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"ExecuteTime": {
"end_time": "2017-10-30T18:03:38.917157Z",
"start_time": "2017-10-30T18:03:37.021890Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"