I downloaded a nighttime light composite satellite image (https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html#AVSLCFC) which uses a scale for light intensity from 0 (no light) to 63 (max light). I cropped it to Chinese boundaries, using a boundary shapefile in QGIS.

Because I want to use the night time light as proxy for econ activity in Chinese cities, the ultimate goal is to have a datset of 3 columns (latitute, longitude, light intensity) for each raster of my satellite image, so that I can then allocate light aggregates to Chinese cities (I have a shapefile of all Chinese city-level admin regions).

So far, I have tried the following in R:


i <- image_read('satellite.tif') # read in the image
i_array <- as.integer(i[[1]]) # save as array (list)

I then converted the array into csv format to have a better look at it.

write.csv(china_sat_array, file = "china_sat_list.csv")

enter image description here Unfortunately, the csv file only contains zeros and sometime the value 255 (which is the value for missing data in the original tif image).

Does anyone have an idea what I'm doing wrong or an alternative how to end up with the desired three columns?

  • Which of the files on that page did you use? Do I need to download all of them? Did you try reading the data into R using the raster package?
    – Spacedman
    May 7 '19 at 15:39
  • I used the F10 Nighttime Lights Composite for 1992 only. I haven't tried the raster package yet but will try now. Is this package only for reading in the image, or also for converting it to a dataset? Many thanks
    – Timo K
    May 7 '19 at 15:48
  • Not sure what you mean by "converting it to a dataset". It reads it into an R object and you can work with that object - crop to Chinese boundaries, aggregate or average over polygons defined in another R object, sample raster values at locations, resample, filter, smooth... All the raster operations.
    – Spacedman
    May 7 '19 at 16:59

As suggested by Spacedman, the library raster has everything you need. you can load, crop and check whether there are raster values other than 0 or 255. I downloaded the night light composite of 1996 and a shapefile of Chinese boundaries here. Then, you might want to try

r <- raster('F121996.v4b_web.stable_lights.avg_vis.tif')
# if you want to plot it for visual inspection:
plot(r, col=grey(1:100/100))

# loading Chinese boundaries and cropping raster:
china <- readOGR('gadm36_CHN_shp/', 'gadm36_CHN_1')
# setting the projection of the boundaries to the raster's, so that they match:
china <- spTransform(china, CRS(as.character(r@crs)))
# cropping it:
r <- crop(r, china)
# plotting again:
plot(r, col=grey(1:100/100))
plot(china, border = 'red', col=NA, add=T)
# moreover, you can see that the cropped raster has values between 0 and 63 (last row):
# class       : RasterLayer 
# dimensions  : 4248, 7346, 31205808  (nrow, ncol, ncell)
# resolution  : 0.008333333, 0.008333333  (x, y)
# extent      : 73.55417, 134.7708, 18.1625, 53.5625  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : F121996.v4b_web.stable_lights.avg_vis 
# values      : 0, 63  (min, max)

If you still want to obtain a vector with the values of all pixels of the cropped raster, you can use the function getValues(). You can try mean(getValues(r)>0 & getValues(r)<64) to see that around 11% of the pixels have the values you were concerned about.

Finally, if your goal is to calculate zonal statistics (average nightlight intensity within Chinese boundaries -- cities in your case), you can check the function extract() from the raster package that performs this task for you. You might want to try the velox package, which also perform zonal statistics much faster than the former.

  • Dear Bruno, thanks a lot! This was a great help already, it all worked perfectly. Do you have some guidance on how to use the extract function in this case, if I want to get either the mean value or the accumulated sum of light intensity for each Chinese city? I tried finding a solution in the RDocumentation (rdocumentation.org/packages/raster/versions/2.8-19/topics/…), but couldn't get it to work. I assume I have to use the "S4 method for Raster,SpatialPolygons"? Sorry, I'm still very new to geospatial analysis. Thanks for your time and help.
    – Timo K
    May 8 '19 at 22:45

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