[R-sig-Geo] converting the projection

[Michael Sumner]

Great thanks, that is good detail -but we still need the extents of both
data to be able to see where the problem is.

You gave us the proj4string for both, but we also need the extent - you can
print each object for a full summary eg.

print(bb)

print(tmp)

and with those anyone can reproduce enough (without any real data) to
recreate your experience. (I am so used to working with this stuff I forget
to remind others that print(theraster) is usually more than enough
information to work through a reproducible example)

But, a guess is - unless the projection metadata is wrong, raster is
probably telling you the truth - but it should be easy to see why there's a
problem with the print output (essentially the proj4string and extent
provide enough information for an experienced problem solver).

Cheers, Mike.



On Wed, 20 Sep 2017 at 11:05 Kyongho Son <kkyong77 at hotmail.com> wrote:

> Sorry and thanks, Mike
> I can not share all codes and data since the original data are huge.
> I copied the scripts that I used to reproject one of rasters and then make
> sure that they have same projection
>
>  Then clip a large raster with a small raster. However, I got the error
> messages, the two raster maps are not overlap.
>
> They should overlap if the projection was correct.
>
>
>
> # east US map tree phenology map
> tmp<-raster("../../obs/eMODIS_East_SOSTAvg0114_v1/sostavg0114")
>
> proj4string(tmp)
> [1] "+proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +datum=WGS84
> +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
>
> # small watershed boundary map in NY state
> bb<-raster("../../obs/gis/bb_soil_def_map.tif")
> proj4string(bb)
> [1] "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
> +towgs84=0,0,0"
>
> ex <- projectExtent(bb,  projection(tmp))
>
> > ex
> class       : RasterLayer
> dimensions  : 598, 515, 307970  (nrow, ncol, ncell)
> resolution  : 30.12012, 14.94084  (x, y)
> extent      : 6691756, 6707267, 9941473, 9950408  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0
> +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
>
>
> c<-crop(tmp,  ex)
>
> Error in .local(x, y, ...) : extents do not overlap
>
>
>
> *Kyongho Son*
> Postdoctoral Fellow
> Research Foundation of The City University of New York
>
>
>
> ------------------------------
> *From:* Michael Sumner <mdsumner at gmail.com>
> *Sent:* Tuesday, September 19, 2017 8:32 PM
> *To:* Kyongho Son; RsigGeo
>
> *Subject:* Re: [R-sig-Geo] converting the projection
> Keep on list please.
>
> Share code to demonstrate so we can see the problem.  Use projectRaster(a,
> b)  for alignment,  try crop(b,  ex,  snap = "out")  for inclusive crop.
>
> And apologies I wrote clip instead of crop earlier.
>
> Cheers,  Mike
>
> On Wed, 20 Sep 2017, 07:42 Kyongho Son <kkyong77 at hotmail.com> wrote:
>
>> I did the exact same way, but the output of the reprojection of 'a'
>> raster did not properly align with 'b' raster.
>>
>>
>>
>>
>>
>> *Kyongho Son*
>> Postdoctoral Fellow
>> Research Foundation of The City University of New York
>>
>>
>>
>> ------------------------------
>> *From:* Michael Sumner <mdsumner at gmail.com>
>> *Sent:* Tuesday, September 19, 2017 5:03 PM
>> *To:* Kyongho Son; r-sig-geo
>> *Subject:* Re: [R-sig-Geo] converting the projection
>>
>> Try projecting just the extent
>> library(raster)
>> ex <- projectExtent(a,  projection(b))
>> clip(b,  ex)
>>
>> where a,  b are your two raster.
>>
>> Cheers,  Mike
>>
>> On Wed, 20 Sep 2017, 06:18 Kyongho Son <kkyong77 at hotmail.com> wrote:
>>
>>>
>>> Hi I have two raster data sets but they have two datum and projection.
>>>
>>> I would like to convert one to the other since I would like to clip the
>>> one raster using the other raster.
>>>
>>> they should have same projections.
>>>
>>>
>>> proj4string(bb)
>>> [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
>>> proj4string(SOSTavg)
>>>
>>>  proj4string(SOSTavg)
>>> [1] "+proj=laea +lat_0=-100 +lon_0=6370997 +x_0=45 +y_0=0 +ellps=WGS84
>>> +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
>>>
>>>
>>> Thanks,
>>>
>>> Kyongho
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>> --
>>
>> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia+%3Chttps://maps.google.com/?q%3D203%2BChannel%2BHighway%2B%250D%2BKingston%2BTasmania%2B7050%2BAustralia%26entry%3Dgmail%26source%3Dg%3E&entry=gmail&source=g>
>> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>> Dr. Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> 203 Channel Highway
>> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>> Kingston Tasmania 7050 Australia
>> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
>>
>> --
>
> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia&entry=gmail&source=g>
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia+%3Chttps://maps.google.com/?q%3D203%2BChannel%2BHighway%2B%250D%2BKingston%2BTasmania%2B7050%2BAustralia%26entry%3Dgmail%26source%3Dg%3E&entry=gmail&source=g>
> Kingston Tasmania 7050 Australia
> <https://maps.google.com/?q=203+Channel+Highway+%0D+Kingston+Tasmania+7050+Australia+%3Chttps://maps.google.com/?q%3D203%2BChannel%2BHighway%2B%250D%2BKingston%2BTasmania%2B7050%2BAustralia%26entry%3Dgmail%26source%3Dg%3E&entry=gmail&source=g>
>
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]