Subscribe to this thread
Home - General / All posts - Match Command Resampling
mtreglia
155 post(s)
#03-Jul-12 01:02

Hi All,

I am trying to line up raster layers for some analyses, and having some issues. Basically, I need to end up with all of my layers with the same exact spatial extent and resolution to run in MaxEnt. So, I have layers such as a 10m DEMs (and slope/aspect), 800m climate data, etc.

It seems like the "Match" command would do what I want, so that I can have everything represented with 10m pixels across the entire area. However, there are no options for the resampling procedures. The main problem comes when I am taking my 800m pixels down to 10m, because everything gets interpolated, yet I want every 10m pixel to maintain the same value as the 800m pixel (nearest neighbor sampling, right?). I just cannot justify interpolation of the climate data layers, and I realize that my inference will be somewhat limited, though the smaller-scale topography is likely more important for what I am working with, so to maintain that fine resolution.

Any suggestions?

Thanks,

Mike

dale

630 post(s)
#03-Jul-12 02:52

Mike,

Try Surface/resize to chop your 800m tiles down to 10m tiles. I understand that the output may not be co-incident with the 10m DEM and derived data.

Another thought, why not make a point grid from the 10m DEM, then transfer heights from the 800m surface to the grid. Cut and paste the grid as a surface with 10m pixels.

d

mtreglia
155 post(s)
#03-Jul-12 07:06

Dale,

Thanks for the suggestions. Unfortunately MaxEnt software really needs everything to be lined up, so the first option doesn't work.

The second option looks good, but I'm afraid it would be a cumbersomely large dataset, for most of San Diego County,CA, so a few million pixels :-\

I'm playing around in R now- there is a package called 'climstats' that seems like it does what I need. The package is downloadable from here: https://r-forge.r-project.org/R/?group_id=861 though it seems like you need to install of the dependencies separately. The command is supposedly 'spatial_sync_raster'

I'll let ya'll know what I do in the end... its too bad... Manifold ALMOST easily does what I want... but not quite

Thanks again,

Mike

PS- a cool thing about the R package 'raster' that 'climstat' relies on is that apparently it doesn't pull all of the data into memory, but rather takes the row/column type of data, and then does any processing in smaller chunks. Thus, it is not super limited on memory and such, and can handle pretty large datasets...WooHoo... we'll see I suppose

mdsumner


4,260 post(s)
#03-Jul-12 08:45

Yes, but you have to completely take on the way that raster sees the world and fit your working into that, and then with climstats you will be another step away. It also largely relies on GDAL (via rgdal) to provide that "on-demand" aspect in R. For this kind of job I'd rather use the GDAL tools directly. You are targetting a very specific tool for what could be a pretty simple job with simple tools. It might be the right choice for you of course, but if not . . . ( Python can be a better choice for scripting low-level stuff with GDAL, and you can install the whole she-bang easily with OSGeo4W).


https://github.com/mdsumner

mtreglia
155 post(s)
#03-Jul-12 09:07

Mike,

Thanks for the input. I realize what you are saying, and I thought of straight-up using GDAL, but I just haven't had a chance to really learn it that well. Have you got any suggestion of the basic steps in GDAL?

I am having to do this for a number of files, so scripting this through Python is really appealing to me... and it is something I'm trying to learn. I guess I always feel caught in the trap of needing to get things done immediately- through brute-force , but knowing that learning Python or something would make life a lot easier in the long-run. I guess its time to work my way into it.

Cheers,

Mike

mdsumner


4,260 post(s)
#03-Jul-12 10:32

Essentially it is gdalwarp with arguments -te (target extent) and -ts (target size), with the default resampling method (which is nearest neighbour). There are plenty of guides online, but feel free to email me, username at gmail.

I would like to keep this relevant to Manifold though! Got any sample data sets for us to mess around with? Let us know what you tried and how you know it's not what you want. BTW the method used by Match might be controllable with the numeric type, just a wild guess . . .


https://github.com/mdsumner

mtreglia
155 post(s)
#03-Jul-12 22:25

