1

I would like to create a map of India, showing the population density of India with a GeoTIFF file and a shapefile.

I took the population density from this website https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-count-rev11/data-download as a geotiff file. I choose the year 2020.

And the shapefile of India, I took from this website https://gadm.org/download_country_v3.html

library(terra)
library(raster)
library(sf) 

p_20 <- terra::rast("gpw-v4-population-count-rev11_2020_2pt5_min_tif/gpw_v4_population_count_rev11_2020_2pt5_min.tif")

IND <- sf::read_sf("gadm36_IND_shp/gadm36_IND_1.shp", quiet=TRUE)
IND <- sf::st_transform(IND , crs = "epsg:3148")

This code above works fine. I changed the EPSG, because then it shows correctly for India.

Then I would like again to a crop and later a mask and create a map. So I gave this code

p_20_crop <- terra::crop(p_20, IND)

But with this code I got that error message

Error: [crop] extents do not overlap

I think, it is because I don't have the same CRS for both files, but I cannot convert the GeoTIFF file to another CRS with sf::st_transform. Or is the problem somewhere else? Any ideas?

0

2 Answers 2

2

Use terra::project to warp a raster from one coordinate system to another:

Change the coordinate reference system

Description:

     Change the coordinate reference system ("project") of a
     SpatVector, SpatRaster or a matrix with coordinates.

Note this is a "lossy" process, involving interpolation or approximation, because you are going from one grid system to another grid system that might not match up neatly. Understand all the options if accuracy is important for you. If its just for presentation, its fine, but if you then want to sample values from locations on a warped raster, better to transform the sample locations back to the original CRS and sample from the unwarped raster.

1

You are using "epsg:3148", that is, UTM zone 48. This is not an appropriate coordinate reference system for India (see UTM zones here). You probably should not project your data to begin with. Why do you want to do that?

library(terra)
r <-rast(crs="epsg:3148")
crs(r, describe=TRUE)
#                        name authority code                            area                      extent
#1 Indian 1960 / UTM zone 48N      EPSG 3148 Cambodia; Vietnam west of 108°E 102.14, 108.00, 23.40, 8.33

As you can see, this crs is appropriate for "Cambodia; Vietnam west of 108°E", (and other areas at these longitudes), not for India.

Your Answer

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