R语言中SRTM栅格重投影、分辨率提升及裁剪操作问题求助

Hey there, let's work through your DEM processing issues step by step. I see where you went wrong with the projection, and we'll fix that plus get you to your 1m resolution and cropped area in no time.

First, why your initial projection attempt didn't work #

When you used projection(r) <- "\"+proj=utm +zone=23 +datum=WGS84\"", you only modified the metadata of the raster, not the actual pixel data. The raster was still stored in geographic (lat/lon) coordinates, so the resolution value stayed in degrees (0.0002777778 degrees ≈ 30m), even though the CRS claimed to use meters. That mismatch is why manually setting a meter-based resolution gave you an empty raster—your data wasn't actually in UTM yet.

Also, your study area is in the southern hemisphere (latitude ~-21 to -22), so you need to add the +south parameter to your UTM CRS to avoid projecting to the wrong northern hemisphere zone, which would also cause empty or misaligned rasters.

Step-by-Step Solution #

We'll use the raster package (plus a bonus modern alternative with terra) to fix this properly.

1. Load your data and clean NoData values

First, load the raster and handle the SRTM NoData value (-32768):

library(raster)
# Load the original SRTM raster
r <- raster("s22_w045_1arc_v3.tif")
# Replace SRTM's NoData marker with standard NA
r[r == -32768] <- NA

2. Correctly reproject to UTM (Southern Hemisphere Zone 23)

Define the proper UTM CRS and use projectRaster() to transform the actual pixel data to this coordinate system (this converts units to meters):

# Define UTM Zone 23S CRS (critical for southern hemisphere locations!)
utm_crs <- CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs")
# Reproject the raster to UTM, preserving the original ~30m resolution
r_utm <- projectRaster(r, crs = utm_crs, method = "bilinear")

Check the attributes of r_utm now—you’ll see resolution in meters (~30.75m, matching your initial calculation) and the extent in valid UTM coordinates.

3. Resample to 1x1 meter resolution

Since you don’t care about elevation precision, you can use either bilinear interpolation or nearest neighbor. Here are two options:
Option 1: Upsample from the UTM raster

# Calculate scaling factor (original resolution / target 1m resolution)
scale_factor <- round(res(r_utm)[1] / 1)
# Upsample to 1m resolution
r_1m <- disaggregate(r_utm, fact = scale_factor, method = "bilinear")

Option 2: Reproject directly to 1m resolution (one-step shortcut)

# Skip the intermediate UTM raster and go straight to 1m resolution
r_1m <- projectRaster(r, crs = utm_crs, res = 1, method = "bilinear")

4. Crop to the 2000x2000m square centered on your target coordinates

First, get your center point in UTM coordinates (convert your original lat/lon center to UTM using a tool like QGIS or an online converter). Replace x_center and y_center with your actual values:

# Calculate the crop extent (1000m buffer on all sides of the center)
crop_extent <- extent(
  x_center - 1000,  # xmin
  x_center + 1000,  # xmax
  y_center - 1000,  # ymin
  y_center + 1000   # ymax
)
# Crop the 1m raster to this square
r_cropped <- crop(r_1m, crop_extent)

Bonus: Modern Alternative with terra #

If you want faster, more efficient processing, use the terra package (the successor to raster):

library(terra)
# Load and clean data
r <- rast("s22_w045_1arc_v3.tif")
r[r == -32768] <- NA
# Define UTM CRS
utm_crs <- "+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs"
# Reproject and resample to 1m in one step
r_1m <- project(r, utm_crs, res = 1, method = "bilinear")
# Crop to the 2000x2000m square
crop_extent <- ext(x_center - 1000, x_center + 1000, y_center - 1000, y_center + 1000)
r_cropped <- crop(r_1m, crop_extent)

Key Tips #

• Never modify raster projection metadata directly (like projection(r) <- ...)—always use projectRaster() or terra::project() to transform the actual pixel data.
• The +south parameter is non-negotiable for southern hemisphere UTM zones; omitting it will send your data to the wrong hemisphere, causing empty or misaligned rasters.
• If speed matters more than elevation smoothness, swap method = "bilinear" with method = "nearest" for faster processing.

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