Thanks Mike- good call on keeping it relevant to Manifold. I'll try to post a subsample of data in the next couple of days to play with.

ColinD

2,081 post(s)
#03-Jul-12 22:14

Mike, I do a lot of modelling with Maxent at a grid size of 100m. Rather than resample, I use centroids of the grid cells and spatial overlay or transfer heights to get the data into the grid. I create the model input rasters by pasting the grid as a surface with no interpolation.


Aussie Nature Shots

mtreglia
155 post(s)
#03-Jul-12 22:29

Thanks Colin- I was interested in using 10m cells, just to maintain the small scale topography, which is likely important to my study organism. I'll might give your approach a whirl though- that is similar to what Dale suggested too. Do you mind if I e-mail you off-list? I'd be interested in hearing about your experiences with MaxEnt and what you are doing with it. My e-mail is: mlt35 [a/t] tamu.edu

Cheers,

Mike

ColinD

2,081 post(s)
#03-Jul-12 22:40

No problems, email sent.


Aussie Nature Shots

tjhb
10,094 post(s)
#04-Jul-12 03:47

Mike,

Try this to combine Match with nearest neighbour resampling.

Say your fine surface is S1 and your coarse surface S2.

Take a duplicate of S2, giving S3. Match S3 to S1. (Spatially correct, but with unwanted interpolation.)

Now in S3, open the Surface Transform dialog, and run the formula

MajValue(S2)

to restore exact original values.

Tim

dale

630 post(s)
#05-Jul-12 12:14

Nice one Tim! That's a sublimely simple solution. Eight might be getting on, but here we are, years later still learning functionality.

mtreglia
155 post(s)
#05-Jul-12 18:08

Tim,

I can't get this to work. I've attached an example file with just a DEM and a Temperature dataset. In this very small subset of my original, I've ended up with the corners of pixels matching up so its not as much of an issue here, but it is still useful to experiment with (not positive of what I did differently vs. before).

Anyway, when I do what you suggest, for MajValue(S2), it wants a window size. When I do set a window size, it seems like it is only acting on surface S2 and not doing anything to update S3. Am I missing something?

Another option I thought might work was as follows, but also ended up with interpolation (using the same notation as above):

1) match surface to get pixels to line up and of appropriate size

2) use transform to set all pixels in matched surface (S3) equal to 0 (using formula: [S3]*0)

3) Use transform to add values of original surface to matched surface ([S3]+[S2])

Again, I was surprised that this ended up with interpolation... seems kind of odd.

I have found other ways of doing what I need for now, and will work out something, but this is still a useful problem to solve.

THanks,

Mike

Attachments:
MatchProblem.map

tjhb
10,094 post(s)
#06-Jul-12 00:53

Yes, still learning--in this case at Mike's expense. I went on what I believed I knew, without fresh testing.

I remembered that the transform operators generally use interpolation—Value(), MedValue(), MinValue(), MaxValue(), AvgValue(), the filter functions, and even a plain surface name all do this. (So even a simple formula like [S] will use interpolation if resampling is required.) But I misremembered that MajValue() does not use interpolation. I was wrong. It does too, as you say.

(It's hard to read the manual description while bearing this in mind, but after all it's unlikely that one function would be wired differently from all the others. It would be great to be able to control this behaviour for every operator—perhaps with an optional extra parameter specifying a resampling method. A bit cumbersome.)

So instead we have to use our trusty SQL to restore the exact original values.

UPDATE [Us_tmax_1971_2000.01 2] -- The "Matched" surface

SET [Height (I)] = 

    Height(

        [Us_tmax_1971_2000.01]-- The original surface

        AssignCoordSys(

            NewPoint([Center X (I)][Center Y (I)]),

            CoordSys("Us_tmax_1971_2000.01 2" AS COMPONENT)

            )

        );

That works fine.

Sorry Mike.

Tim

mtreglia
155 post(s)
#06-Jul-12 17:54

Thanks a lot Tim!

I didn't see in the manual that Transform operators generally use interpolation, but its definitely good to know.

