maiac_vApr2018 = "/data-coco/ECHO_PM/MAIAC_20180423/"
example_tile = "h04v04/"
example_year = 2016
maiac_hdf_files = list.files(paste0(maiac_vApr2018, example_tile, example_year), pattern = ".hdf", full.names = TRUE)
maiac_latlon_files = list.files(maiac_vApr2018, pattern = ".hdf", full.names = TRUE)Open the first file in the specified tile and year. It happens to be an Aqua AOT file.
suppressMessages(library(gdalUtils))
gdalinfo(maiac_hdf_files[[1]])## [1] "Driver: HDF4/Hierarchical Data Format Release 4"
## [2] "Files: /data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf"
## [3] "Size is 512, 512"
## [4] "Coordinate System is `'"
## [5] "Metadata:"
## [6] " HDFEOSVersion=HDFEOS_V2.12"
## [7] "Subdatasets:"
## [8] " SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:Optical_Depth_047"
## [9] " SUBDATASET_1_DESC=[1200x1200] Optical_Depth_047 grid1km (16-bit integer)"
## [10] " SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:Optical_Depth_055"
## [11] " SUBDATASET_2_DESC=[1200x1200] Optical_Depth_055 grid1km (16-bit integer)"
## [12] " SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:AOT_Uncertainty"
## [13] " SUBDATASET_3_DESC=[1200x1200] AOT_Uncertainty grid1km (16-bit integer)"
## [14] " SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:FineModeFraction"
## [15] " SUBDATASET_4_DESC=[1200x1200] FineModeFraction grid1km (16-bit integer)"
## [16] " SUBDATASET_5_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:Column_WV"
## [17] " SUBDATASET_5_DESC=[1200x1200] Column_WV grid1km (16-bit integer)"
## [18] " SUBDATASET_6_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:Injection_Height"
## [19] " SUBDATASET_6_DESC=[1200x1200] Injection_Height grid1km (32-bit floating-point)"
## [20] " SUBDATASET_7_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:AOT_QA"
## [21] " SUBDATASET_7_DESC=[1200x1200] AOT_QA grid1km (16-bit unsigned integer)"
## [22] " SUBDATASET_8_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid1km:AOT_MODEL"
## [23] " SUBDATASET_8_DESC=[1200x1200] AOT_MODEL grid1km (8-bit unsigned integer)"
## [24] " SUBDATASET_9_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid5km:cosSZA"
## [25] " SUBDATASET_9_DESC=[240x240] cosSZA grid5km (16-bit integer)"
## [26] " SUBDATASET_10_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid5km:cosVZA"
## [27] " SUBDATASET_10_DESC=[240x240] cosVZA grid5km (16-bit integer)"
## [28] " SUBDATASET_11_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid5km:RelAZ"
## [29] " SUBDATASET_11_DESC=[240x240] RelAZ grid5km (16-bit integer)"
## [30] " SUBDATASET_12_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid5km:Scattering_Angle"
## [31] " SUBDATASET_12_DESC=[240x240] Scattering_Angle grid5km (16-bit integer)"
## [32] " SUBDATASET_13_NAME=HDF4_EOS:EOS_GRID:\"/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf\":grid5km:Glint_Angle"
## [33] " SUBDATASET_13_DESC=[240x240] Glint_Angle grid5km (16-bit integer)"
## [34] "Corner Coordinates:"
## [35] "Upper Left ( 0.0, 0.0)"
## [36] "Lower Left ( 0.0, 512.0)"
## [37] "Upper Right ( 512.0, 0.0)"
## [38] "Lower Right ( 512.0, 512.0)"
## [39] "Center ( 256.0, 256.0)"sd_prefix = "HDF4_EOS:EOS_GRID:"
sd_layer = ":grid1km:AOT_QA"
subdataset_path = paste0(sd_prefix, maiac_hdf_files[[1]], sd_layer)
gdalinfo(subdataset_path)## [1] "Driver: HDF4Image/HDF4 Dataset"
## [2] "Files: /data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf"
## [3] "Size is 1200, 1200"
## [4] "Coordinate System is:"
## [5] "PROJCS[\"unnamed\","
## [6] " GEOGCS[\"Unknown datum based upon the WGS 84 ellipsoid\","
## [7] " DATUM[\"Not specified (based on WGS 84 spheroid)\","
## [8] " SPHEROID[\"WGS 84\",6378137,298.257223563,"
## [9] " AUTHORITY[\"EPSG\",\"7030\"]]],"
## [10] " PRIMEM[\"Greenwich\",0],"
## [11] " UNIT[\"degree\",0.0174532925199433]],"
## [12] " PROJECTION[\"Albers_Conic_Equal_Area\"],"
## [13] " PARAMETER[\"standard_parallel_1\",1145915590.261647],"
## [14] " PARAMETER[\"standard_parallel_2\",3437746770.784939],"
## [15] " PARAMETER[\"latitude_of_center\",1317802928.800894],"
## [16] " PARAMETER[\"longitude_of_center\",-5500394833.255903],"
## [17] " PARAMETER[\"false_easting\",0],"
## [18] " PARAMETER[\"false_northing\",0],"
## [19] " UNIT[\"Meter\",1]]"
## [20] "Origin = (800000.000000000000000,2400000.000000000000000)"
## [21] "Pixel Size = (1000.000000000000000,-1000.000000000000000)"
## [22] "Metadata:"
## [23] " data description=Bits\tDefinition"
## [24] "0-2 Cloud Mask"
## [25] " 000 --- Undefined"
## [26] " 001 --- Clear"
## [27] " 010 --- Possible Cloudy"
## [28] " 011 --- Cloudy "
## [29] " 101 --- Cloud shadow"
## [30] " 110 --- Fire hotspot"
## [31] " 111 --- Water Sediments"
## [32] "3-4 Land Water Snow/ice Mask"
## [33] " 00 --- Land"
## [34] " 01 --- Water"
## [35] " 10 --- Snow"
## [36] " 11 --- Ice"
## [37] "5-7 Adjacency Mask"
## [38] " 000 --- Normal condition"
## [39] " 001 --- Adjacent to cloud"
## [40] " 010 --- Surrounded by more than 8 cloudy pixels"
## [ reached getOption("max.print") -- omitted 41 entries ]See the prefix and postfix data around the path of the subdatasets above
suppressMessages(library(rgdal))
maiac_aot_qa = readGDAL(subdataset_path)## HDF4_EOS:EOS_GRID:/data-coco/ECHO_PM/MAIAC_20180423/h04v04/2016/MAIACAAOT.h04v04.20160011825.hdf:grid1km:AOT_QA has GDAL driver HDF4Image
## and has 1200 rows and 1200 columnsmaiac_aot_qa## Object of class SpatialGridDataFrame
## Object of class SpatialGrid
## Grid topology:
## cellcentre.offset cellsize cells.dim
## x 800500 1000 1200
## y 1200500 1000 1200
## SpatialPoints:
## x y
## [1,] 800500 2399500
## [2,] 801500 2399500
## [3,] 802500 2399500
## [4,] 803500 2399500
## [5,] 804500 2399500
## [6,] 805500 2399500
## [7,] 806500 2399500
## [8,] 807500 2399500
## [9,] 808500 2399500
## [10,] 809500 2399500
## [11,] 810500 2399500
## [12,] 811500 2399500
## [13,] 812500 2399500
## [14,] 813500 2399500
## [15,] 814500 2399500
## [16,] 815500 2399500
## [17,] 816500 2399500
## [18,] 817500 2399500
## [19,] 818500 2399500
## [20,] 819500 2399500
## [ reached getOption("max.print") -- omitted 1439980 rows ]
## Coordinate Reference System (CRS) arguments: +proj=aea
## +lat_1=1145915590.261647 +lat_2=3437746770.784939
## +lat_0=1317802928.800894 +lon_0=-5500394833.255903 +x_0=0 +y_0=0
## +ellps=WGS84 +units=m +no_defs
##
## Data summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 3.00 3.00 6.76 11.00 32770.00 1class(maiac_aot_qa)## [1] "SpatialGridDataFrame"
## attr(,"package")
## [1] "sp"proj4string(maiac_aot_qa)## [1] "+proj=aea +lat_1=1145915590.261647 +lat_2=3437746770.784939 +lat_0=1317802928.800894 +lon_0=-5500394833.255903 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"Based on this CRS, it places a raster that should cover Michigan to South Carolina and places it between New Zealand and Antarctica.
suppressMessages(library(mapview))
mapview(maiac_aot_qa, at =seq(1,12,1))## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 1440000
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 1440000 'Change the proj4string to match the description from Yujie’s email:
You can always read the projection information from any MAIAC HDF file, in the global attribute fields.
For North America, I used Albers equal area projection, the lat of the first standard parallel is 20 degree north and the second is 60 degree north; the central meridian is 96 degree west.
The lat of the projection origin is 23 degree north. The Upper left corner of the map is -4000000 meter east and 7200000 meter north.
You can calculate the boundary of each tile using the tile dimension (1200kmx1200km) and pixel size(1000m).Old proj4string has the parallels and meridians in what initially looks to be meters, except the circumfrence of the Earth is 40m meters, so the values are too large: +proj=aea +lat_1=1145915590.261647 +lat_2=3437746770.784939 +lat_0=1317802928.800894 +lon_0=-5500394833.255903 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
Change the values to match his email, using degrees for units of the standard parallels and meridians: +proj=aea +lat_1=20 +lat_2=60 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
proj4string(maiac_aot_qa) <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"## Warning in ReplProj4string(obj, CRS(value)): A new CRS was assigned to an object with an existing CRS:
## +proj=aea +lat_1=1145915590.261647 +lat_2=3437746770.784939 +lat_0=1317802928.800894 +lon_0=-5500394833.255903 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
## without reprojecting.
## For reprojection, use function spTransformmaiac_aot_qa## Object of class SpatialGridDataFrame
## Object of class SpatialGrid
## Grid topology:
## cellcentre.offset cellsize cells.dim
## x 800500 1000 1200
## y 1200500 1000 1200
## SpatialPoints:
## x y
## [1,] 800500 2399500
## [2,] 801500 2399500
## [3,] 802500 2399500
## [4,] 803500 2399500
## [5,] 804500 2399500
## [6,] 805500 2399500
## [7,] 806500 2399500
## [8,] 807500 2399500
## [9,] 808500 2399500
## [10,] 809500 2399500
## [11,] 810500 2399500
## [12,] 811500 2399500
## [13,] 812500 2399500
## [14,] 813500 2399500
## [15,] 814500 2399500
## [16,] 815500 2399500
## [17,] 816500 2399500
## [18,] 817500 2399500
## [19,] 818500 2399500
## [20,] 819500 2399500
## [ reached getOption("max.print") -- omitted 1439980 rows ]
## Coordinate Reference System (CRS) arguments: +proj=aea +lat_1=20
## +lat_2=60 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=WGS84 +units=m
## +no_defs
##
## Data summary:
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 3.00 3.00 6.76 11.00 32770.00 1mapview(maiac_aot_qa, at =seq(1,12,1))## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 1440000
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 1440000 '