使用R语言的radiocorr函数实现Landsat影像辐射校正技术咨询

Hey there! Let's walk through your Landsat radiometric correction workflow using R's radiocorr function, break down critical technical details, and highlight common issues you might encounter along the way.

Key Technical Parameter Deep Dive #

First, let's verify that each parameter you're using aligns with Landsat specifications (since your scene is from 2006, I'm assuming it's Landsat 5 TM—adjust if it's another sensor):

• Grescale & Brescale: These are the gain and offset values for your specific band, pulled directly from the Landsat metadata file (MTL.txt). Double-check that you're using the correct values for Band 1—your current numbers look plausible for Landsat 5 TM, but always cross-reference with the scene's official metadata to avoid mismatches.
• sunelev: This is the solar elevation angle (in degrees) at the time of acquisition. Make sure this value is pulled from the metadata (look for SUN_ELEVATION)—your 43.99° seems reasonable, but confirm it's not in radians (the function expects degrees).
• satzenith: You set this to 0, which assumes the satellite is looking straight down (nadir). While Landsat scenes are near-nadir, check the metadata's SAT_ZENITH_ANGLE for the actual value—using 0 is a simplification, but might introduce minor errors if the true angle is significantly different.
• edist: The Earth-Sun distance in astronomical units (AU). This is calculated using the Julian day of acquisition, and your 1.009 AU fits for early May (since Earth is closest to the Sun in January). Again, cross-reference with metadata or calculate it via julian() in R to be precise.
• Esun: The mean solar exoatmospheric irradiance for Band 1. For Landsat 5 TM, Band 1's Esun is indeed ~1983 W/m²/μm—great job getting that right!

Workflow Validation & Improvements #

Let's go through your code step by step:

1. Reading the image: You're using readGDAL() from the rgdal package. Note that rgdal is now retired—for newer workflows, use the terra package's rast() function instead, which is more efficient and actively maintained. For example:

   library(terra)
   B1 <- rast("_X20060509_B_1.tif")
   

2. Running radiocorr: Ensure that your input B1 is a raw DN (digital number) raster—if you've already done any preprocessing (like scaling), this will throw off the correction. The function outputs a raster with apparent reflectance values, which should range between 0 and 1 (if not, you might have a parameter error).
3. Writing the output: Your path cuts off, but make sure the output directory exists (C:/Users/Documents/ Reflectance/—note the space after the slash, which might cause issues! Remove that space to avoid path errors). Use writeRaster() if switching to terra:

   writeRaster(B1.ar, "C:/Users/Documents/Reflectance/B1_reflectance.tif", overwrite=TRUE)
   

Common Pitfalls to Watch For #

• Parameter mismatches: Using gain/offset or Esun values from the wrong Landsat sensor (e.g., mixing Landsat 5 and Landsat 7 data) will produce incorrect reflectance values. Always pull parameters directly from your scene's MTL file.
• Data type errors: If your input raster is stored as a floating-point type instead of integer (raw DN), the correction calculation will be wrong. Check with typeof(B1[]) to confirm it's an integer.
• Solar angle units: If sunelev is in radians instead of degrees, the correction will be way off. Convert it with sunelev_rad * (180/pi) if needed.
• Output path issues: Spaces or special characters in your output directory path can cause writeGDAL() or writeRaster() to fail. Stick to simple, space-free paths.
• Ignoring atmospheric effects: The apparentreflectance method corrects for solar irradiance and Earth-Sun distance, but does not account for atmospheric scattering/absorption. If you need surface reflectance, you'll need additional atmospheric correction (e.g., using the 6S model or tools like LaSRC).

内容的提问来源于stack exchange，提问作者Cotimaass