Adding Orthoimagery Support

[KermMartian]

I'm looking to patch in support for fetching orthoimagery along with the landcover and elevation data, and use the data from landcover + orthoimagery combined to pick the blocks to be used for the surface terrain on the final MC map. Do you have any suggestions or comments on whether the USGS services that region.py pulls from offer such data? From looking through their websites, it looks like a qualified "yes", but I'm having a great deal of difficulty finding documentation on the API(s) they offer. Thanks in advance, and I will be more than happy to send you a pull request if and when I get what I'm trying to do working if it's something you'd be interested in having in the TopoMC trunk.

[mathuin]

This is an interesting idea!

The best way to see what the USGS has to offer is to use http://viewer.nationalmap.gov/example/eros_inventory_proxy_json_v2.html. For example, the Provincetown region has a product with the key "P19" and the title "NAIP (4 Band) UTM Zone 19N". Select the product, and you can choose one of several slices of the Provincetown area to download. Once you've downloaded the file, you will have to mess about with it to figure out exactly how it is formatted. The metadata (available at http://extract.cr.usgs.gov/distmeta/servlet/gov.usgs.edc.MetaBuilder?TYPE=HTML&DATASET=NAIP_Z19_G4 for that type of file) will give you some information, but you'll need to use some GIS software or the GDAL binaries to really understand the file: things like what projection and scale are being used, what the pixel values mean for each band, all that stuff. Once you've done that, you can reproject the file to match the elevation and landcover files and figure out what block you want each pixel to represent.