No worries- it wasn't much time on my part, and you still found an elegant solution that overcomes the problem for everybody on the forum.

Cheers,

Mike

sjvanek10 post(s)
#27-Apr-15 22:22

Hi, a related request, can I aggregate the pixels in a finer-scale surface (in this case population with 0.025 degree pixels) by adding all the pixels together from a 0.25x0.25 grid of these finer pixels, and create a coarser dataset which will then match an earlier version of the same dataset - which was only produced at 0.25x0.25 degree resolution?

I can send or post a link to the dataset if needed. Just want to see if anyone listening first...

Steve

tjhb
10,094 post(s)
#28-Apr-15 00:21

Yes you can, and it's not too hard to do, at least if the two surfaces are registered to a common grid--i.e. such that surface A exactly contains 100 cells from surface B. Is that the case?

Taking you literally, you then want to sum, for each cell in A, all the values from the 100 coincident cells in B. Is that right? (And if not by addition, then how else should values be aggregated? E.g. average, maximum, minimum, median, mode.)

sjvanek10 post(s)
#28-Apr-15 18:03

Yes that is correct. They are on exactly the same grid except I was wrong about relative size and 25 of the smaller pixels occupy each larger pixel. (5x5 and so if larger pixels are 0.25 degrees then smaaller are 0.05 degrees -- I assume something complicated happens with convergence of longitude lines but since both are lat long projection (ie unprojected) I'm willing to just go for it.

And yes, I want to sum them. Because each pixel is a population count within that area, so by summing I get the correct pop. count within the larger area.

Thanks for your prompt response. Really makes manifold a great program to have this help. It's so powerful but I would be pretty lost without the help.

Steve

tjhb
10,094 post(s)
#29-Apr-15 00:06

There are at least a few ways to do it.

A.

1. Open the fine-scale surface (which I'll call S1). Open the Surface > Transform dialog and enter the following formula.

SumValue([S1], 2)

Set Scope to All pixels, check the option at the bottom to Save result as new component, and run the transform.

This creates a new surface, the same size as S1, where each pixel gets the sum of all values in the surrounding 5x5 window (expressed here as a window size of 2--think square radius).

2. Now open the resulting new surface, S2. Open the Surface > Resize dialog. Enter Nearest Neighbour for the Method, and divide the current width and height by 5. (I'm hoping that the original dimensions are both exact multiples of 5. Otherwise it may be necessary to ensure this by adding extra edge pixels.)

This shrinks the surface S2 to the size of the coarse dataset. Using nearest neighbour means sampling just the centre pixel in each 5x5 window, discarding the rest; from step 1, this centre pixel contains the sum of all the values in this window.

That's probably the simplest way, and as fast as you'll get in Manifold. (If you have a CUDA-capable graphics card, then step 1 uses GPGPU for greater speed.)

B.

The same can be done via SQL.

1. First make a duplicate S2 of the fine-scale surface S1 or, since this is quicker, make a new surface S2 having the same dimensions and data type as S1, then use Assign Projection to give it the exact projection of S1. (Uncheck "Preserve local values, then use the Load from component button, upper right.)

2. Create and run the following query.

UPDATE [S2]

SET [Height (I)] = HeightSum([S1][X (I)][Y (I)], 2);

Now resize as for step 2 in A.

C.

Again via SQL, but using a coarse-scale surface as the target.

1. Create a duplicate S4 of the coarse scale surface S3, or the alternative as in B step 1 above, for the different surfaces.

2. Create and run the following query.

UPDATE [S4]

SET [Height (I)] = HeightSum(

    [S1],

    BoundingBox(AssignCoordSys(

        NewLine(

            NewPoint([X (I)][Y (I)]),

            NewPoint([X (I)] + 1, [Y (I)] + 1)

            ),

        CoordSys("S4" AS COMPONENT)

        ))

    );

(No need to resize after this method.)

I expect method A will be fastest and C will be slowest.

Manifold User Community Use Agreement Copyright (C) 2007-2021 Manifold Software Limited. All rights reserved.