In the previous session, we used the GDAL library to open HDF files. GDAL is not limited to a single file format, but can actually cope with many different raster file formats semalessly. For vector data (i.e., data that is stored by features, each being made up of fields containing different types of information, one of them being a geometry, such as polygon, line or point), GDAL has a sister library called OGR. The usefulness of these two libraries is that they allow the user to deal with many of the different file formats in a consistent way.
It is important to note that both GDAL and OGR come with a suite of command line tools that let you do many complex tasks on the command line directly. A listing of the GDAL command line tools is available here, but bear in mind that many blogs etc. carry out examples of using GDAL tools in practice. For OGR, most of the library can be accessed using ogr2ogr, but as usual, you might find more useful information on blogs etc.
GDAL provides a very handy way of dealing with raster data in many different formats, not only by making access to the data easy, but also abstracting the nuances and complications of different file formats. In GDAL, a raster file is made up of the actual raster data (i.e., the values of each pixel of LAI that we saw before), and of some metadata. Metadata is data that describes something, and in this case it could be the projection, the location of corners and pixel spacing, etc.
GDAL stores information about the location of each pixel using the GeoTransform. The GeoTransform contains the coordinates (in some projection) of the upper left (UL) corner of the image (taken to be the borders of the pixel in the UL corner, not the center), the pixel spacing and an additional rotation. By knowing these figures, we can calculate the location of each pixel in the image easily. Let's see how this works with an easy example. We have prepared a GeoTIFF file (GeoTIFF is the more ubiquitous file format for EO data) of the MODIS landcover product for the UK. The data has been extracted from the HDF-EOS files that are similar to the LAI product that we have seen before, and converted. The file is lc_h17v03.tif
. We will open the file in Python, and have a look at finding a particular location.
Assume we are interested in locating Kinder Scout, a moorland in the Peak District National Park. Its coordinates are 1.871417W, 53.384726N. In the MODIS integerised sinusoidal projection, the coordinates are (-124114.3, 5936117.4) (you can use the MODLAND tile calculator website to do this calculation yourself).
import gdal # Import GDAL library
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# have to make sure have access to gdal data files
import os
if 'GDAL_DATA' not in os.environ:
os.environ["GDAL_DATA"] = '/opt/anaconda/share/gdal'
g = gdal.Open ( "lc_h17v03.tif" ) # Open the file
if g is None:
print "Could not open the file!"
geo_transform = g.GetGeoTransform ()
print geo_transform
print g.RasterXSize, g.RasterYSize
(-1111950.519667, 463.3127165279167, 0.0, 6671703.118, 0.0, -463.3127165279165) 2400 2400
In the previous code, we open the file (we just use the filename), and then query the object for its GeoTransform, which we then print out. The 6-element tuple comprises
We have also seen that the dataset is of size 2400x2400, using RasterXSize
and RasterYSize
. The goal is to find the pixel number $(i,j), \;\;0\le i,j < 2400$ that corresponds to Kinder Scout. To do this, we use the following calculations:
pixel_x = (-124114.3 - geo_transform[0])/geo_transform[1] \
# The difference in distance between the UL corner (geot[0] \
#and point of interest. Scaled by geot[1] to get pixel number
pixel_y = (5936117.4 - geo_transform[3])/(geo_transform[5]) # Like for pixel_x, \
#but in vertical direction. Note the different elements of geot \
#being used
print pixel_x, pixel_y
2132.11549009 1587.66572071
If you can't be bothered doing those simple calculations, GDAL offers you a couple of methods to do the same thing already: you can use ApplyGeoTransform
, a method that takes:
Using the standard geo-transform, you can go from pixel location to projected units:
print gdal.ApplyGeoTransform(geo_transform, pixel_x, pixel_y)
[-124114.30000000016, 5936117.4]
If we want to go the other way (from projected units to pixel row/column pair) we can "invert" the geotransform as follows:
print gdal.ApplyGeoTransform(gdal.InvGeoTransform(geo_transform), -124114.3, 5936117.4)
[2132.115490094644, 1587.6657207091303]
So the pixel number is a floating point number, which we might need to round off as an integer. Let's plot the entire raster map (with minimum value 0 to ignore the ocean) using plt.imshow
and plot the location of Kinder Scout using plt.plot
. We will also use plt.annotate
to add a label with an arrow:
lc = g.ReadAsArray() # Read raster data
# Now plot the raster data using gist_earth palette
plt.figure(figsize=(11,11))
plt.imshow ( lc, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth )
# Plot location of our area of interest as a red dot ('ro')
plt.plot ( pixel_x, pixel_y, 'ro')
# Annotate
plt.annotate('Kinder Scout', xy=(pixel_x, pixel_y), \
xycoords='data', xytext=(-120, -40), \
textcoords='offset points', size=12, \
bbox=dict(boxstyle="round4,pad=.5", fc="0.8"), \
arrowprops=dict(arrowstyle="->", \
connectionstyle="angle,angleA=0,angleB=-90,rad=10", \
color='w'), )
# Remove vertical and horizontal ticks
plt.xticks([])
plt.yticks([])
([], <a list of 0 Text yticklabel objects>)
Park name | Longitude [deg] | Latitude [deg] |
---|---|---|
Dartmoor national park | -3.904 | 50.58 |
New forest national park | -1.595 | 50.86 |
Exmoor national park | -3.651 | 51.14 |
Pembrokeshire coast national park | -4.694 | 51.64 |
Brecon beacons national park | -3.432 | 51.88 |
Pembrokeshire coast national park | -4.79 | 51.99 |
Norfolk and suffolk broads | 1.569 | 52.62 |
Snowdonia national park | -3.898 | 52.9 |
Peak district national park | -1.802 | 53.3 |
Yorkshire dales national park | -2.157 | 54.23 |
North yorkshire moors national park | -0.8855 | 54.37 |
Lake district national park | -3.084 | 54.47 |
Galloway forest park | -4.171 | 54.87 |
Northumberland national park | -2.228 | 55.28 |
Loch lomond and the trossachs national park | -4.593 | 56.24 |
Tay forest park | -4.025 | 56.59 |
Cairngorms national park | -3.545 | 57.08 |
Projections in GDAL objects are stored can be accessed by querying the dataset using the GetProjection()
method. If we do that on the currently opened dataset (stored in variable g
), we get:
print g.GetProjection()
PROJCS["unnamed",GEOGCS["Unknown datum based upon the custom spheroid",DATUM["Not_specified_based_on_custom_spheroid",SPHEROID["Custom spheroid",6371007.181,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Sinusoidal"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
The above is the description of the projection (in this case, MODIS sinusoidal) in WKT
(well-known text) format. There are a number of different ways of specifying projections, the most important being
The site spatialreference.org allows you to search a large collection of projections, and get the representation that you want to use.
So far, we have read data from files, but lets see how we can save raster data to a new file. We will use the previous landcover map as an example. We will write a method to save the data in a format provided by the user. The procedure is fairly straightforward: we get a handler to a driver (e.g. a GeoTIFF or Erdas Imagine format), we create the output file (giving a filename, number of rows, columns, bands, the data type), and then add the relevant metadata (projection, geotransform, ...). We then select a band from the output and copy the array that we want to write to that band.
g = gdal.Open ( "lc_h17v03.tif" ) # Open original file
# Get the x, y and number of bands from the original file
x_size, y_size, n_bands = g.RasterXSize, g.RasterYSize, g.RasterCount
data = g.ReadAsArray ()
driver = gdal.GetDriverByName ( "HFA" ) # Get a handler to a driver
# Maybe try "GeoTIFF" here
# Next line creates the output dataset with
# 1. The filename ("test_lc_h17v03.img")
# 2. The raster size (x_size, y_size)
# 3. The number of bands
# 4.The data type (in this case, Byte.
# Other typical values might be gdal.GDT_Int16
# or gdal.GDT_Float32)
dataset_out = driver.Create ( "test_lc_h17v03.img", x_size, y_size, n_bands,
gdal.GDT_Byte )
# Set the output geotransform by reading the input one
dataset_out.SetGeoTransform ( g.GetGeoTransform() )
# Set the output projection by reading the input one
dataset_out.SetProjection ( g.GetProjectionRef() )
# Now, get band # 1, and write our data array.
# Note that the data array needs to have the same type
# as the one specified for dataset_out
dataset_out.GetRasterBand ( 1 ).WriteArray ( data )
# This bit forces GDAL to close the file and write to it
dataset_out = None
The output file should hopefully exist in this directory. Let's use gdalinfo
to find out about it
!gdalinfo test_lc_h17v03.img
Driver: HFA/Erdas Imagine Images (.img) Files: test_lc_h17v03.img Size is 2400, 2400 Coordinate System is: PROJCS["Sinusoidal", GEOGCS["GCS_Unknown datum based upon the custom spheroid", DATUM["Not_specified_based_on_custom_spheroid", SPHEROID["Custom_spheroid",6371007.181,0], TOWGS84[0,0,0,-0,-0,-0,0]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["Sinusoidal"], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["Meter",1]] Origin = (-1111950.519667000044137,6671703.117999999783933) Pixel Size = (463.312716527916677,-463.312716527916507) Corner Coordinates: Upper Left (-1111950.520, 6671703.118) ( 20d 0' 0.00"W, 60d 0' 0.00"N) Lower Left (-1111950.520, 5559752.598) ( 15d33'26.06"W, 50d 0' 0.00"N) Upper Right ( 0.000, 6671703.118) ( 0d 0' 0.01"E, 60d 0' 0.00"N) Lower Right ( 0.000, 5559752.598) ( 0d 0' 0.01"E, 50d 0' 0.00"N) Center ( -555975.260, 6115727.858) ( 8d43' 2.04"W, 55d 0' 0.00"N) Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined Description = Layer_1 Metadata: LAYER_TYPE=athematic
So the previous code works. Since this is something we typically do (read some data from one or more files, manipulate it and save the result in output files), it makes a lot of sense to try to put this code in a method that is more or less generic, that we can test and then re-use. Here's a first attempt at it:
def save_raster ( output_name, raster_data, dataset, driver="GTiff" ):
"""
A function to save a 1-band raster using GDAL to the file indicated
by ``output_name``. It requires a GDAL-accesible dataset to collect
the projection and geotransform.
Parameters
----------
output_name: str
The output filename, with full path and extension if required
raster_data: array
The array that we want to save
dataset: str
Filename of a GDAL-friendly dataset that we want to use to
read geotransform & projection information
driver: str
A GDAL driver string, like GTiff or HFA.
"""
# Open the reference dataset
g = gdal.Open ( dataset )
# Get the Geotransform vector
geo_transform = g.GetGeoTransform ()
x_size = g.RasterXSize # Raster xsize
y_size = g.RasterYSize # Raster ysize
srs = g.GetProjectionRef () # Projection
# Need a driver object. By default, we use GeoTIFF
driver = gdal.GetDriverByName ( driver )
dataset_out = driver.Create ( output_name, x_size, y_size, 1, \
gdal.GDT_Float32 )
dataset_out.SetGeoTransform ( geo_transform )
dataset_out.SetProjection ( srs )
dataset_out.GetRasterBand ( 1 ).WriteArray ( \
raster_data.astype(np.float32) )
dataset_out = None
Previously, we have used the MODLAND grid converter 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?
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.
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
, ...).
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.
Let's see how this works by converting some latitude/longitude pairs to the Ordnance Survey's National Grid projection. The projection is also available in spatialreference.org, where we can gleam its EPSG code (27700). The EPSG code for longitude latitude is 4326. Let's see this in practice:
from osgeo import osr, ogr
# Define the source projection, WGS84 lat/lon.
wgs84 = osr.SpatialReference( ) # Define a SpatialReference object
wgs84.ImportFromEPSG( 4326 ) # And set it to WGS84 using the EPSG code
# Now for the target projection, Ordnance Survey's British National Grid
osng = osr.SpatialReference() # define the SpatialReference object
# In this case, we get the projection from a Proj4 string
osng.ImportFromEPSG( 27700)
# or, if using the proj4 representation
#osng.ImportFromProj4 ( "+proj=tmerc +lat_0=49 +lon_0=-2 " + \
# "+k=0.9996012717 +x_0=400000 +y_0=-100000 " + \
# "+ellps=airy +datum=OSGB36 +units=m +no_defs" )
# Now, we define a coordinate transformtion object, *from* wgs84 *to* OSNG
tx = osr.CoordinateTransformation( wgs84, osng)
# We loop over the lines of park_data,
# using the split method to split by newline characters
park_name, lon, lat = "Snowdonia national park", -3.898,52.9
# create a geometry from coordinates
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(lon, lat)
# Actually do the transformation using the TransformPoint method
point.Transform ( tx )
osng_x = point.GetX()
osng_y = point.GetY()
osng_z = point.GetZ()
# Print out
print "Site: {:s} Lon/Lat: {:g},{:g}, OS x,y,z: {:g},{:g},{:g}".format(park_name,
lon, lat, osng_x, osng_y, osng_z)
Site: Snowdonia national park Lon/Lat: -3.898,52.9, OS x,y,z: 272430,335305,-51.8129
You can test the result is reasonable by feeding the data for osng_x
and osng_y
in the OS own coordinate conversion website and making sure that the calculated longitude latitude pair is the same as the one you started with.
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).
gg = gdal.Warp("lc_h17v03.img", "lc_h17v03_wgs84.tif",
format="GTiff", width=2400, height=2400, dstSRS="EPSG:4326")
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:
import gdal
gdal.WarpOptions?
The easiest way to reproject a raster file is to use GDAL's gdalwarp
tool. As an example, let's say we want to reproject the landcover file from earlier on into latitude/longitude (WGS84):
%%bash
# in case you don't have library path set
# use 'locate libnetcdf` or similar if its not in here
LD_LIBRARY_PATH=/opt/anaconda/lib:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH
gdalwarp -of GTiff -t_srs "EPSG:4326" -ts 2400 2400 test_lc_h17v03.img lc_h17v03_wgs84.tif
ERROR 1: Output dataset lc_h17v03_wgs84.tif exists, but some command line options were provided indicating a new dataset should be created. Please delete existing dataset and run again.
We see here that the command takes a number of arguments:
-of GTiff
is the outut format (in this case GeoTIFF)-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-ts 2400 2400
instructs gdalwarp
to use an output of size 2400*2400
.test_lc_h17v03.img
is the input datasetlc_h17v03_wgs84.tif
is the output datasetNote 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.
%%bash
# in case you don't have library path set
# use 'locate libnetcdf` or similar if its not in here
LD_LIBRARY_PATH=/opt/anaconda/lib:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH
pwd
gdalinfo test_lc_h17v03.img
/home/ucfajlg/python/geogg122-1/Chapter4_GDAL Driver: HFA/Erdas Imagine Images (.img) Files: test_lc_h17v03.img Size is 2400, 2400 Coordinate System is: PROJCS["Sinusoidal", GEOGCS["GCS_Unknown datum based upon the custom spheroid", DATUM["Not_specified_based_on_custom_spheroid", SPHEROID["Custom_spheroid",6371007.181,0], TOWGS84[0,0,0,-0,-0,-0,0]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["Sinusoidal"], PARAMETER["longitude_of_center",0], PARAMETER["false_easting",0], PARAMETER["false_northing",0], UNIT["Meter",1]] Origin = (-1111950.519667000044137,6671703.117999999783933) Pixel Size = (463.312716527916677,-463.312716527916507) Corner Coordinates: Upper Left (-1111950.520, 6671703.118) ( 20d 0' 0.00"W, 60d 0' 0.00"N) Lower Left (-1111950.520, 5559752.598) ( 15d33'26.06"W, 50d 0' 0.00"N) Upper Right ( 0.000, 6671703.118) ( 0d 0' 0.01"E, 60d 0' 0.00"N) Lower Right ( 0.000, 5559752.598) ( 0d 0' 0.01"E, 50d 0' 0.00"N) Center ( -555975.260, 6115727.858) ( 8d43' 2.04"W, 55d 0' 0.00"N) Band 1 Block=64x64 Type=Byte, ColorInterp=Undefined Description = Layer_1 Metadata: LAYER_TYPE=athematic
%%bash
# in case you don't have library path set
# use 'locate libnetcdf` or similar if its not in here
LD_LIBRARY_PATH=/opt/anaconda/lib:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH
gdalinfo lc_h17v03_wgs84.tif
Driver: GTiff/GeoTIFF Files: lc_h17v03_wgs84.tif Size is 2400, 2400 Coordinate System is: GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433], AUTHORITY["EPSG","4326"]] Origin = (-19.999999994952233,59.999999994611805) Pixel Size = (0.008333919248404,-0.004166959624202) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( -20.0000000, 60.0000000) ( 20d 0' 0.00"W, 60d 0' 0.00"N) Lower Left ( -20.0000000, 49.9992969) ( 20d 0' 0.00"W, 49d59'57.47"N) Upper Right ( 0.0014062, 60.0000000) ( 0d 0' 5.06"E, 60d 0' 0.00"N) Lower Right ( 0.0014062, 49.9992969) ( 0d 0' 5.06"E, 49d59'57.47"N) Center ( -9.9992969, 54.9996484) ( 9d59'57.47"W, 54d59'58.73"N) Band 1 Block=2400x3 Type=Byte, ColorInterp=Gray Description = Layer_1 Metadata: LAYER_TYPE=athematic
Let's see how different these two datasets are:
g = gdal.Open ( "lc_h17v03_wgs84.tif" )
wgs84 = g.ReadAsArray()
g = gdal.Open("test_lc_h17v03.img")
modis = g.ReadAsArray()
fig, axs = plt.subplots(nrows=1, ncols=2,
sharex=True, sharey=True, figsize=(13, 5))
axs[0].imshow ( modis, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth )
axs[1].imshow ( wgs84, interpolation='nearest', vmin=0, cmap=plt.cm.gist_earth )
<matplotlib.image.AxesImage at 0x7f1b56b0e410>
Suppose in the previous example that one is only interested in Ireland. How could you ignore the pixels that don't believe to the UK? One way to go around that would be if we had a vector with the country boundaries, then we could select the pixels that fall within Ireland, and create a binary mask of the same size as our map with True
within Ireland and False
otherwise.
Predictably, GDAL can do this. We will use a global map of countries provided as a zipped shapefile. You can download that file using requests as here:
import requests
r = requests.get ("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/" +
"cultural/ne_110m_admin_0_countries.zip")
if r.ok:
with open("data/ne_110m_admin_0_countries.zip", 'w') as fp:
fp.write(r.content)
But if you can't be bothered, you can also open remote zipped files with GDAL. Given this is a zipfile, we can use the command line tool ogrinfo
to see the contents of the file and, in particular, the layers it contains. Basically, we build the filename as
/vsizip/vsicurl/<remote_filename.zip>/<file_inside_zipfile>
We get a few warnings, but...
!ogrinfo -al -so \
/vsizip/vsicurl/http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp
ERROR 1: Random access not supported for /vsizip Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". ERROR 1: Random access not supported for /vsizip Warning 4: Failed to open /vsizip/vsicurl/http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp, No such file or directory. Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Had to open data source read-only. INFO: Open of `/vsizip/vsicurl/http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp' using driver `ESRI Shapefile' successful. Layer name: ne_110m_admin_0_countries Metadata: DBF_DATE_LAST_UPDATE=2017-10-04 Geometry: Polygon Feature Count: 177 Extent: (-180.000000, -90.000000) - (180.000000, 83.645130) Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Layer SRS WKT: GEOGCS["WGS84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["degree",0.017453292519943295], AUTHORITY["EPSG","4326"]] scalerank: Integer (4.0) featurecla: String (30.0) LABELRANK: Real (6.2) SOVEREIGNT: String (32.0) SOV_A3: String (3.0) ADM0_DIF: Real (4.2) LEVEL: Real (4.2) TYPE: String (17.0) ADMIN: String (36.0) ADM0_A3: String (3.0) GEOU_DIF: Real (4.2) GEOUNIT: String (36.0) GU_A3: String (3.0) SU_DIF: Real (4.2) SUBUNIT: String (36.0) SU_A3: String (3.0) BRK_DIFF: Real (4.2) NAME: String (36.0) NAME_LONG: String (36.0) BRK_A3: String (3.0) BRK_NAME: String (36.0) BRK_GROUP: String (30.0) ABBREV: String (13.0) POSTAL: String (4.0) FORMAL_EN: String (52.0) FORMAL_FR: String (35.0) NAME_CIAWF: String (45.0) NOTE_ADM0: String (22.0) NOTE_BRK: String (164.0) NAME_SORT: String (36.0) NAME_ALT: String (38.0) MAPCOLOR7: Real (4.2) MAPCOLOR8: Real (4.2) MAPCOLOR9: Real (4.2) MAPCOLOR13: Real (6.2) POP_EST: Real (13.2) POP_RANK: Real (5.2) GDP_MD_EST: Real (11.2) POP_YEAR: Real (7.2) LASTCENSUS: Real (7.2) GDP_YEAR: Real (7.2) ECONOMY: String (26.0) INCOME_GRP: String (23.0) WIKIPEDIA: Real (6.2) FIPS_10_: String (3.0) ISO_A2: String (5.0) ISO_A3: String (3.0) ISO_A3_EH: String (3.0) ISO_N3: String (3.0) UN_A3: String (4.0) WB_A2: String (3.0) WB_A3: String (3.0) WOE_ID: Real (11.2) WOE_ID_EH: Real (11.2) WOE_NOTE: String (190.0) ADM0_A3_IS: String (3.0) ADM0_A3_US: String (3.0) ADM0_A3_UN: Real (6.2) ADM0_A3_WB: Real (6.2) CONTINENT: String (23.0) REGION_UN: String (23.0) SUBREGION: String (25.0) REGION_WB: String (26.0) NAME_LEN: Real (5.2) LONG_LEN: Real (5.2) ABBREV_LEN: Real (5.2) TINY: Real (6.2) HOMEPART: Real (6.2) MIN_ZOOM: Real (6.2) MIN_LABEL: Real (6.2) MAX_LABEL: Real (6.2)
We see that we can basically filter when the NAME
attribute is equal to the country name. On the shell, we can see what information is available on ireland by filtering:
!ogrinfo -al -where "NAME='IRELAND'" \
/vsizip/vsicurl/http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp
ERROR 1: Random access not supported for /vsizip Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". ERROR 1: Random access not supported for /vsizip Warning 4: Failed to open /vsizip/vsicurl/http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp, No such file or directory. Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Had to open data source read-only. INFO: Open of `/vsizip/vsicurl/http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp' using driver `ESRI Shapefile' successful. Layer name: ne_110m_admin_0_countries Metadata: DBF_DATE_LAST_UPDATE=2017-10-04 Geometry: Polygon Feature Count: 1 Extent: (-180.000000, -90.000000) - (180.000000, 83.645130) Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Warning 1: Recode from CP437 to UTF-8 failed with the error: "Invalid argument". Layer SRS WKT: GEOGCS["WGS84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["degree",0.017453292519943295], AUTHORITY["EPSG","4326"]] scalerank: Integer (4.0) featurecla: String (30.0) LABELRANK: Real (6.2) SOVEREIGNT: String (32.0) SOV_A3: String (3.0) ADM0_DIF: Real (4.2) LEVEL: Real (4.2) TYPE: String (17.0) ADMIN: String (36.0) ADM0_A3: String (3.0) GEOU_DIF: Real (4.2) GEOUNIT: String (36.0) GU_A3: String (3.0) SU_DIF: Real (4.2) SUBUNIT: String (36.0) SU_A3: String (3.0) BRK_DIFF: Real (4.2) NAME: String (36.0) NAME_LONG: String (36.0) BRK_A3: String (3.0) BRK_NAME: String (36.0) BRK_GROUP: String (30.0) ABBREV: String (13.0) POSTAL: String (4.0) FORMAL_EN: String (52.0) FORMAL_FR: String (35.0) NAME_CIAWF: String (45.0) NOTE_ADM0: String (22.0) NOTE_BRK: String (164.0) NAME_SORT: String (36.0) NAME_ALT: String (38.0) MAPCOLOR7: Real (4.2) MAPCOLOR8: Real (4.2) MAPCOLOR9: Real (4.2) MAPCOLOR13: Real (6.2) POP_EST: Real (13.2) POP_RANK: Real (5.2) GDP_MD_EST: Real (11.2) POP_YEAR: Real (7.2) LASTCENSUS: Real (7.2) GDP_YEAR: Real (7.2) ECONOMY: String (26.0) INCOME_GRP: String (23.0) WIKIPEDIA: Real (6.2) FIPS_10_: String (3.0) ISO_A2: String (5.0) ISO_A3: String (3.0) ISO_A3_EH: String (3.0) ISO_N3: String (3.0) UN_A3: String (4.0) WB_A2: String (3.0) WB_A3: String (3.0) WOE_ID: Real (11.2) WOE_ID_EH: Real (11.2) WOE_NOTE: String (190.0) ADM0_A3_IS: String (3.0) ADM0_A3_US: String (3.0) ADM0_A3_UN: Real (6.2) ADM0_A3_WB: Real (6.2) CONTINENT: String (23.0) REGION_UN: String (23.0) SUBREGION: String (25.0) REGION_WB: String (26.0) NAME_LEN: Real (5.2) LONG_LEN: Real (5.2) ABBREV_LEN: Real (5.2) TINY: Real (6.2) HOMEPART: Real (6.2) MIN_ZOOM: Real (6.2) MIN_LABEL: Real (6.2) MAX_LABEL: Real (6.2) OGRFeature(ne_110m_admin_0_countries):74 scalerank (Integer) = 1 featurecla (String) = Admin-0 country LABELRANK (Real) = 3.00 SOVEREIGNT (String) = Ireland SOV_A3 (String) = IRL ADM0_DIF (Real) = 0.00 LEVEL (Real) = 2.00 TYPE (String) = Sovereign country ADMIN (String) = Ireland ADM0_A3 (String) = IRL GEOU_DIF (Real) = 0.00 GEOUNIT (String) = Ireland GU_A3 (String) = IRL SU_DIF (Real) = 0.00 SUBUNIT (String) = Ireland SU_A3 (String) = IRL BRK_DIFF (Real) = 0.00 NAME (String) = Ireland NAME_LONG (String) = Ireland BRK_A3 (String) = IRL BRK_NAME (String) = Ireland BRK_GROUP (String) = (null) ABBREV (String) = Ire. POSTAL (String) = IRL FORMAL_EN (String) = Ireland FORMAL_FR (String) = (null) NAME_CIAWF (String) = Ireland NOTE_ADM0 (String) = (null) NOTE_BRK (String) = (null) NAME_SORT (String) = Ireland NAME_ALT (String) = (null) MAPCOLOR7 (Real) = 2.00 MAPCOLOR8 (Real) = 3.00 MAPCOLOR9 (Real) = 2.00 MAPCOLOR13 (Real) = 2.00 POP_EST (Real) = 5011102.00 POP_RANK (Real) = 13.00 GDP_MD_EST (Real) = 322000.00 POP_YEAR (Real) = 2017.00 LASTCENSUS (Real) = 2011.00 GDP_YEAR (Real) = 2016.00 ECONOMY (String) = 2. Developed region: nonG7 INCOME_GRP (String) = 1. High income: OECD WIKIPEDIA (Real) = -99.00 FIPS_10_ (String) = EI ISO_A2 (String) = IE ISO_A3 (String) = IRL ISO_A3_EH (String) = IRL ISO_N3 (String) = 372 UN_A3 (String) = 372 WB_A2 (String) = IE WB_A3 (String) = IRL WOE_ID (Real) = 23424803.00 WOE_ID_EH (Real) = 23424803.00 WOE_NOTE (String) = Exact WOE match as country ADM0_A3_IS (String) = IRL ADM0_A3_US (String) = IRL ADM0_A3_UN (Real) = -99.00 ADM0_A3_WB (Real) = -99.00 CONTINENT (String) = Europe REGION_UN (String) = Europe SUBREGION (String) = Northern Europe REGION_WB (String) = Europe & Central Asia NAME_LEN (Real) = 7.00 LONG_LEN (Real) = 7.00 ABBREV_LEN (Real) = 4.00 TINY (Real) = -99.00 HOMEPART (Real) = 1.00 MIN_ZOOM (Real) = 0.00 MIN_LABEL (Real) = 3.00 MAX_LABEL (Real) = 8.00 POLYGON ((-7.57216793459106 55.1316222194549,-7.36603064617879 54.5958409694527,-7.57216793459106 54.059956366586,-6.95373023113807 54.0737022975756,-6.19788489422103 53.8675650091634,-6.03298539877761 53.1531641709444,-6.78885657391085 52.2601179062923,-8.56161658368356 51.6693012558994,-9.97708574059027 51.8204548203531,-9.16628251793078 52.8646288112427,-9.68852454267245 53.8813626165853,-8.32798743329201 54.6645189479686,-7.57216793459106 55.1316222194549))
The associated polygon is a polygon describing the borders of Ireland. We want to rasterise that polygon to provide a mask for e.g. lc_h17v03_wgs84.tif
(the Lat/Long projection you created earlier). To do this, we need to extract information about
We will extract the relevant information in what follows
g = gdal.Open("lc_h17v03_wgs84.tif")
proj = g.GetProjectionRef()
geoT = g.GetGeoTransform()
xs = []
ys = []
for x,y in [ [0, 0], [0, g.RasterYSize], [g.RasterXSize, g.RasterYSize], [g.RasterXSize, 0]]:
xx, yy = gdal.ApplyGeoTransform(geoT, x,y)
xs.append(xx)
ys.append(yy)
extent = [min(xs), min(ys), max(xs), max(ys)]
xRes = geoT[1]
yRes = geoT[-1]
nx = g.RasterXSize
ny = g.RasterYSize
gg = gdal.Rasterize("data/lc_h17v03_wgs84_ireland.tif", "/vsizip/" +
"data/ne_110m_admin_0_countries.zip/ne_110m_admin_0_countries.shp",
format="GTiff", xRes=xRes, yRes=yRes, where="NAME='Ireland'",
outputBounds=[min(xs), min(ys), max(xs), max(ys)],
width=nx, height=ny, noData=0, burnValues=1)
if gg is not None:
print "Done!"
Done!
mask = gg.ReadAsArray()
g = gdal.Open ( "lc_h17v03_wgs84.tif" )
wgs84 = g.ReadAsArray()
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(13,5))
axs[0].imshow(mask, interpolation="nearest")
axs[1].imshow((1.*mask)*wgs84, interpolation="nearest", cmap=plt.cm.gist_earth)
<matplotlib.image.AxesImage at 0x7f1b50c505d0>
That was using the lowest resolution map. Let's use something a bit better... There's a higher resolution map available from
[]http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries.zip%5D(http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries.zip). We'll download it if it's not available in your data/
folder, and do the same as above..
import os
if not os.path.exists("data/ne_50m_admin_0_countries.zip"):
r = requests.get ("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/"
+ "cultural/ne_50m_admin_0_countries.zip")
if r.ok:
with open("data/ne_50m_admin_0_countries.zip", 'w') as fp:
fp.write(r.content)
g = gdal.Open("lc_h17v03_wgs84.tif")
proj = g.GetProjectionRef()
geoT = g.GetGeoTransform()
xs = []
ys = []
for x,y in [ [0, 0], [0, g.RasterYSize], [g.RasterXSize, g.RasterYSize], [g.RasterXSize, 0]]:
xx, yy = gdal.ApplyGeoTransform(geoT, x,y)
xs.append(xx)
ys.append(yy)
extent = [min(xs), min(ys), max(xs), max(ys)]
xRes = geoT[1]
yRes = geoT[-1]
nx = g.RasterXSize
ny = g.RasterYSize
gg = gdal.Rasterize("data/lc_h17v03_wgs84_ireland.tif", "/vsizip/" +
"data/ne_50m_admin_0_countries.zip/ne_50m_admin_0_countries.shp",
format="GTiff", xRes=xRes, yRes=yRes, where="NAME='Ireland'",
outputBounds=[min(xs), min(ys), max(xs), max(ys)],
width=nx, height=ny, noData=0, burnValues=1)
if gg is not None:
print "Done!"
mask = gg.ReadAsArray()
g = gdal.Open ( "lc_h17v03_wgs84.tif" )
wgs84 = g.ReadAsArray()
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(13,5))
axs[0].imshow(mask, interpolation="nearest")
axs[1].imshow((1.*mask)*wgs84, interpolation="nearest", cmap=plt.cm.gist_earth)
Done!
<matplotlib.image.AxesImage at 0x7f1b50ace690>