[GRASS-user] EPSG 32633 and spatial resolution for hourly LST product from Copernicus

Nikos Alexandris nik at nikosalexandris.net
Wed Oct 3 21:55:39 PDT 2018


(For context see previous messages in this thread)

Nikos A:

>> >> g.region raster=g2_BIOPAR_LST_201606220100_GLOBE_GEO_V1.2.nc_LST
>> >> -pag w=-180 e=180

Markus M:

>> >this shifts the grid by half a cell to the east. The -a flag does
>> >not mak sense because 1) you want to force new extends, 2) there is
>> >no resolution given for use with the -a flag.

Eh.. hm, right, using `-a` is nonsensical here :-)

>> More precise, this shifts the computational's region grid by half a
>> cell to the east. Right?
>
>it shifts the raster map to the east, not the computational region.

Just a conceptual clarification, still after all these years:

`g.region` governs the computational region and `r.region` a raster
map's spatial specifications. These are related but independent.
How does `g.region` alter the raster map's location, if it can't modify
a raster map's spatial specifications?  Do you refer to any output that
will base upon the region set?

>> If the default of `-a` is "to align the region resolution to match the
>> region boundaries",

>?
>from the manual:
>"With the *-a* flag all four boundaries are adjusted to be even multiples
>of the resolution, aligning the region to the resolution supplied by the
>user."
>
>i.e. with the -a flag, the region boundaries are modified and the
>resolution is kept. The default (without -a) is to align the region
>resolution to match the region boundaries.

I misfired. The manual is clear, the `-a` flag adjusts the computational
region's boundaries to be even multiples of the spatial resolution
*supplied by the user*.

Two full examples of using `-a` to `g.region` the right way, hopefully:

```
r.in.gdal in=NETCDF:"g2_BIOPAR_LST_201606220100_GLOBE_GEO_V1.2.nc":LST out=import -o --o

WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Over-riding projection check
360 degree EW extent is exceeded by 1.00019 cells
360 degree EW extent is exceeded by 0.000192261 cells
Importing raster map <import>...
 100%
```
and
```
r.info import -g

360 degree EW extent is exceeded by 0.00019226 cells
360 degree EW extent is exceeded by 1.00019 cells
north=80.0223214291667
south=-79.9776823855556
east=179.977687153889
west=-180.022321429167
nsres=0.0446428582072328
ewres=0.0446428582072241
rows=3584
cols=8064
cells=28901376
datatype=CELL
ncats=0
```
and then setting the region via
```
g.region raster=import -ag nsres=0.0446428582072328 ewres=0.0446428582072241

360 degree EW extent is exceeded by 1.00019 cells
360 degree EW extent is exceeded by 0.00019226 cells
360 degree EW extent is exceeded by 0.00019226 cells
360 degree EW extent is exceeded by 1.00019 cells
projection=3
zone=0
n=80.0446447655684
s=-80.0000019073612
w=-180.044647149735
e=180.000004291528
nsres=0.0446428582072328
ewres=0.0446428582072241
rows=3585
cols=8065
cells=28913025
```
or via
```
g.region raster=import -ag res=0.044642858207226

360 degree EW extent is exceeded by 1.00019 cells
360 degree EW extent is exceeded by 0.00019226 cells
360 degree EW extent is exceeded by 0.00019226 cells
360 degree EW extent is exceeded by 1.00019 cells
projection=3
zone=0
n=80.0446447655562
s=-80.000001907349
w=-180.044647149742
e=180.000004291535
nsres=0.044642858207226
ewres=0.044642858207226
rows=3585
cols=8065
cells=28913025
```

all of which looks fine, except for the "exceeding" warnings.

With `-a`
```
r.in.gdal in=NETCDF:"g2_BIOPAR_LST_201606220100_GLOBE_GEO_V1.2.nc":LST out=import_with_a -o -a

WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Over-riding projection check
360 degree EW extent is exceeded by 0.000192261 cells
360 degree EW extent is exceeded by 1.00019 cells
Importing raster map <import_with_a>...
 100%
```
and
```
r.info import_with_a -g
360 degree EW extent is exceeded by 1.00019 cells
north=80.0223333333333
south=-79.9634444444444
east=179.945666666667
west=-180.022333333333
nsres=0.0446388888888889
ewres=0.0446388888888889
rows=3584
cols=8064
cells=28901376
datatype=CELL
ncats=0
```
and then setting the region
```
g.region raster=import_with_a -ag res=0.0446388888888889
360 degree EW extent is exceeded by 0.00019226 cells
360 degree EW extent is exceeded by 0.283136 cells
projection=3
zone=0
n=80.0375277777778
s=-79.9928888888889
w=-180.028638888889
e=179.984
nsres=0.0446388888888889
ewres=0.0446388888888889
rows=3585
cols=8065
cells=28913025
```
or
```
g.region raster=import_with_a -ag res=0.044642858207226
360 degree EW extent is exceeded by 1.00019 cells
360 degree EW extent is exceeded by 0.000192261 cells
projection=3
zone=0
n=80.0446447655562
s=-80.000001907349
w=-180.044647149742
e=179.955361433328
nsres=0.044642858207226
ewres=0.044642858207226
rows=3585
cols=8064
cells=28909440
```

[..]

>> won't the above command try to modify the spatial
>> resolution, of the computational region, so as to perfectly fit inside
>> the currently set, of user provided, boundaries?
>>
>> If, say, the region's boundaries are n=80 s=-80 w=-180 e=180 and the
>> raster map's rows and columns 3584 and 8064 respectively, won't the
>> above command try to adjust the region's resolution so as to fit these
>> cells inside these boundaries?
>
>If you want to "adjust the region's resolution so as to fit these cells
>inside these boundaries", you must not use the -a flag.

I got it wrong.
In short, `-a` attempts to *preserve* the (user-given?) cell resolution.

I'd like to work without the warning messages and, at the same time,
respect the spatial specifications of the input map. Given the
metadata-reported boundaries and the resolution (rows, columns) are
correct, I'd like to get these boundaries and rows by columns to govern
the definition of the computational region. Expectedly, the spatial
resolution should then be close to the one reported in the metadata.

How do we do away the warning messages? By cutting off half a cell from
the map's west limit, that is -180.0223214 + 0.044642858207226/2 = 180.
This means that the computational region will "see" only the requested
part of the map, not that a part of the map is really cut off.

Or should we cut off from the map?

Nikos
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 228 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/grass-user/attachments/20181004/8d7a44e1/attachment.sig>


More information about the grass-user mailing list