After you've got that figured out, and the map turns out the way you want it to, you'll need to sample other parts of the country to determine all the product keys that return your imagery of choice (I imagine there's one per UTM zone, whatever that means), and start hacking region.py. Start by adding a new key-value pair to productIDs in region to support imagery -- the list of product keys you just generated goes here. If there are different types of imagery based on scale, you'll want to write some code to handle that -- use the elevation stuff as a template. Add another line to the end of the init routine like the ones for lclayer and ellayer, and go through the code adding your third layer to all the stuff that calls those two. Translate the pixels if necessary to a desired format, and add another layer to the existing Map.tif (there are currently four -- elevation, landcover, crust, and bathy). Then you'll have to modify terrain.py and friends to place the proper block on the surface to match the pixel of imagery.

It'll be a lot of work, but if you do it right, it'll make a big difference. Let me know if you need more help!

Jack.

[KermMartian]

Thanks for the very thorough feedback! I have gotten started trying to understand the data that EROS provides. I understand now that the TIFFs provide (R, G, B, IR), which should be interesting for things like figuring out where to place trees, but one thing I'm stuck on is figuring out how you got a product key of P19. I see various Zone IDs ranging from Zone 12 to Zone 19; are those the product keys? Or is there somewhere else I should be looking.

[mathuin]

The product key came from hacking region.py actually -- I knew I forgot something! Add 'SEAMTITLE' to the desiredAttributes value in region.py and then run ./GetRegion.py on the region of your choice with --debug enabled. You will get a ton of XML back and forth between you and the USGS, and one of those big XML blobs will contain the title and product key for every product offered by the USGS for that region. P19 is the product key for the product mentioned above.

Sorry about that! For what it's worth, the reason I had such a detailed response is that I'm looking into adding the 'terrestrial ecosystems' layer to properly set the biomes, and the process is the same. I'm not yet sure if that's better or worse than the landcover values, but I'm experimenting. :-)

Jack.

[mathuin]

http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system

Fun fact: there are ten UTM zones across the continental US. No idea about Puerto Rico, Hawaii, or Alaska. The GDAL transformations should be pretty straightforward from the UTM zone projection to the one used by TopoMC.

[KermMartian]

Thanks for your continued help. I got as far as puling the orthoimagery, joining and warping it, and applying that data onto the terrain. The color/materials that appear on the terrain match the orthoimagery TIF, which is great, but unfortunately that imagery is missing several large areas (and also seems to be scaled incorrectly). For the scaling, I may have made coding errors; for the missing areas, is it possible that some data is available in some sub-regions and not in others, and may need to be overlapped with other products? I notice that if I decrease the tile size and the scaling, the process fails entirely, as it seems none of my imagery (P10-P18, OF9, and ONF) is available over the entire shrunken region. Thanks in advance! I'd be happy to provide code if it helps.

[mathuin]

Can you attach some images? Specifically, the TIF and the test terrain would help.

If you hack region.py and add SEAMTITLE, any product you find that way will cover the entire region -- not always in one file of course -- but it may need tweaking with projections and such.

Also, have you looked at the TIF and the elevation/landcover files with something like qgis?

[KermMartian]

I have looked at all the files with MapWindow GIS on Windows, and I'm starting to understand how all the data (and metadata) fits together. I found the reason that the data appeared to be missing (P10-P19 have STATUS=Seamless, not STATUS=Tiled), so I must modify my question slightly: How can I use Seamless-type products, if you know? With the current code, the file list request for P18 (the elided coordinates are entirely within P18):

<REQUEST_SERVICE_INPUT><AOI_GEOMETRY><EXTENT><TOP>xxxx</TOP><BOTTOM>yyyy</BOTTOM><LEFT>zzzz</LEFT><RIGHT>aaaa</RIGHT></EXTENT><SPATIALREFERENCE_WKID/></AOI_GEOMETRY><LAYER_INFORMATION><LAYER_IDS>P18</LAYER_IDS></LAYER_INFORMATION><CHUNK_SIZE>250</CHUNK_SIZE><JSON></JSON></REQUEST_SERVICE_INPUT>

Returns an error:

<?xml version="1.0" encoding="UTF-8" standalone="no"?><REQUEST_SERVICE_RESPONSE><ORIGINATOR>RVS</ORIGINATOR><ERROR><FIELD_NAME>Ipc</FIELD_NAME><ERROR_MSG>Invalid product code: </ERROR_MSG></ERROR></REQUEST_SERVICE_RESPONSE>

TIA!

[mathuin]

Seamless content is a little different than tiled content. Before last December, TopoMC used seamless content, so you can look at the history of region.py to see what I was doing there. I had to change because elevation and landcover were no longer available as seamless content. That will make your work a little more frustrating, I fear.

How exactly did you retrieve the TIF? Any chance you can attach it to the issue, or give me a URL so I can get it?

[KermMartian]

The TIF I had had was from layer OF9, which is shady. P10-P19 are much better, and thanks to pulling in a great deal of download code from your pre-March region.py, I was able to get seamless data to be downloaded and cached. The code is ugly and needs to be cleaned up, but it works. Here's a 50%-scaled terrain using just R/G/B data without bothering with landcover types and taking advantage of the IR channel:

[mathuin]

That is really, really neat! I wonder how it'll look if you add in elevation?

Jack.

[KermMartian]

Actually, that's with all the layers including elevation already there. Here's another pair of renderings (from DynMap this time): http://www.cemetech.net/forum/viewtopic.php?p=206390#206390 . The water in the reservoir got replaced accidentally, but that's a very easy fix. My next challenge is that when I try to fetch the orthoimagery for the entire island of Manhattan, the USGS API tells me that my AOI is too large by roughly a factor of two. Do you know of a cleaner/easier way to resolve this than splitting the rectangular region in two (or N) pieces and fetching two (or N) TIFs using the USGS service and then passing them all to gdalwarp? TIA.

[mathuin]

No, that's pretty much how I'd do it. There's code for the gdalwarp merging in region.py -- it might be in the seamless stuff but it's definitely in the tiled stuff.

Have you considered only using the orthographic data when it corresponds to built-up areas (21-24) and letting landcover do the rest? That would let you get water and grass and beach while focusing the "interesting" coloration on the hardest-to-simulate-by-hand terrain...

Oh, and the dynmap stuff is awesome. Very nice work!

[KermMartian]

What I currently do is apply a different "palette" of blocks to landcover types 22 through 25 than to all the others, so that things like roads and paths in wooded areas (like parks) will still show up properly. I still have a bit of tweaking to do, and I might tweak it further to only include grass, sand, gravel, and a few other entities, in which case I would use the OpenStreetMap API to try to pull in and render streets and paths onto wooded/rural areas. For trees and buildings, I'm hoping to get permission from the Google Earth team to use their building and tree layers, but that's going to be difficult. I'm currently testing my building converter on single models, as you might have seen from that topic, for example with the Chrysler Building. My newest issue, from my work placing models onto the generated terrain today, is that the terrain seems to be rotated from "true" north by slightly under 15 degrees clockwise. Any insight into what might be causing that?

[mathuin]

http://www.charlespetzold.com/etc/avenuesofmanhattan/index.html

I would have expected you to see a rotation of about thirty degrees, based on that post. I wonder if there's some error due to the projection -- it's centered on Kansas, I think. Interesting question.

[KermMartian]

Based on reading up on the North American Datum, it seems to me that all points within the region I'm converting would be offset by roughly an equal amount of the issue was (for example) a NAD27/NAD83 mismatch. I see a rotation instead, so I suspect instead that something within my process of creating or warping the original input files into the VRT in the GetRegion step or warping the data in the PrepRegion step is causing the error. I note that to the extent the landcover array is correct, the bathy, orthoimagery, and elevation information all match up in the reservoir in the middle of the park in my sample data, so whatever error is occuring is applying to all of the layers. I note the rotations indicated below.