[GRASS-user] question on i.nightlights.intercalibration code

Gabriel Cotlier gabiklm01 at gmail.com
Thu Jul 5 06:03:49 PDT 2018


Dear Nikos, Markus, and Vero,


Here is a kind of very small routine I have run on the bases of our last
chat.
I guess most of steps are covered, but I still have two questions:


1. as you can see in the message after running the code the resolution is
finally set to 1, but before this message, it still appears "exceeded by
0.999827".

The question is: does GRASS prints also the  "exceeded by 0.999827"
together with "exceeded by 1 cells" message because one corresponds to the
old resolution and the other to the new one, and it is letting know about
the change or is because other resason?

2. How can I run this same code for many layers instead with one with only
one using the command r.in.gdal -a input=PathToMap output=MapName, since
r.in.gdal -a takes as input one layer at time?

Thanks a lot.

Best regards,

Gabriel


+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Step 1) iport the original raster data layers with the code:
============================================================

> r.in.gdal -a input=C:\Users\Gabriel\Documents\grassdata\light\F121996.tif
output=F121996

360 degree EW extent is exceeded by 0.999827 cells
360 degree EW extent is exceeded by 0.999827 cells
360 degree EW extent is exceeded by 0.999827 cells
360 degree EW extent is exceeded by 1 cells


 use of r.in.gdal -a helps to adjust resolution for lat long maps fixing
the resolution from 0.008333333300000 to 0.008333333333333, i.e. exactly 30
arc-seconds.


Step 2) show metadate to check:
==============================

> r.info map=F121996


360 degree EW extent is exceeded by 1 cells
360 degree EW extent is exceeded by 0.999827 cells
 +----------------------------------------------------------
------------------+
 | Map:      F121996 at PERMANENT              Date: Thu Jul 05 09:19:34 2018
  |
 | Mapset:   PERMANENT                      Login of Creator: Gabriel
   |
 | Location: lightsLoc1
   |
 | DataBase: C:\Users\Gabriel\Documents\grassdata\light
   |
 | Title:    F121996
  |
 | Timestamp: none
  |
 |----------------------------------------------------------
------------------|
 |
  |
 |   Type of Map:  raster               Number of Categories: 0
   |
 |   Data Type:    CELL
   |
 |   Rows:         16801
  |
 |   Columns:      43201
  |
 |   Total Cells:  725820001
  |
 |        Projection: Latitude-Longitude
  |
 |            N:  75:00:15N    S:  65:00:15S   Res: 0:00:30
   |
 |            E: 180:00:15E    W: 180:00:15W   Res: 0:00:30
   |
 |   Range of data:    min = 0  max = 63
  |
 |
  |
 |   Data Description:
  |
 |    generated by r.in.gdal
  |
 |
  |
 |   Comments:
  |
 |    r.in.gdal -a input="C:\Users\Gabriel\Documents\grassdata\light\F1219\
 |
 |    96.tif" output="F121996" memory=300 offset=0 num_digits=0
   |
 |
  |
 +----------------------------------------------------------------------------+


Step 2: run i.nightlights
====================================================================

> i.nightlights.intercalibration image=F121996 suffix=c model=elvidge2014
-t

|i Inter-satellite calibration of DMSP-OLS Nighttime Stable Lights
WARNING: Operating on current region

|> Calibrating average visible Digital Number values
Regression coefficients: (-0.0959, 1.2727, -0.004) | Associated R^2: 0.9319
360 degree EW extent is exceeded by 0.999827 cells
360 degree EW extent is exceeded by 1 cells


Step 3: seting the region to one of the imported maps with the code:
==================================================================

>r.region -a map=F121996

360 degree EW extent is exceeded by 0.999827 cells
360 degree EW extent is exceeded by 1 cells
360 degree EW extent is exceeded by 1 cells
r.region complete.


On Wed, Jun 27, 2018 at 4:55 AM, Nikos Alexandris <nik at nikosalexandris.net>
wrote:

> * Markus Metz <markus.metz.giswork at gmail.com> [2018-06-26 14:40:25 +0200]:
>
>
> On Tue, Jun 26, 2018 at 9:53 AM, Nikos Alexandris <nik at nikosalexandris.net
>> >
>> wrote:
>>
>>>
>>> * Markus Metz <markus.metz.giswork at gmail.com> [2018-06-25 08:29:45
>>> +0200]:
>>>
>>> [..]
>>>
>>> The resolution is a bit wrong, it is 0.008333333300000 but should be
>>>> 0.008333333333333, i.e. exactly 30 arc-seconds. This can be solved with
>>>>
>>> the
>>
>>> -a flag of r.in.gdal, or after import with r.region -a.
>>>>
>>>> The message
>>>>
>>>> 360 degree EW extent is exceeded by 0.999827 cells
>>>>>
>>>>
>>>>
>>>> will change to
>>>>
>>>> 360 degree EW extent is exceeded by 1 cells
>>>>>
>>>>
>>>>
>>>> but will not go away, because 360 degree EW extent is exceeded in the
>>>>
>>> input
>>
>>> data, the first and last column cover the same geographical area. You can
>>>> change your current region to chop of e.g. the first column: set the
>>>>
>>> region
>>
>>> to the raster, then modify the current region with g.region w=179:59:45W
>>>>
>>> -p
>>
>>> and use this region for further processing.
>>>>
>>>
>
> I guess this is worth being documented in the manual of the add-on.
>>>
>>
>> This is a universal problem applying to various raster data in latlong.
>> The
>> first issue, 30 sec represented as 0.008333333300000 instead of
>> 0.008333333333333 is solved by r.in.gdal -a. The second problem, this
>> extra
>> column responsible that 360 degree EW extent is exceeded by 1 cell can be
>> solved by setting the current region accordingly. This is also a universal
>> problem. Maybe the manual of r.in.gdal could include a hint about how to
>> do
>> this. Generally, users are encouraged to inspect the output of r.info
>> after
>> importing raster data to check if everything is as expected.
>>
>
> Danke Markus. We should collect some of these hints.
>
> (
> Oh! I took on a refreshing read-tour: a pixel's anchor point,
> pixel-is-area,
> pixel-is-center, center-to-center extent model, precision of pixel size.
> Yet, I ended up reading about many more issues that people have to deal
> with.
>
> Some interesting entries:
>
> - https://gis.stackexchange.com/q/122670/5256
> - https://gis.stackexchange.com/a/243050/5256
> )
>
> Would
>>> it also make sense to let the module attempt to perform this
>>> "correction"?
>>>
>>
>> If you refer to i.nightlights.intercalibration, I would say no, because
>> it
>> is a more general issue not restricted to DMSP-OLS nightlight data. Even
>> more general, the current region as set by the user is used for raster
>> processing (with a very few exceptions).
>>
>
> Yes, I refer to it. Understood, I won't touch upon this within the
> module.
>
> Nikos
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20180705/1f95a640/attachment.html>


More information about the grass-user mailing list