1

I am working on choropleth map showing particulate matter by some district of an area. The attribute table is as shown below (with geometry):-

enter image description here

The code I used is:-

library(mapview)
library(dplyr)
library(sf)
setwd("C:/Users/lenovo/Desktop/IEM")
BagmatiPM<-st_read("DistrictsPM.shp")
View(BagmatiPM)
summary(BagmatiPM)
mapview(BagmatiPM,zcol="TSPM.May", 
col.region=colorRampPalette(c("blue","green")))

And the map rendered is:- enter image description here

Its looks like CRS of map isn't same with base map. It's because though the base map is rendered but the plotted map (i.e. above map) doesn't lies on the extent of base map as shown below.

enter image description here

The spatial information of the shapefile is given below. enter image description here Thus, how can this be solved? If Reference System is to be changed how can it be done?

6
  • Please, you don't need all those packages to run this. It would be nice if we could have some sample data too.
    – Spacedman
    Commented Jun 16, 2022 at 17:35
  • Have you tried switching between the base layers using the drop-down options from the button in the top left corner?
    – Spacedman
    Commented Jun 16, 2022 at 17:41
  • @Spacedman, Yeah but it doesn't displays the basemap.
    – Walker
    Commented Jun 16, 2022 at 17:42
  • Does it show a list of basemaps, and you don't see any of them if you try selecting any of the others? Are you sure your data is in the right place? We can't tell. At least show something like summary of your data or print it so we can see the extent and the coordinate system. Can you zoom out and see anything? It could be network errors to the basemap servers.
    – Spacedman
    Commented Jun 16, 2022 at 17:46
  • If that scale says "1m / 1ft" then your data is in the wrong place and has shrunk to tiny proportions. Its now only a few meters across, and probably in an empty part of the world. You've got the wrong coordinates or coordinate system.
    – Spacedman
    Commented Jun 16, 2022 at 17:48

1 Answer 1

2

Based on the previous edit (where you shared a reproducible file and that's why reproducible questions are easy to answer):

Chitwan
# Simple feature collection with 1 feature and 16 fields
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 83.91772 ymin: 27.35449 xmax: 84.79493 ymax: 27.87813
# Projected CRS: WGS 84 / UTM zone 44N
#   dist_name shape_area shape_len created_at TSPM.Jan TSPM.Feb TSPM.Mar TSPM.Apr
# 1   Chitwan   0.206977  2.982311 2015-02-21       NA       NA   108.85   201.15
#   TSPM.May TSPM.Jun TSPM.Jul TSPM.Aug TSPM.Sep TSPM.Oct TSPM.Nov TSPM.Dec
# 1    74.35    40.77     20.3    17.13       NA       NA       NA       NA
#                      geometry
# 1 POLYGON ((84.67154 27.81635...

Ok, a polygon representing a district. Right? Look the CRS: it's UTM, so coordinates unit is meter.

Then look the bounding box: 83.92 to 84.79 and 27.34 to 27.88. Since it's UTM CRS, the size of your district is around half m^2.

That can't represent a district, it is smaller than a bed. So probably the CRS is wrong.

Just suposing the original CRS is WGS 84 geographic (it's a guess, I can't know it without getting more details of shapefile's origin):

Chitwan <- st_set_crs(Chitwan, 'epsg:4326')

mapview(Chitwan,zcol="TSPM.May", 
        col.region=colorRampPalette(c("blue","green")))

enter image description here

So probably shape_area and shape_len are in grades, not meters. Compute that property again.

2
  • @aldo_tapia-yeah that's seems the reason. Now, I have data's all transformed to Projected CRS: WGS 84 / UTM zone 44N (initially using QGIS). I worked on one map (with SRS different than above) after your comment and it lies on base map. But it have following spatial information: Geometry type: POLYGON Dimension: XY Bounding box: xmin: 83.20626 ymin: 27.32996 xmax: 83.63592 ymax: 27.76799 Geodetic CRS: GCS_NepalDD . Now will you help me know can I write a code to transform SRS from (not to set) WGS 84 / UTM zone 44N to epsg:4326 or to GCS_NepalDD.
    – Walker
    Commented Jun 16, 2022 at 22:29
  • 1
    First, st_set_crs as I showed to 4327, then st_transform to Gcs_NepalDD
    – aldo_tapia
    Commented Jun 17, 2022 at 0:38

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.