An Introduction to GDAL - Robert Simmon

[Music] [Applause] [Music] Thank You Irene it's definitely good to be back here so this is the Canyonlands in southeastern Utah and it is the last place in the United States that was mapped in detail actually a few tough you separate attempts were made to do that mapping in 1869 and 1871 John John Wesley Powell led an expedition down the green in Colorado rivers a scientific expedition to map and explore this area mostly geologists and surveyors but they were limited to traveling along the rivers because the area is so rugged and inhospitable and remote and even when they could clamber off that the river itself and up onto the canyon edges you know this was their view so hoodoos and mesas but they really didn't get a good perspective which you need to build a map so in 1885 the USGS produced this so it was the first of the official US government maps of the desert Southwest and this went through five different editions over 65 years with only the only difference being the typography so the labeling of features and it wasn't until the mid-1950s that we got these very detailed topographic maps that were similar with today and that really required the development of some new technologies and it also was forced on by our uranium boom so prospectors actually crawled over this entire area looking for material to fuel the us's nuclear program and the technologies required were aircraft big heavy aircraft that could carry these huge mapping cameras the cameras themselves that could take wide area high contrast high-resolution photographs and then the science of photogrammetry and what that is is it allows by looking at overlapping photographs it allows you to precisely plot where or each feature is in relation to the other photographs and eventually onto the surface of the earth and from those we have the wonderful high-res USGS one to 24,000 Maps with these great highly detailed topographic contours so that took 65 years in the intervening 65 years we have gone from not being able to map a place at all to being able to see it almost every day so this is a planet image from March 12 and another one from April camp and in addition to this color imagery there's also an amazing amount of other types of data so the entire wealth of US government satellite data plus you know the data exhaust from our phones geologic maps everything so there's dozens of different data sets of every place on earth and they're all useful separately but they become incredibly powerful when you start combining these different types of datasets together but then that becomes a challenge because they're all in different formats they all might be in different map projections so how do you really stitch everything together to make a greater whole and one answer to that is G tall so G tall is the geospatial data abstraction library and it is in fact pronounced G tall not Google maybe back in the old days it was called doodle but then Frank went to work for Google and there was Google Doodle maybe a little confusion um so it's doodle and it's a library that basically allows you to manipulate geospatial data so data that is tied to a specific place on earth and it's great because it's fast it's free and it's flexible so it supports a vast number of data formats it's runs on multiple platforms it's automatable both through bash any command line scripting and then it has powerful Python bindings so you can bring things in and then do analysis and um pie and things like that unfortunately it's really kind of deep and obscure so it's inconsistently documented a lot of times if you read the detail page on it on a certain function it'll be like well the widget widgets makes widgets you know the widget function makes digits and that's not super useful it also requires a thorough knowledge of map projections and data formats and a certain amount of being able to manage like installing things and unix and managing the dependencies and things like that so hopefully I can show why GDL is cool and then give you some of the tools so that you can then start to learn it on your own afterwards and like even know the right terms to Google so I'll just start out with some first steps showing some very basic commands that are essentially the G dolls hello world and then also introduce a few of the more obscure strings that you might find in options then talk about map projections and how they're necessary because the earth is round not flat despite what you might think if you google my name on YouTube then copy covers some functions jido warp some of the more esoteric options for how to control that and then the concept of a virtual file and then go into local map projection so it turns out Earth not only isn't flat it's not perfectly round either so we have all of these mathematical contracts to sort of model that non-spherical nests and so that matters a little bit for local projections and then some techniques for making large scale which means highly detailed maps from either a big data set or many smaller data sets put together and then because the earth isn't static we can also manipulate satellite data with GDL and so we'll talk about reordering and combining bands scaling data and then getting data from different formats to play well together so for installing GDL Konda is a good route if you're already using it be aware though that it installs a slightly older version that won't read JPEG 2000 so you might need to fig use geo forge to get GW point 1.3 going if you're on OS 10 you can use King Kaos it sounds like a mid 80s hacker collective passing out floppy disk with Commodore 64 games on it but it's not it's these really nicely packaged binaries of gee doll and some of the supporting libraries and some of the things like QGIS that are built on GDL windows you can go through the UCLA has a really good demo and tutorials how to install it and then if you're running Linux you know way more about this than I do and you're on your own so once you have G deho installed and I'll give the URL in a second but I've written a blog that goes in detail in every single step and also has demo data sets that you can download but if you download a geo TIFF go ahead and run a command called GL info and the basic structure is you specify a command you point it at a data set and then you give it options and so you run GL info and you get back this big long string of text but the important things for for us is that a do TIF is just a normal kiss but it's got another header on it which has the information so that we can precisely place every single pixel on an image to a precise coordinate on the Earth's surface and so part of that is the coordinate system you'll also see corner coordinates so you can sort of know where the corners are the edges and place the entire image in context and then it'll give you information about each band individually and that becomes really important when we start working with the satellite data so one of the easiest functions you can do with G doll is G doll translate and that's just from going from hundreds of different data formats into whatever data format you want so for making a image of a map for this slideshow you can say geo translate is the command and then specify an output format it's simply - Oh F there's - Co a Co is a creation option and this is one of those kind of obscure really difficult parts of G tall in that they're specific to each individual file format so they're not really documented all in one place or spread it over the place but quality is the JPEG quality from 0 to 100 just set that to 90 outsize in this case is in pixels so in 1920 great for HD display and then that's 4x width you know I didn't specify height at all so Gino will actually keep the aspect ratio if you set one or the other variable to be zero and so that's sort of an important concept about GDL is that it tries to do a lot of the work for you which can be super convenient when it works right but occasionally will break spectacularly and like trying to figure out exactly what went wrong can be quite tricky then I specified a resampling method bi-linear is a nice way to do smooth it will default the nearest neighbor which gives stair steps and looks pretty terrible so it's a good thing to keep in mind and then just an input and output file and you end up with this map perfectly sized for the screen so so far I haven't really shown anything that's particularly unique to G tall you could do this in pillow you could do it an image magic but as I said before the earth isn't flat and so the power of G tall is really being able to work with things that are directly tied to real positions on a 3d surface but the earth looks flat from ground level it looks flat from an airliner or from the summit of Mount Everest and it's not until you get to about seventy seventy-five thousand feet that you can begin to see that curve and so this is from a balloon expedition launched from South Dakota in 1935 from just above seventy two thousand feet and it's documents the first time humans ever saw the curvature of the earth and it's not until you get to about a million miles away that you see the earth as this perfect disc against a black sky and from this far you're actually getting very close to that orthographic projection that Mike was showing when he's spinning his globe where it's basically seeing everything from ninety north to 90 South projected on to a perfect sphere so each of these views is in effect a different map projection because you're looking at the Earth from a certain perspective in space and then have a 3d object that then you're representing in 2d so the 2d of the screen it's something we've known about as the species and a culture for a long time this is the first globe it was published in 1492 before Columbus had returned from his sedition and so you can see that each of these sort of laws and shape things is what's called the gore and you would cut them out and then paste them on to a spherical sub right substrate to make the globe and so cartographers really have a challenge because there is no perfect representation of the 3d Earth in 2d and so a well-designed map is going to be solving a particular problem and for each different problem there's a specific map projection and G doll is this great tool for actually moving between the different projections so if we start with the natural earth which is this incredible free data set they have different scales of raster data in a bunch of different styles this one's really good for making base maps they also have some vector data but it's in what's known as a cylindrical equirectangular so it's just an even grid of latitude and longitude so it comes out as this two to one ratio and then in G tall we would invoke a function called GL warp to do the map projection and so some of the options are similar to what you would see in GTO translate and there is some consistency between commands unfortunately this syntax actually is quite different from command to command to command so you really have to pay attention and sometimes look it up and then at the bottom is a link to the blog series where again you can go through and I have examples and step through in more detail what each of these things is doing the next is specifying a target spatial reference system and so that's sort of the mathematical magic that allows each pixel to be precisely placed and rather than going into detail about every single component of this I will just point you to a spatial reference org where you would basically look up a map projection and it gives you an option of formats several options for formats and if you just click on the Prague's for link you'll get something that you then pasted in single quotes into that target source of spatial reference source system and do the map projection from there there are some map projections that get super complex mathematically where you can go in one direction but not the other and gdl really wants to be able to go in both ways so there are some map projections that just as a support and somewhat conveniently most of them are not supported by proj for so if you click on that link and get a blank then you probably can't do it in GL or at least not without a lot of extra work the other one of the other interesting flags that have set in this command is DSC alpha and what that does is all of the areas that are off the map because we're going from a rectangular map to a elliptical map get given an alpha channel so it'll default to just filling it with black but it's really nice to have an alpha so you can overlay it on anything and this next one is what's called the warp option and these are starting to really get into the nitty-gritty of GDL so it's manipulating how that very convoluted trigonometry is from going one map projection to the other is working in this case it's it's actually giving an extra border around the edges of the map so we can calculate the edges of the destination map which again isn't a perfect rectangle anymore and this is where again I was saying GDL does a lot of guesswork for you but it can go spectacularly wrong and so what the source extra is doing is just it's giving you a little extra boundary so if you run something and it doesn't look right or there's some missing information on the output try adding this line or looking up some of the other options that you have and then finally I just give an in-and-out file and I get this as the output and this is mole wide it's a eat what's called an equal area projection and so it's great for thematic maps so things like temperature rainfall and population density where you really want them like per square kilometer and if you were saying you know like showing change in temperature in Greenland it would give you a disproportionate view of how the world was changing but like most global projections the poles are highly distorted and shoved off to the edge of the map so what if you wanted a projection that focused right on one of the poles so there's another one called polar stereographic that's optimized for that but it blows up at the optimum opposite poles so you need to pre prepare your data set beforehand and so what I've done here is do a two step process instead of a one step process where I'm chopping out the autum 30 degrees of the globe so basically from the tip of the Antarctic Antarctic Peninsula south and then doing the next step of a jido warp and I'm doing that with what's called a V R T which is a virtual file instead of being a full-blown TIFF with all of your information in it it's a little tiny text file that is telling GDL how to do that transformation and so if you're doing a huge bulk job it's really useful could do things in terms of these V RTS instead of tips to save time and space and processing power and the next important bit was setting the projection window so negative 180 and the you know that's just defining the edges and so I get a little strip of the South Pole and then do the projection and the one interesting thing about how I did the projection versus the MOL wide is that referred to an EPS decode and so there's an Institute of petroleum engineers who have defined predefined a ton of different projections for like sort of well-known areas and so instead of that long string you can just enter in a numerical code and then you can look those up at epsg IO run the command and I get this which is lacking something and that's because I didn't include the source extra command and if you do that it makes sure to calculate the boundaries correctly as we mentioned earlier not only is the earth not flat it's not precisely round either there's a bulge in the middle and in the 19th century well before John Wesley Powell set out to map the canyons of the Colorado the East India Company did this incredibly detailed 70 or 80 year long expedition to survey India and they weren't setting out necessarily to make a highly detailed map of India and they weren't planning on determining the height of the highest mountains in the world to within 30 feet from 100 miles away although they did both of those things they were really trying to determine the shape of the earth the earth has an equatorial bulge and in fact it's pronounced enough that if we measured elevation in terms of distance from the center of the earth Chimborazo in Ecuador which is I think a quite modest 18,000 feet would be the telephone in the world instead of Mount Everest so these effects don't really matter much at a global scale or even a continental scale but they matter quite a bit if you're trying to connect to railroads in the center of a continent or drill a well precisely where you have our permit and this kind of shows up in some really complicated concepts in these projection strings that you'll see you don't really necessarily need to internalize all this but just be aware that they exist and be aware that if you are in a certain datum it's usually best to stay in that datum instead of trying to move it to another one because that may be an undefined problem or might require like a ton of coding to get it to work right so basically there's an ellipsoid or spheroid and this is a idealized mathematical shape there's a geoid which is an approximate shape which is equivalent to the mean sea level and then there's a datum which sort of ties these two things together and so some common datums you might see nad 83 or 27 and then there's individual ones for Great Britain in Europe but for a lot of satellite data it's wgs84 and if you can find data that's in this map projection or in that uses the statum which is confusing the also on map projection use it because it'll solve a lot of it will just prevent a lot of nasty problems so how would we go about making a global map you can either a local map can either start big and go small or start small and go big we'll start with the former this is using a continental scale map of the US a variant on the natural earth data set and then we're going to crop that down and reproject it and move from CONUS Albers projection to what's called the lambert conformal conic as i mentioned there's all these different map projections optimized for different things you really have four options area angle distance and direction and no map can do all for them perfectly and so there's various trade-offs and so to equip a larger larger map you can run this important part from what I've shown you before is again using a project for string to specify the spatial reference system but instead of pulling it from a library I used a really handy tool called projection wizard org go to that site you drag a box around the area you want you pick the parameter you want to optimize and then it'll spit out this string for you and then you can just paste it in to geo warp and then also you need to define the extent so these are the edges of the map I've done that and lat/long generally speaking Gao will default to the coordinate system of the destination that can be really tricky because it's often defined in meters versus some possibly arbitrary centroid so if you specify your own spatial reference system for those target extents in this case EPS epsg 43 26 which is that same lat/long you can just type in the latitude and longitude corners and then from that continental scale map you get this reasonably high resolution map of the canyon lens again but say we want into like a super high res map we can turn to another data set this example is from Naples photography it's done like every one or two or three years over the entire u.s. one meter per pixel or higher but it comes in these little tiny tiles so we're going to use a new program called geo merge this is actually was originally a demo program but it became widely used and is super just commonplace so it got kept in the codebase but unfortunately it's like all the syntax is inconsistent with all the other GDL syntax so again it's just something that you should really keep keep aware of so you can just run this and you specify an output file and then it's really nice because you can just give it a wild card so I'm just going start out TIFF and it takes every tip of my directory and stitches it into a nice local map so it's great just because of how simple it is relative to the rest of GDL and then what I've shown you so far is all pretty much static Maps you know a map is an idealized representation but the earth is always changing so we can turn to satellite data as this extremely powerful way of looking at earth as it evolved and you might have heard California had this extremely wet winter so this was what the Caruso plane which is where the San Andreas runs north of Los Angeles you can see that horizont that straight line running right through the middle the image this is what the area looks like in the winter and then this is what it looked like at the end of March in the spring so all this rain led to this incredible profusion of wildflowers you've probably seen some pictures but satellite data is another way of looking at how widespread this phenomenon really was and so G Dahl is once again a really good tool for working with that that specifically because satellite data is strictly RGB so it needs to be massaged to be useful and plugging my blog endlessly I have links to a sample data set on the blog but you can also sign up for my company's Open California data set where you can get all the data of California that's more than two weeks old under a Creative Commons license so if you download some of that RapidEye imagery and you open it up in a tissue or it would look like this so not only is it super dark the colors are also kind of wrong and that's because the satellite data is packaged in sort of a scientifically sensical way and not a otherwise sensical way so it's blue green red red edge near-infrared instead of red green blue so it's got not only visible channels but also these extra bands that are useful for doing specific measurements so in GTO we need to reorder or to view it for at all we need to reorder it and GDL provides a handy way of doing that we go back to gol translate and most of the command is the same but we specify each band individually with a flag it's just - B and so and then it goes in order so band 3 is red I'm putting that in the first channel Bantu is green and it stays where it is in the second Channel and then band one is blue and it's getting moved to the third channel and then one other thing I need to specify another of these really obscure creation options photometric and this sets the image to be interpreted as a full-color image rather than separate greyscale bands so if you were just to load this into Photoshop without setting this flag it would say like grey and then alpha 1 and then alpha 2 instead of giving you an RGB image we do that still dark scientific data is linear our sensory senses are nonlinear and so we want to do a little bit of scaling on the data to make it approximate how we see we can do that again with GL translate but I've run GL info first to get the image statistics that we need to then know the end point and then after that I can add this - scale flag and what's that's doing is it's taking the minimum and maximum extents and then stretching them into the full range of a 16-bit file so from 0 to 65,000 and then I also do an additional step of scaling exponentially again because we don't see linearly we have this curve that sort of compresses black and expands white lighter images lighter information and if we do that we get a much more viewable picture a satellite data color correction is like a whole nother thing I could go on for days and days about that but this is a really good first approximation so what would happen if I wanted to compare this year to last year to get a sense of you know how extraordinary of an event this is and I went back checked for rapid eye data and there was none that wasn't cloudy but there was some good Landsat data another data set this one's free provided by the US government you can get it off Amazon and I download the data for a day last spring and unlike RapidEye which is multiple bands in one file and Lance at every single band is its own file and so we need to use a different technique to do the combination again using GDR and using the separate flag to specify the red the green in the blue bands and you would probably have to look this up independently for every satellite because every satellite was made by a different team of engineers that wanted to precisely solve their own problem so they're all different it's totally nightmarish but for Landsat four three and two are red green and blue and you just combine them into a single file and then this is the color corrected version of this and you'll notice that it's a much wider area and you can't tell but it's also a lower resolution so how would we match this dataset with the rapid eye which comes in a little tiny tile we use a couple other commands so gelt index is actually walking around the corner of this the file and then creating a second vector file called the shape file which is a really common spatial data set but instead of storing things as pixels it's storing things as points with a latitude and longitude and you do that with the T index and create the shape and then you point G doll at the shape file to create a cut line so it's basically making a box and superimposing it on that larger image and then just tell it to do the cropping so the crop function and then the file you're cropping and then the output and the final little bit of GDL magic is I'm setting this target resolution with TR and this is not in pixels this is actually in real world units so rapid eye data is five meters each pixels five by five meters and so I just set that flag and then it resizes the Landsat data to match so this is March 25th of 2016 March 31st of 2017 and you can see 2015 was a fairly wet winter but not extraordinary 2017 was like biblical and so we got much greener and these like incredible colors all over the place so hopefully I've given you a fairly good introduction and maybe comfortable enough to start looking into this on your own some functions I'd encourage looking at our GTL diem Oh gr2 ogr there's some ways to wrap it our rasterize vector data you can pan sharpen low resolution data with high resolution data and then you can do banned math between like the red and near-infrared bands like addition multiplication subtraction to do some fairly sophisticated analysis purely indeed all references I recommend including Derrick Watkins geo cheat sheet which is a good thing to just keep open because it has specific commands and then syntax and then I'd like to thank all these people for proofreading and sitting down and showing the source code and all sorts of wonderful stuff and as a teaser this is the stamen watercolour map which was pulled out of web tiles and into a geo TIFF with a single command thanks [Applause]