About Me!
This blog is about my musings and thoughts. I hope you find it useful, at most, and entertaining, at least.
Other Pages
Presence Elsewhere
Adding NAIP (and MrSID format in general) to QGIS
Date: 2015-02-12
Tags: qgis gis orthoimagery
I happened across 1-meter National Agriculture Imagery Program (NAIP) imagery on PASDA. Excitedly I downloaded the zip file for my county and uncompressed it. Inside was a shp file and a very large sid file. The shp contains what appear to be the photography tracks and not very useful for my purposes (looking at pretty pictures!). I attempted to add the sid as a raster layer only to be greeted by an error telling me that it’s an unsupported filetype.
What is this format? file was no help, just telling me it’s “data” I hit Google to attempt to figure out what this mystery format was and found out it was called MrSID (Multiresolution Seamless Image Database) an image compression format developed at Los Alamos and patented by LizardTech. (Aside: Why a government agency is giving out data in a proprietary format when suitable open alternatives exist is beyond me.)
The friendly denizens of “#qgis”:irc://irc.freenode.net/#qgis confirmed what I found and pointed me at a command line tool from LizardTech that allows one to convert MrSID files into a more open format. (While I generally don’t run proprietary code on my laptop, I sucked it up and did the conversion).
geoexpress_commandlineutils_linux/linux64/GeoExpressCLUtils-9.1.0.3982/bin/mrsidgeodecode -i ortho_1-1_1n_s_pa003_2013_1.sid -o ortho_1-1_1n_s_pa003_2013_1.tiff
BEWARE: This will increase the size of the file 15 to 20 times.
To get some idea of what’s stored in this file you can run
gdalinfo ortho_1-1_1n_s_pa003_2013_1.tiff
and see that there are 3 bands (Red, Green, and Blue) and actual size of the image.
Driver: GTiff/GeoTIFF
Files: ortho_1-1_1n_s_pa003_2013_1.tiff
Size is 59112, 56449
Coordinate System is `'
Metadata:
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,56449.0)
Upper Right (59112.0, 0.0)
Lower Right (59112.0,56449.0)
Center (29556.0,28224.5)
Band 1 Block=59112x1 Type=Byte, ColorInterp=Red
Band 2 Block=59112x1 Type=Byte, ColorInterp=Green
Band 3 Block=59112x1 Type=Byte, ColorInterp=Blue
This tiff, despite being 10GB loaded just peachy in QGIS. However, I wanted some of that disk space back, so I was told to convert it to a jpeg image. Additionally, the file does not contain it’s own projection (and probably is relying on the prj file in the same directory), so I decided to add that into the jpeg.
gdal_translate -co TILED=YES -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=85 -a_srs ortho_1-1_1n_s_pa003_2013_1.prj ortho_1-1_1n_s_pa003_2013_1.tiff ortho_1-1_1n_s_pa003_2013_1.jpeg
gdalinfo shows that we converted it, it’s the same size, and it’s the correct projection.
Driver: GTiff/GeoTIFF
Files: ortho_1-1_1n_s_pa003_2013_1.jpeg
Size is 59112, 56449
Coordinate System is:
PROJCS["NAD_1983_UTM_Zone_17N",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-81],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Metadata:
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
COMPRESSION=YCbCr JPEG
INTERLEAVE=PIXEL
SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,56449.0)
Upper Right (59112.0, 0.0)
Lower Right (59112.0,56449.0)
Center (29556.0,28224.5)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
In #qgis I was told that it’s usually prudent to precalulate overviews (which are similar to zoom-levels in a TMS). This calculates an image pyramid, but not exactly that generates images that can be used to quickly display the image when zoomed in and out without having to compute it from the base blocks every time.
gdaladdo -r gauss ortho_1-1_1n_s_pa003_2013_1.jpeg 2 4 8 16 32 64 128 256 512
Will compute overviews for half, quarter, eight, &c of the entire image (Note that these are powers of 2, as that works best because it becomes easy to calculate the size of the pyramid). gdalinfo shows that we’ve added overviews.
Driver: GTiff/GeoTIFF
Files: ortho_1-1_1n_s_pa003_2013_1.jpeg
Size is 59112, 56449
Coordinate System is:
PROJCS["NAD_1983_UTM_Zone_17N",
GEOGCS["GCS_North_American_1983",
DATUM["D_North_American_1983",
SPHEROID["GRS_1980",6378137,298.257222101]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-81],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]
Metadata:
AREA_OR_POINT=Area
TIFFTAG_XRESOLUTION=1
TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
COMPRESSION=YCbCr JPEG
INTERLEAVE=PIXEL
SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,56449.0)
Upper Right (59112.0, 0.0)
Lower Right (59112.0,56449.0)
Center (29556.0,28224.5)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Overviews: 29556x28225, 14778x14113, 7389x7057, 3695x3529, 1848x1765, 924x883, 462x442, 231x221, 116x111
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Overviews: 29556x28225, 14778x14113, 7389x7057, 3695x3529, 1848x1765, 924x883, 462x442, 231x221, 116x111
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
Overviews: 29556x28225, 14778x14113, 7389x7057, 3695x3529, 1848x1765, 924x883, 462x442, 231x221, 116x111
Now, how does this compare with the original SRID in terms of size?
| Format | Size | Relative Size |
|---|---|---|
| MrSID | 637M | 1 |
| TIFF | 9.4G | 15 |
| JPEG | 708M | 1.1 |
| JPEG (with overviews) | 1.6G | 2.5 |
So, a little larger, but natively supported by F/OSS tools such as GDAL and QGIS.
In terms of quality, I haven’t noticed a difference.
Happy Orthoimagerying!