Mapping Vietnam
This is about building a web mapserver to share maps of Vietnam
Obtain source maps
If you already know the numbers of the maps you need,
- Go to USGS site http://store.usgs.gov/
- Click Advanced Search
- Enter map number in Old Material
- Click "Search"
- Click on the map image
- Save the PDF file
Samples for testing
- 1501ANC4810
- 1501ANC4811
- 1501ANC4815
You can also use the search features to find what you need.
Working with GeoPDFs
The GeoPDFs downloaded from USGS were scanned from paper maps, converted to a digital form (probably TIFF) and then georeferenced. Then they were converted to GeoPDF format using TerraGo Publisher.
Each scan covers only a relatively small area on the ground. So to create a complete seamless map for use on the web, we must
- Convert GeoPDFs to some format that we can further process.
- Remove the collars (the border containing legend, dates, etc).
- Process the set of resultant files into tiles that can be served.
Step one: Reading GeoPDF
If you are using ArcMap you can use "TerraGo Publisher for ArcGIS" ($2700) to import GeoPDFs directly. At this point, you can treat the GeoPDF data as you would any other ArcMap data source. Unfortunately $2700 is outside my budget.
Another excellent commercial product is Global Mapper (retails for $349). It can read the PDFs and export them as GeoTIFFs. An added plus, it has code to remove the collars and can also reproject the files if needed. I am told it is scriptable too.
You could ignore the geospatial encoding in the source PDFs and convert them into a format (TIFF for example) readable by ArcMap and then georeference it, but that's a lot of work given the PDF is already georeferenced!
The very latest version of GDAL (1.8.0) can read GeoPDFs. It's free and I have worked with it in the past, so I am giving it a try.
Recipe to convert a GeoPDF to a GeoTIFF using GDAL
% gdal_translate -of GTiff -co "TILED=YES" -co "TFW=YES" 1501ANC4815_geo.pdf 1501ANC4815_geo.tif
This generated a 47 MB TIFF that is 4470 x 3375 pixels.
You can learn more about GeoTIFF options here: http://www.gdal.org/frmt_gtiff.html A good one to use might be "-co "COMPRESS=LZW"
Step two: Removing collar
With Global Mapper this is as simple as checking an option box for the layer.
Working with my sample Vietnam maps, what I want to do is to detect the border of the actual map, then trim away the outside.
The actual trimming is easy, here is an expanded GDAL recipe, using pixels (trial and error): -srcwin xoff yoff xsize ysize
% gdal_translate -srcwin 350 0 4120 2750 -of GTiff -co "TILED=YES" 1501ANC4815_geo.pdf 1501ANC4815_geo.tif
According to the "Acrobat Supplement to the ISO 32000" there is a definition in the file to tell me where the map lives. The gdalinfo command reveals the NEATLINE setting:
NEATLINE=POLYGON (( 445035.80158801417565 999277.556482305750251, ul 444993.235997625160962 991902.349403560161591, 444864.738057852664497 884201.712978136376478, ll 452159.580137342331 884208.817744017927907, 585964.54743435990531 884289.527591207530349, 614300.221590773784555 884383.269830735283904, lr 614190.235319046885706 999293.477391553344205, ur 605175.140901734353974 999218.552309990860522, 553978.789515793323517 999300.979716640315019, 505958.457250406034291 999353.427943945163861, 451172.635888226970565 999382.750743734533899, 445035.722449889173731 999360.237513555795886, 445035.80158801417565 999277.556482305750251, 445035.80158801417565 999277.556482305750251, 445035.80158801417565 999277.556482305750251))
The corners are given as:
Corner Coordinates: Upper Left ( 429999.053, 999725.915) (104d21'46.72"E, 9d 2'39.27"N) Lower Left ( 430135.804, 856853.093) (104d21'58.76"E, 7d45' 6.78"N) Upper Right ( 619287.590, 999910.271) (106d 5' 7.80"E, 9d 2'41.48"N) Lower Right ( 619424.341, 857037.449) (106d 4'59.36"E, 7d45' 9.52"N) Center ( 524711.697, 928381.682) (105d13'28.20"E, 8d23'57.68"N)
Size is 4470, 3375
It seems to make sense that I assume that the corners correspond to the edges of the "paper", the latlon numbers look about right. It also seems I should be able to plug the neatline numbers into a gdal_translate command, but why are there so many, why not just four pair?
I put the neatline points into a text file, then I loaded the TIFF into QGIS, and used the "Add delimited text file" plugin to add the text file to the map. The points all fell neatly onto the neatline. So I should be able to select the limits (the corners) and pass them to gdal_translate.
Being basically a lazy human being I wrote a little Python to sort them.
Using map coords (using data from map): -projwin ulx uly lrx lry
% gdal_translate -projwin 104.5 9 106 8 -of GTiff -co "TILED=YES" 1501ANC4815_geo.pdf 1501ANC4815_geo.tif
Step three: make tiles
Turning one TIFF into a map
% gdal2tiles.py -t Vietnam -g ABQIAAAAGRHc9gAG3W9UhaqQD0pvExSVvaIHsSWhBu9dnlvZR7nHQNV82hT1HFUt0HcccsqJ2Z2e8femz7wrsw -u http://wildsong.biz/map/ 1501ANC4811_4326.tif vietnam
This creates a set of PNG files and HTML files in the output directory "vietnam". I then copy them up to the server and you can see them at http://wilkdsong.biz/map/vietnam
But I want to have multiple input files. I can merge them and then tile the output file. I also want to zoom in more than allowed by default so I make more zoom levels.
% export GKEY="ABQIAAAAGRHc9gAG3W9UhaqQD0pvExSVvaIHsSWhBu9dnlvZR7nHQNV82hT1HFUt0HcccsqJ2Z2e8femz7wrsw" % export URL="http://wildsong.biz/map/" % gdal_merge.py -o vietnam.tif -of GTiff 1501ANC4810_4326.tif 1501ANC4811_4326.tif 1501ANC4815_4326.tif' % gdal2tiles.py -t "Vietnam" -g $GKEY -u $URL -z 7-14 vietnam.tif vietnam
More data
Additional map layers from ESRI: http://www.arcgis.com/home/group.html?q=tags:ArcMap931_Base&t=group&owner=esri&title=ESRI%20Maps%20and%20Data&sortField=title&sortOrder=asc&content=all
Notes
GeoPDF is a trademark of TerraGo. The generic term is "geospatial pdf". Please in your mind translate all references to GeoPDF in this document to "geospatial pdf". Thank you.
