1

I'd like to download in GeoTiff format, but I'm not fussy. I'd like to use Julia, and failing that Matlab, but, again, not fussy.

What's important is that I be able to download topological data and finish with a cropped dataset. The cropping should follow the administrative boundaries, as in the below picture: enter image description here

The frustrating thing here is that I did this before. I made that map of Michigan, from data I downloaded online. But I cannot for the life of me remember how I did it, nor can I find appropriate notes (all I've got is some references to R code, but I've never installed R so I know that can't be it).

1 Answer 1

0

I figured it out. I downloaded R console (macOS binaries) and then used the below code (which was heavily adapted from https://github.com/sikli/srtm_country):

(NOTE: I possibly had to install the raster package but don't recall and don't know how to uninstall/reinstall so I can be sure I'm checking correctly.)

The program is easy to use, just modify the countryName, stateName, and countyName. These map to the 0th, 1st, and 2nd levels of administrative boundaries. For instance:

  • AUT for Austria
  • FRA>>Bourgogne-Franche-Comté, for the Burgundy-Franche-Comté region in France
  • USA>>Massachusetts>>Essex for Essex County in Massachusetts, USA

Then be sure to set the boundaryLevel to represent the number of administrative levels required. So in the above example, for Austria set boundaryLevel<-0, Burgundy-Franche-Comté region in France set boundaryLevel<-1, and Essex County in Massachusetts, USA set boundaryLevel<-2.

(NOTE: Nothing bad will happen if you set all three level names, but only use the zeroth or 1st level.)

library(raster)
library(rgeos)
library(rasterVis)


#------------SETTINGS--------------

#Specify target ISO country code and path to downloaded shapefile

shp          <- shapefile("srtm/tiles.shp")       #Path to SRTM Tiles (can be found in subfolder srtm)

countryName <- "USA"
stateName <- "Kentucky"
countyName <- "Clark"


boundaryLevel <- 2

#------------EXECUTE FROM HERE--------------
# Determine if this is a 0th-level administrative boundary
if (boundaryLevel==0){
    message("Configured for zeroth-level subdivide")
    
    # Download the country
    usa0<-getData('GADM', country=countryName, level=0)
    country <- usa0
  
    downloadBorders <- country
    cropBorders <- country
    fileName <- cropBorders$NAME_0
} else {
    
    # Download the first subdivide in the country, usually a "state" or "region"
    usa1<-getData('GADM', country=countryName, level=1)
    state1 <- usa1[usa1$NAME_1 %in% c(stateName), ] 
    downloadBorders <- state1
    
    # Determine if this is a 1st-level or 2nd-level administrative boundary
    if (boundaryLevel==1) {
        message("Configured for first-level subdivide")
        
        cropBorders <- state1
        fileName <- cropBorders$NAME_1
    } else {
        message("Configured for second-level subdivide")
        # Download the second subdivide in the country, usually a "county"
        usa2<-getData('GADM', country=countryName, level=2)
        state2 <- usa2[usa2$NAME_1 %in% c(stateName), ]
        county <- state2[state2$NAME_2 %in% c(countyName), ]

        cropBorders <- county
        fileName <- cropBorders$NAME_2
    }
}

# Some console spew
message(paste("Downloading boundaries for admin area: ", fileName))

#Intersect country geometry with tile grid
intersects <- gIntersects(downloadBorders, shp, byid=T)
tiles      <- shp[intersects[,1],]


#Download tiles
message(paste("Downloading SRTM data for admin area: ", fileName))
srtm_list  <- list()
for(i in 1:length(tiles)) {
  lon <- extent(tiles[i,])[1]  + (extent(tiles[i,])[2] - extent(tiles[i,])[1]) / 2
  lat <- extent(tiles[i,])[3]  + (extent(tiles[i,])[4] - extent(tiles[i,])[3]) / 2
  
  tile <- getData('SRTM', 
                  lon=lon, 
                  lat=lat)
  
  srtm_list[[i]] <- tile
}


#Mosaic
message("Mosaicing...")
srtm_list$fun <- mean 
srtm_mosaic   <- do.call(mosaic, srtm_list)

#Crop to country borders
message("Cropping...")
srtm_crop     <- mask(srtm_mosaic, cropBorders)

message("Saving to disk...")
outFile <- paste(fileName, ".tiff", sep="")
writeRaster(srtm_crop, outFile, format="GTiff", datatype='FLT4S', overwrite=TRUE)

# Final console spew
message("")
message("**************")
message(paste("Successfully saved to: ", outFile))
message("**************")
message("")
message("")

(After downloading, the GeoTiff can be imported into a program like QGIS and then extracted to an STL with, e.g., DEMto3D.)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.