Mapchete - parallelized batch geoprocessing using Python
CC Attribution 3.0 Germany:
You are free to use, adapt and copy, distribute and transmit the work or content in adapted or unchanged form for any legal purpose as long as the work is attributed to the author in the manner specified by the author or licensor.
Ungar, Joachim (EOX IT Services)
Processing geodata can be fairly simple until the input data reaches a certain size. Creating a hillshade or extracting contour lines from a DEM can be done quickly, but if you want to do this with e.g. the global SRTM dataset (1296001 x 417601 pixel), the process will crash (unless you are visiting from the future). Besides, if there are additional steps required like clipping the data to the 400MB landpolygon behemoths from OSM or applying custom filters, you probably find yourself starting to write your own tool chunking the data. mapchete tries to solve this issue by helping you to focus on developing your geoprocess written in Python and applying this process to the data. It does so by automatically reprojecting and chunking the input datasets into tiles (based on the “WMTS simple profile”) and running your Python process for each tile individually and in parallel on all available CPU cores. mapchete offers two command line tools. mapchete_execute runs the process on the full dataset, similar to tile pyramid seeding for map caches. mapchete_serve hosts an OpenLayers interface and processes only the data in areas and zoom levels you are currently inspecting. This allows you to test and assess your process on the full dataset on the server, instead of clipping and downloading subsets on your laptop. mapchete is used as the data preprocessing backbone of EOX Maps, a service which provides background maps for example to the European Space Agency.
 Hello? Okay, I've just checked in. Hey, welcome to the sessions this afternoon. I'm the chair for your next three sessions. I'm Johann. Sean, who was supposed to chair the sessions, sadly couldn't make it, but I'm here. The next talk is going to be a chat by Johann. Okay, thank you very much, Johann. Hello, everybody. My name is Johann. I work as a cartographer at EOX. We are a Viennese-based company. And we are doing mainly stuff with the European Space Agency. So, yeah. So, we are working mainly in the Earth observation field. But we're also doing maps. Great. Ah, yeah. Okay. So, we're also offering global background maps. So, for our own services and for some services for ESA, you can get them there with a WMS endpoint. And currently we're working on a new version. And, yeah, I want to talk about our tool we've been developing in the last year to create, to pre-process our map data. As you might know, if you ever start a map in the recent years, it became fairly easy now with tools like TimeMill, for example. So, there are a lot of great tools out there where you can just say, okay, here's my data. Here are my style files. Then these tools would provide you with maps. However, if you, for example, want to make a terrain layer and don't have a hill shading ready, then you somehow have to get this hill shading. And, of course, there's GDAL, or the command line tools from GDAL. They can do it. But it's probably not a good idea if you applied to a whole SRTM dataset. Because SRTM is really huge. It's 30 meters globally. And, yeah, it wouldn't work. The same would apply if you want to compute contour lines, for example. So, it wouldn't work because, first, it's not really viable because the input rusters are so big. And a second point, there is no customization available. So, if you want to tweak your hill shading afterwards by applying some image filters, for example, or doing some clipping by the coastlines, then it wouldn't work. Well, when generating maps, tile permits were introduced a couple of years ago. So, if you go to OpenStreetMaps, basically, the global map is cut into little tiles, and you just get served the area in the scale you requested, which is a pretty good approach. So, they dared with it by simply chunking the data into smaller parts, which are more easily digestible for processing. Yeah, and it works basically like this. Here, we have our data, then we define how our map looked like, and then the tool would provide the map. So, we thought, I mean, it would be great if we had, again, our input data and then some kind of a process, so not a styling, but any kind of geo process, and we would get our process output tile or whatever. And, yeah, we thought, yeah, this approach worked so well for maps, so there must be something out there that would do it with data, but, yeah, nothing really to be found. So, we started to create our own, and it's called MapShady. It's a tool written in Python. I mean, I could tell you a long list why Python is great and why we are using it. The main reason we started with it, because, or I started with it, is I'm a cartographer and not really used to programming, and Python is the only thing I'm comfortable with yet, so, yeah. Therefore, it was Python. But, on the other hand, Python has a really great ecosystem, so you have a lot of modules that will help you with any kind of geo processing tasks. So you have a SIP package, for example, or Numpy, both for raster or array processing, or Shapely, if you need vector processing. Well, MapShady, basically what it does, it chops your input data into tiles, and then applies any kind of user-defined process, so you write the process in Python, and it applies this process to every single tile. Under the hood, we are using Rasterio and Fiona to read raster and the vector data, which are great tools, both, by the way, and we support both the spherical and mercator protection, so the mostly used protection, and also the WGS84 let-long tile grids, which we use a little bit more because we rather have to show global data, and then this projection is slightly better because you can also cover the poles, for example. MapShady comes with two command line utilities. The first one is MapShady's surf. So what it does is it starts a web server, it's written in Flask, and then you can open your browser and go to port 500, and then you have an open layers instance, and as you zoom through your map and pan through your map, the tiles are being processed on demand. There's also an override flex, so if you turn this on, then you can edit your code, then hit the reload button, and you can immediately see what effects your editing had on your code. So it's very easy and great to immediately test what you have, because usually if you have big amounts of data, then you would either put it on the server and download regions of interest and process them locally on your laptop, or you do the processing on the server and then download the outputs afterwards, which is really messy. So with this MapShady surf, you're much faster in your development cycle. And there's a second tool. It's basically the seed command. It's called MapShady execute. It takes the same process, and it can process the full tile pyramid, or you can limit it to certain zoom levels, so you can say, I just want zoom level 0 to 8, or simply process a specific tile if you know the tile index, or of course you can do a spatial subset as well, so you can tell just render me this bounding box. So this is the configuration that it should look like. It's a YAML file, and you need some certain input values, of course the path to your Python file, then the input files you're using, so any kind of raster or vector data, and some information on the output format. So the path where it has to be stored in, and the format we're currently supporting PNG, then a special PNG version, so we developed it for HealShading. What it does here is that it writes the HealShade information into the alpha channel and makes the other channels black, because then you can overlay the HealShade over any kind of other data or other map layers. We also support of course GOTIF, so we use this to store our global elevation models, or use GOTAS and our post-GS for vector data, and we also experimented in dumping the raw arrays as numpy blobs on the file system. Here you specify your output protection, so it's either a Geodetic or MacCator, and the data type and bands and no data you will need for raster data saving. Then there are also optional commands, so you can specify the zoom level validity, so to say, for example, your map-shake process can be only available between zoom level 0 and 8, for example, or you can do a spatial boundary here as well. We also support metatiling, which means that you can say it should not... Metatiling basically just increases the tile size. So on your web client, it will still be the small tiles, but in the back end, it will render these tiles by 8 by 80 tiles, so 64 times, which in some cases really makes the processing faster. And of course, you can specify down here your process-specific parameters. We also developed a basic configuration parser which would allow you to specify zoom-dependent, zoom-level-dependent parameters. For example, for the heel shading, if you want a higher exaggeration of the elevation model before calculating the heel shade, and zoom-levels below 5, for example, then you can specify it here. So it always takes... So when processing a tile, it always takes the configuration view of the zoom level, which comes really handy if you're dealing with differently simplified vector data, for example. Yeah, this is how your process file should look like. We tried to use the same paradigm as WPS does. So basically, you have your process class up here, and then your execute function, and you can write basically anything here and you can import any kind of Python module you need for your processing. To open data, there's an open function, which can be used both for vector and for raster data. It also has a function which will return true if the tile of your source data is empty or not, which comes handy if your process is depending on this data, and it immediately quits, makes it faster. And yeah, you can, for example, to read RGB data, you simply write this. It's very similar to Rasterio and how they implemented it. That should not have happened. And yeah, there's also a global write function, which will either get a tuple of your bands if you're using raster data or a single array if you have just one band, or you can provide it with a list of features and then it would simply write it as vector data, so depending on how you define your output. So we really tried to make it very easy so that you don't have to care about reading or writing and simply focus on what you want to do with your data. And there are also some common functions we introduced. So there's a hillshade function, which comes handy, then a contour extraction function and a clip. So the clip function you can provide it with raster data and a vector data set, and it will clip the data. Also, you can specify some buffers here, so you can say, I don't want the exact coastline, but give me a buffer around the coastline and so it would add this buffer to a clip. This list could be extended, but for our purposes, at least these three functions were really used a lot. Yeah, what did we use it for? As I told you, for our maps. So we used the list of elevation models collected by the open terrain community and created our own global elevation models based on ASTAR and SRTM and some regional data sets. We used MapShady to merge this elevation data and to clip it to the coastline boundaries and also combine it with some bathymetry data. You see it here. It's from CleanTopper. It's a sea surface elevation. We also created our EOX-specific hillshading. So we worked on this rather a lot. We tried to get as close as possible to the great Swiss Topo hillshading. Yeah, it's not that easy to get it, but we are quite satisfied with it. But yeah, to get this hillshading, the normal hillshade function wouldn't work because we use hillshades from different resolutions and combine them and then applied some median filters to get the best out of hillshades. We also used it to show you the vector functionalities. I don't know whether it can be seen, but we are actually having multicolor contour lines on our maps. So we make the contour lines blue where glaciers are, so like the great topographic maps have. And the contour lines are brown on land. And we also used the tool that does the clipping and the intersecting with the glacier mask. And we also dug into Sentinel-2 data. So for this example, we collected or we downloaded data for, I think I guess it was a time span from two or three months and sorted the pixels by brightness and then picked out the pixel that's most likely to be not cloudy. Yeah, and for this version, we had to experiment with the numpy dumping because when you dump the numpy array, you are not bound to two or three dimensions, but you can basically dump a time series of arrays as well. Yeah, this is another area in Austria. It's very cloudy, especially during the time, but we are, I mean, for our first version, we are quite satisfied with the result. So this is basically a pixel by pixel mossoic within a time frame of three months. Yeah, Mapchetti is available on PyPy, so you can use PIP to install it. And yeah, we have GitHub as well, and it also has a Twitter account. And yes, well done. Thank you very much. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you. Thank you. Cool, let's go into questions. Are there questions? I'm going to run to the back first. Cool, thanks. That was a great talk. I had a question around, I think I read in the blurb, that the Mapchetti can run in parallel, or it's parallel. Yes, if you. And kind of, is that much faster than running it serial, or is it, because I know that Python is sometimes not particularly great at parallel processing, because of the global interpreter lot? Yes, it's definitely faster. I mean, it depends on the kind of process you have. If the bottleneck is this guy, for example, and it doesn't make much sense. But for most of our use cases, it worked well, although you have to tweak a lot with the metatile size and with the numbers of processes. Right, okay. Cool, thanks. You were speaking about the projection, and right now it's supported only StereoCAD and VGS.84. Is it possible also to add the other projection system, or not? You mean for as output format? Yeah. Okay. Well, it should be theoretically possible. But for our use cases, we just need these two projections. So we, yeah, maybe we do it one day. What I forgot to tell you, what we did as well, for the global projections is if you read the tile with a buffer near the anti meridian, then it will take the data from the other side as well and stick it together. And this would make it a little bit difficult to define your regional tile metrics. I mean, it would be possible, but it's pretty hard coded right now. Well, thank you. Very interesting address, very interesting tool. I'm just wondering why didn't you use directly PyWPS and code a lot of stuff that is already there? Yeah, I mean, I found out about it a little bit later. I didn't think about it, that PyWPS could do this. And then I was in the middle and thought, ah. It does everything. Really? Another question, which technology do you use for parallelizing? Multiprocessing? Yeah, the multiprocessing module from Python. So I will just leave you the challenge to integrate this with PyWPS4 if you can. Okay. A bit more programming question. The thing is that you have multiple tiles and you do processing parallel, and then you go looping through the tiles or apply a map, or... And also another question is that, have you tried running this on other weird Python interpreters like PyPy? No. Okay, no. Yeah, and of the tiles, in the beginning I calculate the tiles that will intersect with my input data, and I don't run over tiles that are not needed. So I do a little bit of a preselection before. And then I simply loop through these tiles. Yeah. Probably if you start doing loops, probably if you use PyPy interpreter, you're going to get to know this faster. Ah, okay. Because it has just an interpreter. Okay, I didn't tuck in into that deep, but I'll have a look at it. Thanks. I have time for more questions, by the way. So three, four more questions or remarks are okay. Thank you. How did you get cloudless images? Did you just... select by hand or... No. I took all the images I could get, and then there's a great blog post from Mapbox about it, so they did it the same with Lanza data. So you basically build up pixel stacks, and then you do pixel by pixel sorting by brightness, because the brightest pixel would most probably be cloudy, and the darkest one would most probably have cloud shadows in it. And then, yeah, you have to tweak around. It's some pixel in the middle that it's almost fine. Okay, thank you. If I understand correctly, it's mostly for mapping that you're doing this. Pre-processing for mapping. Yeah, but if you go a bit further in terms of analysis, how do you handle borders? Do you have... Because you can have types of processes where you can have border effects if you don't work with overlaps and... Ah, I forgot to mention that you can provide a buffer value around the tiles. It's required for hill shading, for example, because the hill shading process, it always leaves you no data border, so it handles this. And if you're at the anti-marillion, it would read the data from the other side and get it. So this is... Okay, anyone want to ask a question? Okay, then we're finishing the session. Thank you.
 Hello?  Okay, I've just checked in.  Hey, welcome to the sessions this afternoon.  I'm the chair for your next three sessions.  I'm Johann.  Sean, who was supposed to chair the sessions, sadly couldn't make it, but I'm here.  The next talk is going to be a chat by Johann.  Okay, thank you very much, Johann.  Hello, everybody.  My name is Johann.  I work as a cartographer at EOX.  We are a Viennese-based company.  And we are doing mainly stuff with the European Space Agency.  So, yeah.  So, we are working mainly in the Earth observation field.  But we're also doing maps.  Great.  Ah, yeah.  Okay.  So, we're also offering global background maps.  So, for our own services and for some services for ESA, you can get them there with a WMS endpoint.  And currently we're working on a new version.  And, yeah, I want to talk about our tool we've been developing in the last year to create,  to pre-process our map data.  As you might know, if you ever start a map in the recent years, it became fairly easy now with tools like TimeMill, for example.  So, there are a lot of great tools out there where you can just say, okay, here's my data.  Here are my style files.  Then these tools would provide you with maps.  However, if you, for example, want to make a terrain layer and don't have a hill shading ready, then you somehow have to get this hill shading.  And, of course, there's GDAL, or the command line tools from GDAL.  They can do it.  But it's probably not a good idea if you applied to a whole SRTM dataset.  Because SRTM is really huge.  It's 30 meters globally.  And, yeah, it wouldn't work.  The same would apply if you want to compute contour lines, for example.  So, it wouldn't work because, first, it's not really viable because the input rusters are so big.  And a second point, there is no customization available.  So, if you want to tweak your hill shading afterwards by applying some image filters, for example,  or doing some clipping by the coastlines, then it wouldn't work.  Well, when generating maps, tile permits were introduced a couple of years ago.  So, if you go to OpenStreetMaps, basically, the global map is cut into little tiles,  and you just get served the area in the scale you requested,  which is a pretty good approach.  So, they dared with it by simply chunking the data into smaller parts,  which are more easily digestible for processing.  Yeah, and it works basically like this.  Here, we have our data, then we define how our map looked like,  and then the tool would provide the map.  So, we thought, I mean, it would be great if we had, again, our input data  and then some kind of a process, so not a styling, but any kind of geo process,  and we would get our process output tile or whatever.  And, yeah, we thought, yeah, this approach worked so well for maps,  so there must be something out there that would do it with data,  but, yeah, nothing really to be found.  So, we started to create our own, and it's called MapShady.  It's a tool written in Python.  I mean, I could tell you a long list why Python is great  and why we are using it.  The main reason we started with it, because, or I started with it,  is I'm a cartographer and not really used to programming,  and Python is the only thing I'm comfortable with yet, so, yeah.  Therefore, it was Python.  But, on the other hand, Python has a really great ecosystem,  so you have a lot of modules that will help you with any kind of geo processing tasks.  So you have a SIP package, for example, or Numpy, both for raster or array processing,  or Shapely, if you need vector processing.  Well, MapShady, basically what it does, it chops your input data into tiles,  and then applies any kind of user-defined process,  so you write the process in Python,  and it applies this process to every single tile.  Under the hood, we are using Rasterio and Fiona to read raster and the vector data,  which are great tools, both, by the way,  and we support both the spherical and mercator protection,  so the mostly used protection, and also the WGS84 let-long tile grids,  which we use a little bit more because we rather have to show global data,  and then this projection is slightly better because you can also cover the poles, for example.  MapShady comes with two command line utilities.  The first one is MapShady's surf.  So what it does is it starts a web server, it's written in Flask,  and then you can open your browser and go to port 500,  and then you have an open layers instance,  and as you zoom through your map and pan through your map, the tiles are being processed on demand.  There's also an override flex, so if you turn this on,  then you can edit your code, then hit the reload button,  and you can immediately see what effects your editing had on your code.  So it's very easy and great to immediately test what you have,  because usually if you have big amounts of data,  then you would either put it on the server and download regions of interest  and process them locally on your laptop,  or you do the processing on the server and then download the outputs afterwards,  which is really messy. So with this MapShady surf,  you're much faster in your development cycle.  And there's a second tool. It's basically the seed command.  It's called MapShady execute. It takes the same process,  and it can process the full tile pyramid,  or you can limit it to certain zoom levels,  so you can say, I just want zoom level 0 to 8,  or simply process a specific tile if you know the tile index,  or of course you can do a spatial subset as well,  so you can tell just render me this bounding box.  So this is the configuration that it should look like.  It's a YAML file, and you need some certain input values,  of course the path to your Python file,  then the input files you're using, so any kind of raster or vector data,  and some information on the output format.  So the path where it has to be stored in,  and the format we're currently supporting PNG,  then a special PNG version, so we developed it for HealShading.  What it does here is that it writes the HealShade information  into the alpha channel and makes the other channels black,  because then you can overlay the HealShade  over any kind of other data or other map layers.  We also support of course GOTIF,  so we use this to store our global elevation models,  or use GOTAS and our post-GS for vector data,  and we also experimented in dumping the raw arrays  as numpy blobs on the file system.  Here you specify your output protection,  so it's either a Geodetic or MacCator,  and the data type and bands and no data  you will need for raster data saving.  Then there are also optional commands,  so you can specify the zoom level validity,  so to say, for example, your map-shake process  can be only available between zoom level 0 and 8, for example,  or you can do a spatial boundary here as well.  We also support metatiling,  which means that you can say it should not...  Metatiling basically just increases the tile size.  So on your web client, it will still be the small tiles,  but in the back end, it will render these tiles  by 8 by 80 tiles, so 64 times,  which in some cases really makes the processing faster.  And of course, you can specify down here  your process-specific parameters.  We also developed a basic configuration parser  which would allow you to specify zoom-dependent,  zoom-level-dependent parameters.  For example, for the heel shading, if you want  a higher exaggeration of the elevation model  before calculating the heel shade,  and zoom-levels below 5, for example,  then you can specify it here.  So it always takes...  So when processing a tile,  it always takes the configuration view of the zoom level,  which comes really handy if you're dealing with  differently simplified vector data, for example.  Yeah, this is how your process file should look like.  We tried to use the same paradigm as WPS does.  So basically, you have your process class up here,  and then your execute function,  and you can write basically anything here  and you can import any kind of Python module  you need for your processing.  To open data, there's an open function,  which can be used both for vector and for raster data.  It also has a function which will return true  if the tile of your source data is empty or not,  which comes handy if your process is depending on this data,  and it immediately quits, makes it faster.  And yeah, you can, for example,  to read RGB data, you simply write this.  It's very similar to Rasterio and how they implemented it.  That should not have happened.  And yeah, there's also a global write function,  which will either get a tuple of your bands  if you're using raster data or a single array  if you have just one band,  or you can provide it with a list of features  and then it would simply write it as vector data,  so depending on how you define your output.  So we really tried to make it very easy  so that you don't have to care about reading or writing  and simply focus on what you want to do with your data.  And there are also some common functions we introduced.  So there's a hillshade function,  which comes handy,  then a contour extraction function and a clip.  So the clip function you can provide it with raster data  and a vector data set, and it will clip the data.  Also, you can specify some buffers here,  so you can say, I don't want the exact coastline,  but give me a buffer around the coastline  and so it would add this buffer to a clip.  This list could be extended, but for our purposes,  at least these three functions were really used a lot.  Yeah, what did we use it for?  As I told you, for our maps.  So we used the list of elevation models  collected by the open terrain community  and created our own global elevation models  based on ASTAR and SRTM and some regional data sets.  We used MapShady to merge this elevation data  and to clip it to the coastline boundaries  and also combine it with some bathymetry data.  You see it here.  It's from CleanTopper. It's a sea surface elevation.  We also created our EOX-specific hillshading.  So we worked on this rather a lot.  We tried to get as close as possible  to the great Swiss Topo hillshading.  Yeah, it's not that easy to get it,  but we are quite satisfied with it.  But yeah, to get this hillshading,  the normal hillshade function wouldn't work  because we use hillshades from different resolutions  and combine them and then applied some median filters  to get the best out of hillshades.  We also used it to show you the vector functionalities.  I don't know whether it can be seen,  but we are actually having multicolor contour lines on our maps.  So we make the contour lines blue where glaciers are,  so like the great topographic maps have.  And the contour lines are brown on land.  And we also used the tool that does the clipping  and the intersecting with the glacier mask.  And we also dug into Sentinel-2 data.  So for this example, we collected or we downloaded data  for, I think I guess it was a time span from two or three months  and sorted the pixels by brightness  and then picked out the pixel that's most likely to be not cloudy.  Yeah, and for this version, we had to experiment  with the numpy dumping because when you dump the numpy array,  you are not bound to two or three dimensions,  but you can basically dump a time series of arrays as well.  Yeah, this is another area in Austria.  It's very cloudy, especially during the time,  but we are, I mean, for our first version,  we are quite satisfied with the result.  So this is basically a pixel by pixel mossoic  within a time frame of three months.  Yeah, Mapchetti is available on PyPy,  so you can use PIP to install it.  And yeah, we have GitHub as well,  and it also has a Twitter account.  And yes, well done.  Thank you very much.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Thank you.  Cool, let's go into questions.  Are there questions?  I'm going to run to the back first.  Cool, thanks.  That was a great talk.  I had a question around, I think I read in the blurb,  that the Mapchetti can run in parallel, or it's parallel.  Yes, if you.  And kind of, is that much faster than running it serial,  or is it, because I know that Python is sometimes  not particularly great at parallel processing,  because of the global interpreter lot?  Yes, it's definitely faster.  I mean, it depends on the kind of process you have.  If the bottleneck is this guy, for example,  and it doesn't make much sense.  But for most of our use cases, it worked well,  although you have to tweak a lot with the metatile size  and with the numbers of processes.  Right, okay.  Cool, thanks.  You were speaking about the projection,  and right now it's supported only  StereoCAD and VGS.84.  Is it possible also to add the other projection system,  or not?  You mean for as output format?  Yeah.  Okay.  Well, it should be theoretically possible.  But for our use cases, we just need these two projections.  So we, yeah, maybe we do it one day.  What I forgot to tell you, what we did as well,  for the global projections is if you read the tile  with a buffer near the anti meridian,  then it will take the data from the other side as well  and stick it together.  And this would make it a little bit difficult  to define your regional tile metrics.  I mean, it would be possible, but it's pretty hard coded right now.  Well, thank you.  Very interesting address, very interesting tool.  I'm just wondering why didn't you use directly  PyWPS and code a lot of stuff that is already there?  Yeah, I mean, I found out about it a little bit later.  I didn't think about it, that PyWPS could do this.  And then I was in the middle and thought, ah.  It does everything.  Really?  Another question, which technology do you use for parallelizing?  Multiprocessing?  Yeah, the multiprocessing module from Python.  So I will just leave you the challenge to integrate this  with PyWPS4 if you can.  Okay.  A bit more programming question.  The thing is that you have multiple tiles  and you do processing parallel,  and then you go looping through the tiles or apply a map,  or...  And also another question is that,  have you tried running this on other weird Python interpreters  like PyPy?  No.  Okay, no.  Yeah, and of the tiles,  in the beginning I calculate the tiles  that will intersect with my input data,  and I don't run over tiles that are not needed.  So I do a little bit of a preselection before.  And then I simply loop through these tiles.  Yeah.  Probably if you start doing loops,  probably if you use PyPy interpreter,  you're going to get to know this faster.  Ah, okay.  Because it has just an interpreter.  Okay, I didn't tuck in into that deep,  but I'll have a look at it.  Thanks.  I have time for more questions, by the way.  So three, four more questions or remarks are okay.  Thank you.  How did you get cloudless images?  Did you just...  select by hand or...  No.  I took all the images I could get,  and then there's a great blog post from Mapbox about it,  so they did it the same with Lanza data.  So you basically build up pixel stacks,  and then you do pixel by pixel sorting by brightness,  because the brightest pixel would most probably be cloudy,  and the darkest one would most probably have cloud shadows in it.  And then, yeah, you have to tweak around.  It's some pixel in the middle that it's almost fine.  Okay, thank you.  If I understand correctly, it's mostly for mapping that you're doing this.  Pre-processing for mapping.  Yeah, but if you go a bit further in terms of analysis,  how do you handle borders?  Do you have...  Because you can have types of processes where you can have border effects  if you don't work with overlaps and...  Ah, I forgot to mention that you can provide a buffer value around the tiles.  It's required for hill shading, for example,  because the hill shading process, it always leaves you no data border,  so it handles this.  And if you're at the anti-marillion,  it would read the data from the other side and get it.  So this is...  Okay, anyone want to ask a question?  Okay, then we're finishing the session.  Thank you.
foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_0_407.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_408_647.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_648_719.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_720_887.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_888_6179.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_6180_8939.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_8940_21383.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_21384_22067.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_22068_24455.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_24456_24539.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_24540_24743.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_24744_24887.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_24888_25931.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_25932_26639.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_26640_27443.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_27444_28151.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_28152_28631.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_28632_29207.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_29208_29435.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_29436_29555.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_29556_30251.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_30252_31511.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_31512_32051.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_32052_32399.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_32400_32747.png foss4g-2016-1252-mapchete_-_parallelized_batch_geoprocessing_using_python-hd_32748_32937.png