adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

classic Classic list List threaded Threaded
22 messages Options
12
Reply | Threaded
Open this post in threaded view
|

adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Dear all,

Im trying to overlap these points:
https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0

and a wrld_simpl object:
library(maptools)
data(wrld_simpl)

Over this raster layer
https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0

This rastr comes from a .nc file without a reference system. The author of
that .nc file gave me the following data about the .nc.

The projection is *Lambert conformal conic* projection
CEN_LAT = 38.0
CEN_LON = -100.0
TRUELAT1 = 25.
TRUELAT2 = 45.

However, despite i have gone through many sites in the internet, i cant
figure it out:

a) if that is all the data i need to set a reference system for my points
and the wrld_simp object.

b) how to change a typical CRS object with such data

Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")

Where do i enter the TRUELAT and CENLAT values?
Are there any site that explains easily what the fields in the CRS mean and
how to change them?

Thanks in advance.
--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Alex Mandel-2
On 02/22/2016 09:50 AM, Agus Camacho wrote:

> Dear all,
>
> Im trying to overlap these points:
> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>
> and a wrld_simpl object:
> library(maptools)
> data(wrld_simpl)
>
> Over this raster layer
> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>
> This rastr comes from a .nc file without a reference system. The author of
> that .nc file gave me the following data about the .nc.
>
> The projection is *Lambert conformal conic* projection
> CEN_LAT = 38.0
> CEN_LON = -100.0
> TRUELAT1 = 25.
> TRUELAT2 = 45.
>
> However, despite i have gone through many sites in the internet, i cant
> figure it out:
>
> a) if that is all the data i need to set a reference system for my points
> and the wrld_simp object.
>
> b) how to change a typical CRS object with such data
>
> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>
> Where do i enter the TRUELAT and CENLAT values?
> Are there any site that explains easily what the fields in the CRS mean and
> how to change them?
>
> Thanks in advance.
>

https://github.com/OSGeo/proj.4/wiki/GenParms
https://trac.osgeo.org/proj/wiki/GenParms

I believe:
+lat_0  = CEN_LAT   Latitude of origin
+lat_1  = TRUELAT1   Latitude of first standard parallel
+lat_2  = TRUELAT2   Latitude of second standard parallel
+lon_0  = CEN_LON   Central meridian

proj strings are defined by the proj4 libary. It's website listed above
and the associated mailing lists or gis stackexchange would be the
places to get help on it.
https://lists.osgeo.org/mailman/listinfo/metacrs

It often helps to browse similar projections on
http://spatialreference.org/
http://epsg.io/

Enjoy,
Alex

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Thanks Alex, but the locations still fall in the sea when i plot them using
your recommended Solution. I looked at the sites you proposed and they have
other values for lat_1, lat_0, etc..

2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:

> On 02/22/2016 09:50 AM, Agus Camacho wrote:
> > Dear all,
> >
> > Im trying to overlap these points:
> >
> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
> >
> > and a wrld_simpl object:
> > library(maptools)
> > data(wrld_simpl)
> >
> > Over this raster layer
> >
> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
> >
> > This rastr comes from a .nc file without a reference system. The author
> of
> > that .nc file gave me the following data about the .nc.
> >
> > The projection is *Lambert conformal conic* projection
> > CEN_LAT = 38.0
> > CEN_LON = -100.0
> > TRUELAT1 = 25.
> > TRUELAT2 = 45.
> >
> > However, despite i have gone through many sites in the internet, i cant
> > figure it out:
> >
> > a) if that is all the data i need to set a reference system for my points
> > and the wrld_simp object.
> >
> > b) how to change a typical CRS object with such data
> >
> > Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
> >
> > Where do i enter the TRUELAT and CENLAT values?
> > Are there any site that explains easily what the fields in the CRS mean
> and
> > how to change them?
> >
> > Thanks in advance.
> >
>
> https://github.com/OSGeo/proj.4/wiki/GenParms
> https://trac.osgeo.org/proj/wiki/GenParms
>
> I believe:
> +lat_0  = CEN_LAT   Latitude of origin
> +lat_1  = TRUELAT1   Latitude of first standard parallel
> +lat_2  = TRUELAT2   Latitude of second standard parallel
> +lon_0  = CEN_LON   Central meridian
>
> proj strings are defined by the proj4 libary. It's website listed above
> and the associated mailing lists or gis stackexchange would be the
> places to get help on it.
> https://lists.osgeo.org/mailman/listinfo/metacrs
>
> It often helps to browse similar projections on
> http://spatialreference.org/
> http://epsg.io/
>
> Enjoy,
> Alex
>
>


--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
but its ok, ill ask in another mail list.

2016-02-22 11:19 GMT-07:00 Agus Camacho <[hidden email]>:

> Thanks Alex, but the locations still fall in the sea when i plot them
> using your recommended Solution. I looked at the sites you proposed and
> they have other values for lat_1, lat_0, etc..
>
> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>
>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>> > Dear all,
>> >
>> > Im trying to overlap these points:
>> >
>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>> >
>> > and a wrld_simpl object:
>> > library(maptools)
>> > data(wrld_simpl)
>> >
>> > Over this raster layer
>> >
>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>> >
>> > This rastr comes from a .nc file without a reference system. The author
>> of
>> > that .nc file gave me the following data about the .nc.
>> >
>> > The projection is *Lambert conformal conic* projection
>> > CEN_LAT = 38.0
>> > CEN_LON = -100.0
>> > TRUELAT1 = 25.
>> > TRUELAT2 = 45.
>> >
>> > However, despite i have gone through many sites in the internet, i cant
>> > figure it out:
>> >
>> > a) if that is all the data i need to set a reference system for my
>> points
>> > and the wrld_simp object.
>> >
>> > b) how to change a typical CRS object with such data
>> >
>> > Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>> >
>> > Where do i enter the TRUELAT and CENLAT values?
>> > Are there any site that explains easily what the fields in the CRS mean
>> and
>> > how to change them?
>> >
>> > Thanks in advance.
>> >
>>
>> https://github.com/OSGeo/proj.4/wiki/GenParms
>> https://trac.osgeo.org/proj/wiki/GenParms
>>
>> I believe:
>> +lat_0  = CEN_LAT   Latitude of origin
>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>> +lon_0  = CEN_LON   Central meridian
>>
>> proj strings are defined by the proj4 libary. It's website listed above
>> and the associated mailing lists or gis stackexchange would be the
>> places to get help on it.
>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>
>> It often helps to browse similar projections on
>> http://spatialreference.org/
>> http://epsg.io/
>>
>> Enjoy,
>> Alex
>>
>>
>
>
> --
> Agustín Camacho Guerrero.
> Doutor em Zoologia.
> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
> Biociências, USP.
> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
> São Paulo - SP, CEP: 05508-090, Brasil.
>



--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Alex Mandel-2
In reply to this post by agus
I made an attempt at it too. Investigating the original data, I'm not
sure that the projection information supplied is correct for the data
linked. When I load up the data in a unprojected space, the coordinates
don't look at all similar to any Lambert projected data I have, they
actually look like Lat/Lon in some unprojected coordinate system,
perhaps a different spheroid than expected.

-Alex

On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:

> Hi
>
> I tried to make it work but I had to give up. I wanted to reproject the Lamberth conformal conic coordinates to long-lat but it didn't work.
>
> Perhaps someone can see what I did wrong. Here is what I did (data in R binary format and figure in png format both attached):
>
> library(raster)
> library(maptools)
> data(wrld_simpl)
>
> r <- raster("raster.grd")
> projection(r)
> ## > NA
>
> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = ",")
> coordinates(uro) <- ~lon+lat
>
> ## Set projections for the 3 data sets
>
> ## Lamberth's confocal conic projection with given parameters
> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 +lat_2=45.0 +ellps=WGS84"
> projection(r)
>
> ## Assume that lon, lat are geographical coordinates (degrees decimal)
> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>
> ## wrld_simpl is in geographical coordinates
> proj4string(wrld_simpl)
>
> ## Make figure in png format
> ## Of course plotting data with 2 different projections will give
> ## some distortions
> pdf("uro.png")
>
> plot(r)
> points(uro)
> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of 'r'
>
> dev.off()
>
> extent(r)
> ## class       : Extent
> ## xmin        : -131.4368
> ## xmax        : -68.56323
> ## ymin        : 12.35567
> ## ymax        : 50.26619
>
> ## Reproject the raster to long-lat
> ## This doesn't work (collapsed domain)
> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
>
> ## Because
>> extent(rp)
> ## class       : Extent
> ## xmin        : -100.0015
> ## xmax        : -99.68557
> ## ymin        : 37.70658
> ## ymax        : 38.00046
>
> ## Save data in R binary format
> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>
>
> Yours sincerely / Med venlig hilsen
>
> Frede Aakmann Tøgersen
> Specialist, M.Sc., Ph.D.
> Plant Performance & Modeling
>
> Technology & Service Solutions
> T +45 9730 5135
> M +45 2547 6050
> [hidden email]
> http://www.vestas.com
>
> Company reg. name: Vestas Wind Systems A/S
> This e-mail is subject to our e-mail disclaimer statement.
> Please refer to www.vestas.com/legal/notice
> If you have received this e-mail in error please contact the sender.
>
>
>
> -----Original Message-----
> From: R-sig-Geo [mailto:[hidden email]] On Behalf Of Agus Camacho
> Sent: 22. februar 2016 19:20
> To: [hidden email]
> Cc: r-sig-geo
> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
>
> Thanks Alex, but the locations still fall in the sea when i plot them using
> your recommended Solution. I looked at the sites you proposed and they have
> other values for lat_1, lat_0, etc..
>
> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>
>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>> Dear all,
>>>
>>> Im trying to overlap these points:
>>>
>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>
>>> and a wrld_simpl object:
>>> library(maptools)
>>> data(wrld_simpl)
>>>
>>> Over this raster layer
>>>
>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>
>>> This rastr comes from a .nc file without a reference system. The author
>> of
>>> that .nc file gave me the following data about the .nc.
>>>
>>> The projection is *Lambert conformal conic* projection
>>> CEN_LAT = 38.0
>>> CEN_LON = -100.0
>>> TRUELAT1 = 25.
>>> TRUELAT2 = 45.
>>>
>>> However, despite i have gone through many sites in the internet, i cant
>>> figure it out:
>>>
>>> a) if that is all the data i need to set a reference system for my points
>>> and the wrld_simp object.
>>>
>>> b) how to change a typical CRS object with such data
>>>
>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>
>>> Where do i enter the TRUELAT and CENLAT values?
>>> Are there any site that explains easily what the fields in the CRS mean
>> and
>>> how to change them?
>>>
>>> Thanks in advance.
>>>
>>
>> https://github.com/OSGeo/proj.4/wiki/GenParms
>> https://trac.osgeo.org/proj/wiki/GenParms
>>
>> I believe:
>> +lat_0  = CEN_LAT   Latitude of origin
>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>> +lon_0  = CEN_LON   Central meridian
>>
>> proj strings are defined by the proj4 libary. It's website listed above
>> and the associated mailing lists or gis stackexchange would be the
>> places to get help on it.
>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>
>> It often helps to browse similar projections on
>> http://spatialreference.org/
>> http://epsg.io/
>>
>> Enjoy,
>> Alex
>>
>>
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Roger Bivand
Administrator
On Tue, 23 Feb 2016, Alex Mandel wrote:

> I made an attempt at it too. Investigating the original data, I'm not
> sure that the projection information supplied is correct for the data
> linked. When I load up the data in a unprojected space, the coordinates
> don't look at all similar to any Lambert projected data I have, they
> actually look like Lat/Lon in some unprojected coordinate system,
> perhaps a different spheroid than expected.

Does anyone have a link to the original data? Is is possible that this is
the General Oblique Transformation used by modellers - that is something
that feels like longlat but is recentred and oblique? Example at the very
end of my GEOSTAT talk last year (slides 81-83):

http://geostat-course.org/system/files/geostat_talk_150817.pdf

Roger

>
> -Alex
>
> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>> Hi
>>
>> I tried to make it work but I had to give up. I wanted to reproject the Lamberth conformal conic coordinates to long-lat but it didn't work.
>>
>> Perhaps someone can see what I did wrong. Here is what I did (data in R binary format and figure in png format both attached):
>>
>> library(raster)
>> library(maptools)
>> data(wrld_simpl)
>>
>> r <- raster("raster.grd")
>> projection(r)
>> ## > NA
>>
>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = ",")
>> coordinates(uro) <- ~lon+lat
>>
>> ## Set projections for the 3 data sets
>>
>> ## Lamberth's confocal conic projection with given parameters
>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 +lat_2=45.0 +ellps=WGS84"
>> projection(r)
>>
>> ## Assume that lon, lat are geographical coordinates (degrees decimal)
>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>
>> ## wrld_simpl is in geographical coordinates
>> proj4string(wrld_simpl)
>>
>> ## Make figure in png format
>> ## Of course plotting data with 2 different projections will give
>> ## some distortions
>> pdf("uro.png")
>>
>> plot(r)
>> points(uro)
>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of 'r'
>>
>> dev.off()
>>
>> extent(r)
>> ## class       : Extent
>> ## xmin        : -131.4368
>> ## xmax        : -68.56323
>> ## ymin        : 12.35567
>> ## ymax        : 50.26619
>>
>> ## Reproject the raster to long-lat
>> ## This doesn't work (collapsed domain)
>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
>>
>> ## Because
>>> extent(rp)
>> ## class       : Extent
>> ## xmin        : -100.0015
>> ## xmax        : -99.68557
>> ## ymin        : 37.70658
>> ## ymax        : 38.00046
>>
>> ## Save data in R binary format
>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>
>>
>> Yours sincerely / Med venlig hilsen
>>
>> Frede Aakmann Tøgersen
>> Specialist, M.Sc., Ph.D.
>> Plant Performance & Modeling
>>
>> Technology & Service Solutions
>> T +45 9730 5135
>> M +45 2547 6050
>> [hidden email]
>> http://www.vestas.com
>>
>> Company reg. name: Vestas Wind Systems A/S
>> This e-mail is subject to our e-mail disclaimer statement.
>> Please refer to www.vestas.com/legal/notice
>> If you have received this e-mail in error please contact the sender.
>>
>>
>>
>> -----Original Message-----
>> From: R-sig-Geo [mailto:[hidden email]] On Behalf Of Agus Camacho
>> Sent: 22. februar 2016 19:20
>> To: [hidden email]
>> Cc: r-sig-geo
>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a reference system implicit in a .nc file
>>
>> Thanks Alex, but the locations still fall in the sea when i plot them using
>> your recommended Solution. I looked at the sites you proposed and they have
>> other values for lat_1, lat_0, etc..
>>
>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>
>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>> Dear all,
>>>>
>>>> Im trying to overlap these points:
>>>>
>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>>
>>>> and a wrld_simpl object:
>>>> library(maptools)
>>>> data(wrld_simpl)
>>>>
>>>> Over this raster layer
>>>>
>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>>
>>>> This rastr comes from a .nc file without a reference system. The author
>>> of
>>>> that .nc file gave me the following data about the .nc.
>>>>
>>>> The projection is *Lambert conformal conic* projection
>>>> CEN_LAT = 38.0
>>>> CEN_LON = -100.0
>>>> TRUELAT1 = 25.
>>>> TRUELAT2 = 45.
>>>>
>>>> However, despite i have gone through many sites in the internet, i cant
>>>> figure it out:
>>>>
>>>> a) if that is all the data i need to set a reference system for my points
>>>> and the wrld_simp object.
>>>>
>>>> b) how to change a typical CRS object with such data
>>>>
>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>
>>>> Where do i enter the TRUELAT and CENLAT values?
>>>> Are there any site that explains easily what the fields in the CRS mean
>>> and
>>>> how to change them?
>>>>
>>>> Thanks in advance.
>>>>
>>>
>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>
>>> I believe:
>>> +lat_0  = CEN_LAT   Latitude of origin
>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>> +lon_0  = CEN_LON   Central meridian
>>>
>>> proj strings are defined by the proj4 libary. It's website listed above
>>> and the associated mailing lists or gis stackexchange would be the
>>> places to get help on it.
>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>
>>> It often helps to browse similar projections on
>>> http://spatialreference.org/
>>> http://epsg.io/
>>>
>>> Enjoy,
>>> Alex
>>>
>>>
>>
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Michael Sumner-2
On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:

> On Tue, 23 Feb 2016, Alex Mandel wrote:
>
> > I made an attempt at it too. Investigating the original data, I'm not
> > sure that the projection information supplied is correct for the data
> > linked. When I load up the data in a unprojected space, the coordinates
> > don't look at all similar to any Lambert projected data I have, they
> > actually look like Lat/Lon in some unprojected coordinate system,
> > perhaps a different spheroid than expected.
>
> Does anyone have a link to the original data? Is is possible that this is
> the General Oblique Transformation used by modellers - that is something
> that feels like longlat but is recentred and oblique? Example at the very
> end of my GEOSTAT talk last year (slides 81-83):
>
> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>
> Roger
>
>
For what it is worth, the General Oblique Transformation is not the only
example - it's very common for modellers to have a mesh that has the
"mostly-properties" of a projection, but is not actually describable with
standard transform + affine parameters. The main cases that I've seen are
polar stereographic, equal area or oblique Mercator. Often they really are
simple transforms and you can reconstruct without loss, but it's not
usually possible to tell without exploration. It's an interesting
dis-connect to see code that builds a mesh with certain properties, then
only stores longitudes and latitudes - when it could be done with standard
tools and be stored and used much more efficiently.

(I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
terminology conflated in this context too. )

I'm also interested to explore the original data.

Cheers, Mike.



> >
> > -Alex
> >
> > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
> >> Hi
> >>
> >> I tried to make it work but I had to give up. I wanted to reproject the
> Lamberth conformal conic coordinates to long-lat but it didn't work.
> >>
> >> Perhaps someone can see what I did wrong. Here is what I did (data in R
> binary format and figure in png format both attached):
> >>
> >> library(raster)
> >> library(maptools)
> >> data(wrld_simpl)
> >>
> >> r <- raster("raster.grd")
> >> projection(r)
> >> ## > NA
> >>
> >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = ",")
> >> coordinates(uro) <- ~lon+lat
> >>
> >> ## Set projections for the 3 data sets
> >>
> >> ## Lamberth's confocal conic projection with given parameters
> >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 +lat_2=45.0
> +ellps=WGS84"
> >> projection(r)
> >>
> >> ## Assume that lon, lat are geographical coordinates (degrees decimal)
> >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
> >>
> >> ## wrld_simpl is in geographical coordinates
> >> proj4string(wrld_simpl)
> >>
> >> ## Make figure in png format
> >> ## Of course plotting data with 2 different projections will give
> >> ## some distortions
> >> pdf("uro.png")
> >>
> >> plot(r)
> >> points(uro)
> >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of 'r'
> >>
> >> dev.off()
> >>
> >> extent(r)
> >> ## class       : Extent
> >> ## xmin        : -131.4368
> >> ## xmax        : -68.56323
> >> ## ymin        : 12.35567
> >> ## ymax        : 50.26619
> >>
> >> ## Reproject the raster to long-lat
> >> ## This doesn't work (collapsed domain)
> >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
> +ellps=WGS84 +towgs84=0,0,0")
> >>
> >> ## Because
> >>> extent(rp)
> >> ## class       : Extent
> >> ## xmin        : -100.0015
> >> ## xmax        : -99.68557
> >> ## ymin        : 37.70658
> >> ## ymax        : 38.00046
> >>
> >> ## Save data in R binary format
> >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
> >>
> >>
> >> Yours sincerely / Med venlig hilsen
> >>
> >> Frede Aakmann Tøgersen
> >> Specialist, M.Sc., Ph.D.
> >> Plant Performance & Modeling
> >>
> >> Technology & Service Solutions
> >> T +45 9730 5135
> >> M +45 2547 6050
> >> [hidden email]
> >> http://www.vestas.com
> >>
> >> Company reg. name: Vestas Wind Systems A/S
> >> This e-mail is subject to our e-mail disclaimer statement.
> >> Please refer to www.vestas.com/legal/notice
> >> If you have received this e-mail in error please contact the sender.
> >>
> >>
> >>
> >> -----Original Message-----
> >> From: R-sig-Geo [mailto:[hidden email]] On Behalf Of
> Agus Camacho
> >> Sent: 22. februar 2016 19:20
> >> To: [hidden email]
> >> Cc: r-sig-geo
> >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
> reference system implicit in a .nc file
> >>
> >> Thanks Alex, but the locations still fall in the sea when i plot them
> using
> >> your recommended Solution. I looked at the sites you proposed and they
> have
> >> other values for lat_1, lat_0, etc..
> >>
> >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
> >>
> >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
> >>>> Dear all,
> >>>>
> >>>> Im trying to overlap these points:
> >>>>
> >>>
> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
> >>>>
> >>>> and a wrld_simpl object:
> >>>> library(maptools)
> >>>> data(wrld_simpl)
> >>>>
> >>>> Over this raster layer
> >>>>
> >>>
> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
> >>>>
> >>>> This rastr comes from a .nc file without a reference system. The
> author
> >>> of
> >>>> that .nc file gave me the following data about the .nc.
> >>>>
> >>>> The projection is *Lambert conformal conic* projection
> >>>> CEN_LAT = 38.0
> >>>> CEN_LON = -100.0
> >>>> TRUELAT1 = 25.
> >>>> TRUELAT2 = 45.
> >>>>
> >>>> However, despite i have gone through many sites in the internet, i
> cant
> >>>> figure it out:
> >>>>
> >>>> a) if that is all the data i need to set a reference system for my
> points
> >>>> and the wrld_simp object.
> >>>>
> >>>> b) how to change a typical CRS object with such data
> >>>>
> >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
> >>>>
> >>>> Where do i enter the TRUELAT and CENLAT values?
> >>>> Are there any site that explains easily what the fields in the CRS
> mean
> >>> and
> >>>> how to change them?
> >>>>
> >>>> Thanks in advance.
> >>>>
> >>>
> >>> https://github.com/OSGeo/proj.4/wiki/GenParms
> >>> https://trac.osgeo.org/proj/wiki/GenParms
> >>>
> >>> I believe:
> >>> +lat_0  = CEN_LAT   Latitude of origin
> >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
> >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
> >>> +lon_0  = CEN_LON   Central meridian
> >>>
> >>> proj strings are defined by the proj4 libary. It's website listed above
> >>> and the associated mailing lists or gis stackexchange would be the
> >>> places to get help on it.
> >>> https://lists.osgeo.org/mailman/listinfo/metacrs
> >>>
> >>> It often helps to browse similar projections on
> >>> http://spatialreference.org/
> >>> http://epsg.io/
> >>>
> >>> Enjoy,
> >>> Alex
> >>>
> >>>
> >>
> >>
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 91 00
> e-mail: [hidden email]
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
> http://depsy.org/person/434412
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

--
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
more giant crab this time.

The creator of the .nc file also looked at this webpage:
http://www.pkrc.net/wrf-lambert.html
It seemed like the right proj4 string might be this one:

+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
    +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs

However this projection also does not allow me to adequately plot the
locations on the raster.

Here is the .nc file. it contains several layers.

https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0



2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:

> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>
> > On Tue, 23 Feb 2016, Alex Mandel wrote:
> >
> > > I made an attempt at it too. Investigating the original data, I'm not
> > > sure that the projection information supplied is correct for the data
> > > linked. When I load up the data in a unprojected space, the coordinates
> > > don't look at all similar to any Lambert projected data I have, they
> > > actually look like Lat/Lon in some unprojected coordinate system,
> > > perhaps a different spheroid than expected.
> >
> > Does anyone have a link to the original data? Is is possible that this is
> > the General Oblique Transformation used by modellers - that is something
> > that feels like longlat but is recentred and oblique? Example at the very
> > end of my GEOSTAT talk last year (slides 81-83):
> >
> > http://geostat-course.org/system/files/geostat_talk_150817.pdf
> >
> > Roger
> >
> >
> For what it is worth, the General Oblique Transformation is not the only
> example - it's very common for modellers to have a mesh that has the
> "mostly-properties" of a projection, but is not actually describable with
> standard transform + affine parameters. The main cases that I've seen are
> polar stereographic, equal area or oblique Mercator. Often they really are
> simple transforms and you can reconstruct without loss, but it's not
> usually possible to tell without exploration. It's an interesting
> dis-connect to see code that builds a mesh with certain properties, then
> only stores longitudes and latitudes - when it could be done with standard
> tools and be stored and used much more efficiently.
>
> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
> terminology conflated in this context too. )
>
> I'm also interested to explore the original data.
>
> Cheers, Mike.
>
>
>
> > >
> > > -Alex
> > >
> > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
> > >> Hi
> > >>
> > >> I tried to make it work but I had to give up. I wanted to reproject
> the
> > Lamberth conformal conic coordinates to long-lat but it didn't work.
> > >>
> > >> Perhaps someone can see what I did wrong. Here is what I did (data in
> R
> > binary format and figure in png format both attached):
> > >>
> > >> library(raster)
> > >> library(maptools)
> > >> data(wrld_simpl)
> > >>
> > >> r <- raster("raster.grd")
> > >> projection(r)
> > >> ## > NA
> > >>
> > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = ",")
> > >> coordinates(uro) <- ~lon+lat
> > >>
> > >> ## Set projections for the 3 data sets
> > >>
> > >> ## Lamberth's confocal conic projection with given parameters
> > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 +lat_2=45.0
> > +ellps=WGS84"
> > >> projection(r)
> > >>
> > >> ## Assume that lon, lat are geographical coordinates (degrees decimal)
> > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
> > >>
> > >> ## wrld_simpl is in geographical coordinates
> > >> proj4string(wrld_simpl)
> > >>
> > >> ## Make figure in png format
> > >> ## Of course plotting data with 2 different projections will give
> > >> ## some distortions
> > >> pdf("uro.png")
> > >>
> > >> plot(r)
> > >> points(uro)
> > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of 'r'
> > >>
> > >> dev.off()
> > >>
> > >> extent(r)
> > >> ## class       : Extent
> > >> ## xmin        : -131.4368
> > >> ## xmax        : -68.56323
> > >> ## ymin        : 12.35567
> > >> ## ymax        : 50.26619
> > >>
> > >> ## Reproject the raster to long-lat
> > >> ## This doesn't work (collapsed domain)
> > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
> > +ellps=WGS84 +towgs84=0,0,0")
> > >>
> > >> ## Because
> > >>> extent(rp)
> > >> ## class       : Extent
> > >> ## xmin        : -100.0015
> > >> ## xmax        : -99.68557
> > >> ## ymin        : 37.70658
> > >> ## ymax        : 38.00046
> > >>
> > >> ## Save data in R binary format
> > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
> > >>
> > >>
> > >> Yours sincerely / Med venlig hilsen
> > >>
> > >> Frede Aakmann Tøgersen
> > >> Specialist, M.Sc., Ph.D.
> > >> Plant Performance & Modeling
> > >>
> > >> Technology & Service Solutions
> > >> T +45 9730 5135
> > >> M +45 2547 6050
> > >> [hidden email]
> > >> http://www.vestas.com
> > >>
> > >> Company reg. name: Vestas Wind Systems A/S
> > >> This e-mail is subject to our e-mail disclaimer statement.
> > >> Please refer to www.vestas.com/legal/notice
> > >> If you have received this e-mail in error please contact the sender.
> > >>
> > >>
> > >>
> > >> -----Original Message-----
> > >> From: R-sig-Geo [mailto:[hidden email]] On Behalf Of
> > Agus Camacho
> > >> Sent: 22. februar 2016 19:20
> > >> To: [hidden email]
> > >> Cc: r-sig-geo
> > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
> > reference system implicit in a .nc file
> > >>
> > >> Thanks Alex, but the locations still fall in the sea when i plot them
> > using
> > >> your recommended Solution. I looked at the sites you proposed and they
> > have
> > >> other values for lat_1, lat_0, etc..
> > >>
> > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
> > >>
> > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
> > >>>> Dear all,
> > >>>>
> > >>>> Im trying to overlap these points:
> > >>>>
> > >>>
> >
> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
> > >>>>
> > >>>> and a wrld_simpl object:
> > >>>> library(maptools)
> > >>>> data(wrld_simpl)
> > >>>>
> > >>>> Over this raster layer
> > >>>>
> > >>>
> >
> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
> > >>>>
> > >>>> This rastr comes from a .nc file without a reference system. The
> > author
> > >>> of
> > >>>> that .nc file gave me the following data about the .nc.
> > >>>>
> > >>>> The projection is *Lambert conformal conic* projection
> > >>>> CEN_LAT = 38.0
> > >>>> CEN_LON = -100.0
> > >>>> TRUELAT1 = 25.
> > >>>> TRUELAT2 = 45.
> > >>>>
> > >>>> However, despite i have gone through many sites in the internet, i
> > cant
> > >>>> figure it out:
> > >>>>
> > >>>> a) if that is all the data i need to set a reference system for my
> > points
> > >>>> and the wrld_simp object.
> > >>>>
> > >>>> b) how to change a typical CRS object with such data
> > >>>>
> > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
> > >>>>
> > >>>> Where do i enter the TRUELAT and CENLAT values?
> > >>>> Are there any site that explains easily what the fields in the CRS
> > mean
> > >>> and
> > >>>> how to change them?
> > >>>>
> > >>>> Thanks in advance.
> > >>>>
> > >>>
> > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
> > >>> https://trac.osgeo.org/proj/wiki/GenParms
> > >>>
> > >>> I believe:
> > >>> +lat_0  = CEN_LAT   Latitude of origin
> > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
> > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
> > >>> +lon_0  = CEN_LON   Central meridian
> > >>>
> > >>> proj strings are defined by the proj4 libary. It's website listed
> above
> > >>> and the associated mailing lists or gis stackexchange would be the
> > >>> places to get help on it.
> > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
> > >>>
> > >>> It often helps to browse similar projections on
> > >>> http://spatialreference.org/
> > >>> http://epsg.io/
> > >>>
> > >>> Enjoy,
> > >>> Alex
> > >>>
> > >>>
> > >>
> > >>
> > >
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > [hidden email]
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> > --
> > Roger Bivand
> > Department of Economics, Norwegian School of Economics,
> > Helleveien 30, N-5045 Bergen, Norway.
> > voice: +47 55 95 93 55; fax +47 55 95 91 00
> > e-mail: [hidden email]
> > http://orcid.org/0000-0003-2392-6140
> > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
> > http://depsy.org/person/434412
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Roger Bivand
Administrator
On Tue, 23 Feb 2016, Agus Camacho wrote:

> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
> more giant crab this time.
>
> The creator of the .nc file also looked at this webpage:
> http://www.pkrc.net/wrf-lambert.html
> It seemed like the right proj4 string might be this one:
>
> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>    +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>
> However this projection also does not allow me to adequately plot the
> locations on the raster.
>
> Here is the .nc file. it contains several layers.
>
> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>
Thanks for the link. The proj4 string is missing +units=km (the spherical
axes are certainly not in metres), but is wrong anyway.

This says:

> gI <- GDALinfo("results_us_future_output_none_0.nc")
Error in .local(.Object, ...) :
   No UNIDATA NC_GLOBAL:Conventions attribute

The values of the coordinate variables are:

library(ncdf4)
r <- nc_open("results_us_future_output_none_0.nc")
pts <- cbind(c(ncvar_get(r, "lon")), c(ncvar_get(r, "lat")))

> summary(pts)
        V1                V2
  Min.   :-151.30   Min.   :12.36
  1st Qu.:-120.44   1st Qu.:25.45
  Median :-100.00   Median :35.72
  Mean   :-100.00   Mean   :35.71
  3rd Qu.: -79.56   3rd Qu.:45.91
  Max.   : -48.70   Max.   :58.33

which are obviously not a grid:

pts1 <- SpatialPoints(pts)
plot(pts1)
# gridded(pts1) <- TRUE fails miserably

This is definitely not lcc, irrespective of what the originator says.

Mike - any forensics?

Roger

>
>
> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>
>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>
>>> On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>
>>>> I made an attempt at it too. Investigating the original data, I'm not
>>>> sure that the projection information supplied is correct for the data
>>>> linked. When I load up the data in a unprojected space, the coordinates
>>>> don't look at all similar to any Lambert projected data I have, they
>>>> actually look like Lat/Lon in some unprojected coordinate system,
>>>> perhaps a different spheroid than expected.
>>>
>>> Does anyone have a link to the original data? Is is possible that this is
>>> the General Oblique Transformation used by modellers - that is something
>>> that feels like longlat but is recentred and oblique? Example at the very
>>> end of my GEOSTAT talk last year (slides 81-83):
>>>
>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>
>>> Roger
>>>
>>>
>> For what it is worth, the General Oblique Transformation is not the only
>> example - it's very common for modellers to have a mesh that has the
>> "mostly-properties" of a projection, but is not actually describable with
>> standard transform + affine parameters. The main cases that I've seen are
>> polar stereographic, equal area or oblique Mercator. Often they really are
>> simple transforms and you can reconstruct without loss, but it's not
>> usually possible to tell without exploration. It's an interesting
>> dis-connect to see code that builds a mesh with certain properties, then
>> only stores longitudes and latitudes - when it could be done with standard
>> tools and be stored and used much more efficiently.
>>
>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>> terminology conflated in this context too. )
>>
>> I'm also interested to explore the original data.
>>
>> Cheers, Mike.
>>
>>
>>
>>>>
>>>> -Alex
>>>>
>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>>> Hi
>>>>>
>>>>> I tried to make it work but I had to give up. I wanted to reproject
>> the
>>> Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>>>
>>>>> Perhaps someone can see what I did wrong. Here is what I did (data in
>> R
>>> binary format and figure in png format both attached):
>>>>>
>>>>> library(raster)
>>>>> library(maptools)
>>>>> data(wrld_simpl)
>>>>>
>>>>> r <- raster("raster.grd")
>>>>> projection(r)
>>>>> ## > NA
>>>>>
>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep = ",")
>>>>> coordinates(uro) <- ~lon+lat
>>>>>
>>>>> ## Set projections for the 3 data sets
>>>>>
>>>>> ## Lamberth's confocal conic projection with given parameters
>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 +lat_2=45.0
>>> +ellps=WGS84"
>>>>> projection(r)
>>>>>
>>>>> ## Assume that lon, lat are geographical coordinates (degrees decimal)
>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>>>
>>>>> ## wrld_simpl is in geographical coordinates
>>>>> proj4string(wrld_simpl)
>>>>>
>>>>> ## Make figure in png format
>>>>> ## Of course plotting data with 2 different projections will give
>>>>> ## some distortions
>>>>> pdf("uro.png")
>>>>>
>>>>> plot(r)
>>>>> points(uro)
>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of 'r'
>>>>>
>>>>> dev.off()
>>>>>
>>>>> extent(r)
>>>>> ## class       : Extent
>>>>> ## xmin        : -131.4368
>>>>> ## xmax        : -68.56323
>>>>> ## ymin        : 12.35567
>>>>> ## ymax        : 50.26619
>>>>>
>>>>> ## Reproject the raster to long-lat
>>>>> ## This doesn't work (collapsed domain)
>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>> +ellps=WGS84 +towgs84=0,0,0")
>>>>>
>>>>> ## Because
>>>>>> extent(rp)
>>>>> ## class       : Extent
>>>>> ## xmin        : -100.0015
>>>>> ## xmax        : -99.68557
>>>>> ## ymin        : 37.70658
>>>>> ## ymax        : 38.00046
>>>>>
>>>>> ## Save data in R binary format
>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>>>
>>>>>
>>>>> Yours sincerely / Med venlig hilsen
>>>>>
>>>>> Frede Aakmann Tøgersen
>>>>> Specialist, M.Sc., Ph.D.
>>>>> Plant Performance & Modeling
>>>>>
>>>>> Technology & Service Solutions
>>>>> T +45 9730 5135
>>>>> M +45 2547 6050
>>>>> [hidden email]
>>>>> http://www.vestas.com
>>>>>
>>>>> Company reg. name: Vestas Wind Systems A/S
>>>>> This e-mail is subject to our e-mail disclaimer statement.
>>>>> Please refer to www.vestas.com/legal/notice
>>>>> If you have received this e-mail in error please contact the sender.
>>>>>
>>>>>
>>>>>
>>>>> -----Original Message-----
>>>>> From: R-sig-Geo [mailto:[hidden email]] On Behalf Of
>>> Agus Camacho
>>>>> Sent: 22. februar 2016 19:20
>>>>> To: [hidden email]
>>>>> Cc: r-sig-geo
>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
>>> reference system implicit in a .nc file
>>>>>
>>>>> Thanks Alex, but the locations still fall in the sea when i plot them
>>> using
>>>>> your recommended Solution. I looked at the sites you proposed and they
>>> have
>>>>> other values for lat_1, lat_0, etc..
>>>>>
>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>>>
>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>>>>> Dear all,
>>>>>>>
>>>>>>> Im trying to overlap these points:
>>>>>>>
>>>>>>
>>>
>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>>>>>
>>>>>>> and a wrld_simpl object:
>>>>>>> library(maptools)
>>>>>>> data(wrld_simpl)
>>>>>>>
>>>>>>> Over this raster layer
>>>>>>>
>>>>>>
>>>
>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>>>>>
>>>>>>> This rastr comes from a .nc file without a reference system. The
>>> author
>>>>>> of
>>>>>>> that .nc file gave me the following data about the .nc.
>>>>>>>
>>>>>>> The projection is *Lambert conformal conic* projection
>>>>>>> CEN_LAT = 38.0
>>>>>>> CEN_LON = -100.0
>>>>>>> TRUELAT1 = 25.
>>>>>>> TRUELAT2 = 45.
>>>>>>>
>>>>>>> However, despite i have gone through many sites in the internet, i
>>> cant
>>>>>>> figure it out:
>>>>>>>
>>>>>>> a) if that is all the data i need to set a reference system for my
>>> points
>>>>>>> and the wrld_simp object.
>>>>>>>
>>>>>>> b) how to change a typical CRS object with such data
>>>>>>>
>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>>>>
>>>>>>> Where do i enter the TRUELAT and CENLAT values?
>>>>>>> Are there any site that explains easily what the fields in the CRS
>>> mean
>>>>>> and
>>>>>>> how to change them?
>>>>>>>
>>>>>>> Thanks in advance.
>>>>>>>
>>>>>>
>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>>>>
>>>>>> I believe:
>>>>>> +lat_0  = CEN_LAT   Latitude of origin
>>>>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>>>> +lon_0  = CEN_LON   Central meridian
>>>>>>
>>>>>> proj strings are defined by the proj4 libary. It's website listed
>> above
>>>>>> and the associated mailing lists or gis stackexchange would be the
>>>>>> places to get help on it.
>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>>>>
>>>>>> It often helps to browse similar projections on
>>>>>> http://spatialreference.org/
>>>>>> http://epsg.io/
>>>>>>
>>>>>> Enjoy,
>>>>>> Alex
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>> --
>>> Roger Bivand
>>> Department of Economics, Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>> e-mail: [hidden email]
>>> http://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>> http://depsy.org/person/434412
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> --
>> Dr. Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> 203 Channel Highway
>> Kingston Tasmania 7050 Australia
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

dschneiderch
In reply to this post by agus
This looks like WRF <http://www.wrf-model.org/index.php> data. I just dealt
with this.
The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
+a=6370000 +b=6370000 +units=m

+proj=lcc which is usually what wrf is run with.
The tricky part is:
+lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
because every WRF run is different (the WRF Preprocessing System optimizes
the projection for the domain).
and then there is probably no shift so you need(?) +x_0=0 +y_0=0

This gives:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
+a=6370000 +b=6370000 +units=m +no_defs

But, wrf users like to give out lat and  long so you need to assign it:
+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs

and then reproject the lat/long to lcc coordinates using this string:
+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
+a=6370000 +b=6370000 +units=m +no_defs

One word of caution, make sure you received the correct coordinates. Some
variables are run cell center while some are run at cell edge. It looks
like from your .nc file it was made by your collaborator so I assume they
are right.

That said, another word of caution, I found that the XLAT and XLONG
variables from WRF output aren't very precise. There is a "geogrid" file
from the preprocessing system that has the domain corners, resolution, nrow
and ncol from which you can make a better grid using the native projection
system (in my case it was a 4km grid). I suggest you try to get those.

I hope this helps... I have to run but wanted to save people too much head
scratching. I can get you running with more help tonight if you need.
Dominik


On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
wrote:

> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
> more giant crab this time.
>
> The creator of the .nc file also looked at this webpage:
> http://www.pkrc.net/wrf-lambert.html
> It seemed like the right proj4 string might be this one:
>
> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>     +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>
> However this projection also does not allow me to adequately plot the
> locations on the raster.
>
> Here is the .nc file. it contains several layers.
>
>
> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>
>
>
> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>
> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
> >
> > > On Tue, 23 Feb 2016, Alex Mandel wrote:
> > >
> > > > I made an attempt at it too. Investigating the original data, I'm not
> > > > sure that the projection information supplied is correct for the data
> > > > linked. When I load up the data in a unprojected space, the
> coordinates
> > > > don't look at all similar to any Lambert projected data I have, they
> > > > actually look like Lat/Lon in some unprojected coordinate system,
> > > > perhaps a different spheroid than expected.
> > >
> > > Does anyone have a link to the original data? Is is possible that this
> is
> > > the General Oblique Transformation used by modellers - that is
> something
> > > that feels like longlat but is recentred and oblique? Example at the
> very
> > > end of my GEOSTAT talk last year (slides 81-83):
> > >
> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf
> > >
> > > Roger
> > >
> > >
> > For what it is worth, the General Oblique Transformation is not the only
> > example - it's very common for modellers to have a mesh that has the
> > "mostly-properties" of a projection, but is not actually describable with
> > standard transform + affine parameters. The main cases that I've seen are
> > polar stereographic, equal area or oblique Mercator. Often they really
> are
> > simple transforms and you can reconstruct without loss, but it's not
> > usually possible to tell without exploration. It's an interesting
> > dis-connect to see code that builds a mesh with certain properties, then
> > only stores longitudes and latitudes - when it could be done with
> standard
> > tools and be stored and used much more efficiently.
> >
> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
> > terminology conflated in this context too. )
> >
> > I'm also interested to explore the original data.
> >
> > Cheers, Mike.
> >
> >
> >
> > > >
> > > > -Alex
> > > >
> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
> > > >> Hi
> > > >>
> > > >> I tried to make it work but I had to give up. I wanted to reproject
> > the
> > > Lamberth conformal conic coordinates to long-lat but it didn't work.
> > > >>
> > > >> Perhaps someone can see what I did wrong. Here is what I did (data
> in
> > R
> > > binary format and figure in png format both attached):
> > > >>
> > > >> library(raster)
> > > >> library(maptools)
> > > >> data(wrld_simpl)
> > > >>
> > > >> r <- raster("raster.grd")
> > > >> projection(r)
> > > >> ## > NA
> > > >>
> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
> ",")
> > > >> coordinates(uro) <- ~lon+lat
> > > >>
> > > >> ## Set projections for the 3 data sets
> > > >>
> > > >> ## Lamberth's confocal conic projection with given parameters
> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0 +lat_2=45.0
> > > +ellps=WGS84"
> > > >> projection(r)
> > > >>
> > > >> ## Assume that lon, lat are geographical coordinates (degrees
> decimal)
> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
> > > >>
> > > >> ## wrld_simpl is in geographical coordinates
> > > >> proj4string(wrld_simpl)
> > > >>
> > > >> ## Make figure in png format
> > > >> ## Of course plotting data with 2 different projections will give
> > > >> ## some distortions
> > > >> pdf("uro.png")
> > > >>
> > > >> plot(r)
> > > >> points(uro)
> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
> 'r'
> > > >>
> > > >> dev.off()
> > > >>
> > > >> extent(r)
> > > >> ## class       : Extent
> > > >> ## xmin        : -131.4368
> > > >> ## xmax        : -68.56323
> > > >> ## ymin        : 12.35567
> > > >> ## ymax        : 50.26619
> > > >>
> > > >> ## Reproject the raster to long-lat
> > > >> ## This doesn't work (collapsed domain)
> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
> > > +ellps=WGS84 +towgs84=0,0,0")
> > > >>
> > > >> ## Because
> > > >>> extent(rp)
> > > >> ## class       : Extent
> > > >> ## xmin        : -100.0015
> > > >> ## xmax        : -99.68557
> > > >> ## ymin        : 37.70658
> > > >> ## ymax        : 38.00046
> > > >>
> > > >> ## Save data in R binary format
> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
> > > >>
> > > >>
> > > >> Yours sincerely / Med venlig hilsen
> > > >>
> > > >> Frede Aakmann Tøgersen
> > > >> Specialist, M.Sc., Ph.D.
> > > >> Plant Performance & Modeling
> > > >>
> > > >> Technology & Service Solutions
> > > >> T +45 9730 5135
> > > >> M +45 2547 6050
> > > >> [hidden email]
> > > >> http://www.vestas.com
> > > >>
> > > >> Company reg. name: Vestas Wind Systems A/S
> > > >> This e-mail is subject to our e-mail disclaimer statement.
> > > >> Please refer to www.vestas.com/legal/notice
> > > >> If you have received this e-mail in error please contact the sender.
> > > >>
> > > >>
> > > >>
> > > >> -----Original Message-----
> > > >> From: R-sig-Geo [mailto:[hidden email]] On Behalf
> Of
> > > Agus Camacho
> > > >> Sent: 22. februar 2016 19:20
> > > >> To: [hidden email]
> > > >> Cc: r-sig-geo
> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
> > > reference system implicit in a .nc file
> > > >>
> > > >> Thanks Alex, but the locations still fall in the sea when i plot
> them
> > > using
> > > >> your recommended Solution. I looked at the sites you proposed and
> they
> > > have
> > > >> other values for lat_1, lat_0, etc..
> > > >>
> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
> > > >>
> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
> > > >>>> Dear all,
> > > >>>>
> > > >>>> Im trying to overlap these points:
> > > >>>>
> > > >>>
> > >
> >
> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
> > > >>>>
> > > >>>> and a wrld_simpl object:
> > > >>>> library(maptools)
> > > >>>> data(wrld_simpl)
> > > >>>>
> > > >>>> Over this raster layer
> > > >>>>
> > > >>>
> > >
> >
> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
> > > >>>>
> > > >>>> This rastr comes from a .nc file without a reference system. The
> > > author
> > > >>> of
> > > >>>> that .nc file gave me the following data about the .nc.
> > > >>>>
> > > >>>> The projection is *Lambert conformal conic* projection
> > > >>>> CEN_LAT = 38.0
> > > >>>> CEN_LON = -100.0
> > > >>>> TRUELAT1 = 25.
> > > >>>> TRUELAT2 = 45.
> > > >>>>
> > > >>>> However, despite i have gone through many sites in the internet, i
> > > cant
> > > >>>> figure it out:
> > > >>>>
> > > >>>> a) if that is all the data i need to set a reference system for my
> > > points
> > > >>>> and the wrld_simp object.
> > > >>>>
> > > >>>> b) how to change a typical CRS object with such data
> > > >>>>
> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
> > > >>>>
> > > >>>> Where do i enter the TRUELAT and CENLAT values?
> > > >>>> Are there any site that explains easily what the fields in the CRS
> > > mean
> > > >>> and
> > > >>>> how to change them?
> > > >>>>
> > > >>>> Thanks in advance.
> > > >>>>
> > > >>>
> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
> > > >>> https://trac.osgeo.org/proj/wiki/GenParms
> > > >>>
> > > >>> I believe:
> > > >>> +lat_0  = CEN_LAT   Latitude of origin
> > > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
> > > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
> > > >>> +lon_0  = CEN_LON   Central meridian
> > > >>>
> > > >>> proj strings are defined by the proj4 libary. It's website listed
> > above
> > > >>> and the associated mailing lists or gis stackexchange would be the
> > > >>> places to get help on it.
> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
> > > >>>
> > > >>> It often helps to browse similar projections on
> > > >>> http://spatialreference.org/
> > > >>> http://epsg.io/
> > > >>>
> > > >>> Enjoy,
> > > >>> Alex
> > > >>>
> > > >>>
> > > >>
> > > >>
> > > >
> > > > _______________________________________________
> > > > R-sig-Geo mailing list
> > > > [hidden email]
> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
> > > --
> > > Roger Bivand
> > > Department of Economics, Norwegian School of Economics,
> > > Helleveien 30, N-5045 Bergen, Norway.
> > > voice: +47 55 95 93 55; fax +47 55 95 91 00
> > > e-mail: [hidden email]
> > > http://orcid.org/0000-0003-2392-6140
> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
> > > http://depsy.org/person/434412
> > > _______________________________________________
> > > R-sig-Geo mailing list
> > > [hidden email]
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> > --
> > Dr. Michael Sumner
> > Software and Database Engineer
> > Australian Antarctic Division
> > 203 Channel Highway
> > Kingston Tasmania 7050 Australia
> >
> >         [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
>
> --
> Agustín Camacho Guerrero.
> Doutor em Zoologia.
> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
> Biociências, USP.
> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
> São Paulo - SP, CEP: 05508-090, Brasil.
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Thanks for that Dominik,

Giving that projection to either the locations, the raster layer generated
from the .nc file, or both, still did not work. I keep having locations
that should be on land falling far on the sea. Might this be a problem
derived from using raster with a file whose original grid distances are not
constant?

Here is a link with the original file which has the original coordinate
data.

https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0


2016-02-23 12:07 GMT-07:00 Dominik Schneider <[hidden email]
>:

> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
> dealt with this.
> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
> +a=6370000 +b=6370000 +units=m
>
> +proj=lcc which is usually what wrf is run with.
> The tricky part is:
> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
> because every WRF run is different (the WRF Preprocessing System optimizes
> the projection for the domain).
> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>
> This gives:
> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
> +a=6370000 +b=6370000 +units=m +no_defs
>
> But, wrf users like to give out lat and  long so you need to assign it:
> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>
> and then reproject the lat/long to lcc coordinates using this string:
> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
> +a=6370000 +b=6370000 +units=m +no_defs
>
> One word of caution, make sure you received the correct coordinates. Some
> variables are run cell center while some are run at cell edge. It looks
> like from your .nc file it was made by your collaborator so I assume they
> are right.
>
> That said, another word of caution, I found that the XLAT and XLONG
> variables from WRF output aren't very precise. There is a "geogrid" file
> from the preprocessing system that has the domain corners, resolution, nrow
> and ncol from which you can make a better grid using the native projection
> system (in my case it was a 4km grid). I suggest you try to get those.
>
> I hope this helps... I have to run but wanted to save people too much head
> scratching. I can get you running with more help tonight if you need.
> Dominik
>
>
> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
> wrote:
>
>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>> more giant crab this time.
>>
>> The creator of the .nc file also looked at this webpage:
>> http://www.pkrc.net/wrf-lambert.html
>> It seemed like the right proj4 string might be this one:
>>
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>     +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>
>> However this projection also does not allow me to adequately plot the
>> locations on the raster.
>>
>> Here is the .nc file. it contains several layers.
>>
>>
>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>
>>
>>
>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>
>> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>> >
>> > > On Tue, 23 Feb 2016, Alex Mandel wrote:
>> > >
>> > > > I made an attempt at it too. Investigating the original data, I'm
>> not
>> > > > sure that the projection information supplied is correct for the
>> data
>> > > > linked. When I load up the data in a unprojected space, the
>> coordinates
>> > > > don't look at all similar to any Lambert projected data I have, they
>> > > > actually look like Lat/Lon in some unprojected coordinate system,
>> > > > perhaps a different spheroid than expected.
>> > >
>> > > Does anyone have a link to the original data? Is is possible that
>> this is
>> > > the General Oblique Transformation used by modellers - that is
>> something
>> > > that feels like longlat but is recentred and oblique? Example at the
>> very
>> > > end of my GEOSTAT talk last year (slides 81-83):
>> > >
>> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf
>> > >
>> > > Roger
>> > >
>> > >
>> > For what it is worth, the General Oblique Transformation is not the only
>> > example - it's very common for modellers to have a mesh that has the
>> > "mostly-properties" of a projection, but is not actually describable
>> with
>> > standard transform + affine parameters. The main cases that I've seen
>> are
>> > polar stereographic, equal area or oblique Mercator. Often they really
>> are
>> > simple transforms and you can reconstruct without loss, but it's not
>> > usually possible to tell without exploration. It's an interesting
>> > dis-connect to see code that builds a mesh with certain properties, then
>> > only stores longitudes and latitudes - when it could be done with
>> standard
>> > tools and be stored and used much more efficiently.
>> >
>> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>> > terminology conflated in this context too. )
>> >
>> > I'm also interested to explore the original data.
>> >
>> > Cheers, Mike.
>> >
>> >
>> >
>> > > >
>> > > > -Alex
>> > > >
>> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>> > > >> Hi
>> > > >>
>> > > >> I tried to make it work but I had to give up. I wanted to reproject
>> > the
>> > > Lamberth conformal conic coordinates to long-lat but it didn't work.
>> > > >>
>> > > >> Perhaps someone can see what I did wrong. Here is what I did (data
>> in
>> > R
>> > > binary format and figure in png format both attached):
>> > > >>
>> > > >> library(raster)
>> > > >> library(maptools)
>> > > >> data(wrld_simpl)
>> > > >>
>> > > >> r <- raster("raster.grd")
>> > > >> projection(r)
>> > > >> ## > NA
>> > > >>
>> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>> ",")
>> > > >> coordinates(uro) <- ~lon+lat
>> > > >>
>> > > >> ## Set projections for the 3 data sets
>> > > >>
>> > > >> ## Lamberth's confocal conic projection with given parameters
>> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>> +lat_2=45.0
>> > > +ellps=WGS84"
>> > > >> projection(r)
>> > > >>
>> > > >> ## Assume that lon, lat are geographical coordinates (degrees
>> decimal)
>> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>> > > >>
>> > > >> ## wrld_simpl is in geographical coordinates
>> > > >> proj4string(wrld_simpl)
>> > > >>
>> > > >> ## Make figure in png format
>> > > >> ## Of course plotting data with 2 different projections will give
>> > > >> ## some distortions
>> > > >> pdf("uro.png")
>> > > >>
>> > > >> plot(r)
>> > > >> points(uro)
>> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>> 'r'
>> > > >>
>> > > >> dev.off()
>> > > >>
>> > > >> extent(r)
>> > > >> ## class       : Extent
>> > > >> ## xmin        : -131.4368
>> > > >> ## xmax        : -68.56323
>> > > >> ## ymin        : 12.35567
>> > > >> ## ymax        : 50.26619
>> > > >>
>> > > >> ## Reproject the raster to long-lat
>> > > >> ## This doesn't work (collapsed domain)
>> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>> > > +ellps=WGS84 +towgs84=0,0,0")
>> > > >>
>> > > >> ## Because
>> > > >>> extent(rp)
>> > > >> ## class       : Extent
>> > > >> ## xmin        : -100.0015
>> > > >> ## xmax        : -99.68557
>> > > >> ## ymin        : 37.70658
>> > > >> ## ymax        : 38.00046
>> > > >>
>> > > >> ## Save data in R binary format
>> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>> > > >>
>> > > >>
>> > > >> Yours sincerely / Med venlig hilsen
>> > > >>
>> > > >> Frede Aakmann Tøgersen
>> > > >> Specialist, M.Sc., Ph.D.
>> > > >> Plant Performance & Modeling
>> > > >>
>> > > >> Technology & Service Solutions
>> > > >> T +45 9730 5135
>> > > >> M +45 2547 6050
>> > > >> [hidden email]
>> > > >> http://www.vestas.com
>> > > >>
>> > > >> Company reg. name: Vestas Wind Systems A/S
>> > > >> This e-mail is subject to our e-mail disclaimer statement.
>> > > >> Please refer to www.vestas.com/legal/notice
>> > > >> If you have received this e-mail in error please contact the
>> sender.
>> > > >>
>> > > >>
>> > > >>
>> > > >> -----Original Message-----
>> > > >> From: R-sig-Geo [mailto:[hidden email]] On
>> Behalf Of
>> > > Agus Camacho
>> > > >> Sent: 22. februar 2016 19:20
>> > > >> To: [hidden email]
>> > > >> Cc: r-sig-geo
>> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
>> > > reference system implicit in a .nc file
>> > > >>
>> > > >> Thanks Alex, but the locations still fall in the sea when i plot
>> them
>> > > using
>> > > >> your recommended Solution. I looked at the sites you proposed and
>> they
>> > > have
>> > > >> other values for lat_1, lat_0, etc..
>> > > >>
>> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>> > > >>
>> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>> > > >>>> Dear all,
>> > > >>>>
>> > > >>>> Im trying to overlap these points:
>> > > >>>>
>> > > >>>
>> > >
>> >
>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>> > > >>>>
>> > > >>>> and a wrld_simpl object:
>> > > >>>> library(maptools)
>> > > >>>> data(wrld_simpl)
>> > > >>>>
>> > > >>>> Over this raster layer
>> > > >>>>
>> > > >>>
>> > >
>> >
>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>> > > >>>>
>> > > >>>> This rastr comes from a .nc file without a reference system. The
>> > > author
>> > > >>> of
>> > > >>>> that .nc file gave me the following data about the .nc.
>> > > >>>>
>> > > >>>> The projection is *Lambert conformal conic* projection
>> > > >>>> CEN_LAT = 38.0
>> > > >>>> CEN_LON = -100.0
>> > > >>>> TRUELAT1 = 25.
>> > > >>>> TRUELAT2 = 45.
>> > > >>>>
>> > > >>>> However, despite i have gone through many sites in the internet,
>> i
>> > > cant
>> > > >>>> figure it out:
>> > > >>>>
>> > > >>>> a) if that is all the data i need to set a reference system for
>> my
>> > > points
>> > > >>>> and the wrld_simp object.
>> > > >>>>
>> > > >>>> b) how to change a typical CRS object with such data
>> > > >>>>
>> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>> > > >>>>
>> > > >>>> Where do i enter the TRUELAT and CENLAT values?
>> > > >>>> Are there any site that explains easily what the fields in the
>> CRS
>> > > mean
>> > > >>> and
>> > > >>>> how to change them?
>> > > >>>>
>> > > >>>> Thanks in advance.
>> > > >>>>
>> > > >>>
>> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
>> > > >>> https://trac.osgeo.org/proj/wiki/GenParms
>> > > >>>
>> > > >>> I believe:
>> > > >>> +lat_0  = CEN_LAT   Latitude of origin
>> > > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>> > > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>> > > >>> +lon_0  = CEN_LON   Central meridian
>> > > >>>
>> > > >>> proj strings are defined by the proj4 libary. It's website listed
>> > above
>> > > >>> and the associated mailing lists or gis stackexchange would be the
>> > > >>> places to get help on it.
>> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
>> > > >>>
>> > > >>> It often helps to browse similar projections on
>> > > >>> http://spatialreference.org/
>> > > >>> http://epsg.io/
>> > > >>>
>> > > >>> Enjoy,
>> > > >>> Alex
>> > > >>>
>> > > >>>
>> > > >>
>> > > >>
>> > > >
>> > > > _______________________________________________
>> > > > R-sig-Geo mailing list
>> > > > [hidden email]
>> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >
>> > > --
>> > > Roger Bivand
>> > > Department of Economics, Norwegian School of Economics,
>> > > Helleveien 30, N-5045 Bergen, Norway.
>> > > voice: +47 55 95 93 55; fax +47 55 95 91 00
>> > > e-mail: [hidden email]
>> > > http://orcid.org/0000-0003-2392-6140
>> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>> > > http://depsy.org/person/434412
>> > > _______________________________________________
>> > > R-sig-Geo mailing list
>> > > [hidden email]
>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>> > --
>> > Dr. Michael Sumner
>> > Software and Database Engineer
>> > Australian Antarctic Division
>> > 203 Channel Highway
>> > Kingston Tasmania 7050 Australia
>> >
>> >         [[alternative HTML version deleted]]
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > [hidden email]
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>>
>>
>> --
>> Agustín Camacho Guerrero.
>> Doutor em Zoologia.
>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>> Biociências, USP.
>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>> São Paulo - SP, CEP: 05508-090, Brasil.
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>


--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Chris Reudenbach
Augustin

just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are
getting the information about the corner coordinates the subdatasets and
so on. Together with Dominiks suggestion you can do something like this:

library(gdalUtils)
library(raster)
Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES")
wrffake<- "+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m
+no_defs"
x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
                   of="GTiff",
                   ot="Byte",
                   output_Raster=TRUE,
                   verbose=TRUE)
wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619)
extent(x) <- extent(wrfExt)
projection(x) <- wrffake
plot(x)

Some remarks:
(1) I just took the first pair of  coordinates as derived from
gdalinfo("geo_em.d01.nc")
you will find 4 different coordinate pairs (i did not proof which one is
right

The data is staggered (as outlined by Dominik) So some of the corner
coordinates belongs to the staggered data and the others coordinates  to
the unstaggered ones. You will find them marked

If you have installed the netcdf libs you easily can use ncview
geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more
information of the header.

Hope this helps

cheers Chris



Am 23.02.2016 um 21:11 schrieb Agus Camacho:

> Thanks for that Dominik,
>
> Giving that projection to either the locations, the raster layer generated
> from the .nc file, or both, still did not work. I keep having locations
> that should be on land falling far on the sea. Might this be a problem
> derived from using raster with a file whose original grid distances are not
> constant?
>
> Here is a link with the original file which has the original coordinate
> data.
>
> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
>
>
> 2016-02-23 12:07 GMT-07:00 Dominik Schneider <[hidden email]
>> :
>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>> dealt with this.
>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>> +a=6370000 +b=6370000 +units=m
>>
>> +proj=lcc which is usually what wrf is run with.
>> The tricky part is:
>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> because every WRF run is different (the WRF Preprocessing System optimizes
>> the projection for the domain).
>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>
>> This gives:
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
>> +a=6370000 +b=6370000 +units=m +no_defs
>>
>> But, wrf users like to give out lat and  long so you need to assign it:
>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>> and then reproject the lat/long to lcc coordinates using this string:
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
>> +a=6370000 +b=6370000 +units=m +no_defs
>>
>> One word of caution, make sure you received the correct coordinates. Some
>> variables are run cell center while some are run at cell edge. It looks
>> like from your .nc file it was made by your collaborator so I assume they
>> are right.
>>
>> That said, another word of caution, I found that the XLAT and XLONG
>> variables from WRF output aren't very precise. There is a "geogrid" file
>> from the preprocessing system that has the domain corners, resolution, nrow
>> and ncol from which you can make a better grid using the native projection
>> system (in my case it was a 4km grid). I suggest you try to get those.
>>
>> I hope this helps... I have to run but wanted to save people too much head
>> scratching. I can get you running with more help tonight if you need.
>> Dominik
>>
>>
>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>> wrote:
>>
>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>>> more giant crab this time.
>>>
>>> The creator of the .nc file also looked at this webpage:
>>> http://www.pkrc.net/wrf-lambert.html
>>> It seemed like the right proj4 string might be this one:
>>>
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>      +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>
>>> However this projection also does not allow me to adequately plot the
>>> locations on the raster.
>>>
>>> Here is the .nc file. it contains several layers.
>>>
>>>
>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>
>>>
>>>
>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>
>>>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>>>
>>>>> On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>>>
>>>>>> I made an attempt at it too. Investigating the original data, I'm
>>> not
>>>>>> sure that the projection information supplied is correct for the
>>> data
>>>>>> linked. When I load up the data in a unprojected space, the
>>> coordinates
>>>>>> don't look at all similar to any Lambert projected data I have, they
>>>>>> actually look like Lat/Lon in some unprojected coordinate system,
>>>>>> perhaps a different spheroid than expected.
>>>>> Does anyone have a link to the original data? Is is possible that
>>> this is
>>>>> the General Oblique Transformation used by modellers - that is
>>> something
>>>>> that feels like longlat but is recentred and oblique? Example at the
>>> very
>>>>> end of my GEOSTAT talk last year (slides 81-83):
>>>>>
>>>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>>>
>>>>> Roger
>>>>>
>>>>>
>>>> For what it is worth, the General Oblique Transformation is not the only
>>>> example - it's very common for modellers to have a mesh that has the
>>>> "mostly-properties" of a projection, but is not actually describable
>>> with
>>>> standard transform + affine parameters. The main cases that I've seen
>>> are
>>>> polar stereographic, equal area or oblique Mercator. Often they really
>>> are
>>>> simple transforms and you can reconstruct without loss, but it's not
>>>> usually possible to tell without exploration. It's an interesting
>>>> dis-connect to see code that builds a mesh with certain properties, then
>>>> only stores longitudes and latitudes - when it could be done with
>>> standard
>>>> tools and be stored and used much more efficiently.
>>>>
>>>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>>> terminology conflated in this context too. )
>>>>
>>>> I'm also interested to explore the original data.
>>>>
>>>> Cheers, Mike.
>>>>
>>>>
>>>>
>>>>>> -Alex
>>>>>>
>>>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>>>>> Hi
>>>>>>>
>>>>>>> I tried to make it work but I had to give up. I wanted to reproject
>>>> the
>>>>> Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>>>>> Perhaps someone can see what I did wrong. Here is what I did (data
>>> in
>>>> R
>>>>> binary format and figure in png format both attached):
>>>>>>> library(raster)
>>>>>>> library(maptools)
>>>>>>> data(wrld_simpl)
>>>>>>>
>>>>>>> r <- raster("raster.grd")
>>>>>>> projection(r)
>>>>>>> ## > NA
>>>>>>>
>>>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>> ",")
>>>>>>> coordinates(uro) <- ~lon+lat
>>>>>>>
>>>>>>> ## Set projections for the 3 data sets
>>>>>>>
>>>>>>> ## Lamberth's confocal conic projection with given parameters
>>>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>> +lat_2=45.0
>>>>> +ellps=WGS84"
>>>>>>> projection(r)
>>>>>>>
>>>>>>> ## Assume that lon, lat are geographical coordinates (degrees
>>> decimal)
>>>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>>>>>
>>>>>>> ## wrld_simpl is in geographical coordinates
>>>>>>> proj4string(wrld_simpl)
>>>>>>>
>>>>>>> ## Make figure in png format
>>>>>>> ## Of course plotting data with 2 different projections will give
>>>>>>> ## some distortions
>>>>>>> pdf("uro.png")
>>>>>>>
>>>>>>> plot(r)
>>>>>>> points(uro)
>>>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>> 'r'
>>>>>>> dev.off()
>>>>>>>
>>>>>>> extent(r)
>>>>>>> ## class       : Extent
>>>>>>> ## xmin        : -131.4368
>>>>>>> ## xmax        : -68.56323
>>>>>>> ## ymin        : 12.35567
>>>>>>> ## ymax        : 50.26619
>>>>>>>
>>>>>>> ## Reproject the raster to long-lat
>>>>>>> ## This doesn't work (collapsed domain)
>>>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>>>> +ellps=WGS84 +towgs84=0,0,0")
>>>>>>> ## Because
>>>>>>>> extent(rp)
>>>>>>> ## class       : Extent
>>>>>>> ## xmin        : -100.0015
>>>>>>> ## xmax        : -99.68557
>>>>>>> ## ymin        : 37.70658
>>>>>>> ## ymax        : 38.00046
>>>>>>>
>>>>>>> ## Save data in R binary format
>>>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>>>>>
>>>>>>>
>>>>>>> Yours sincerely / Med venlig hilsen
>>>>>>>
>>>>>>> Frede Aakmann Tøgersen
>>>>>>> Specialist, M.Sc., Ph.D.
>>>>>>> Plant Performance & Modeling
>>>>>>>
>>>>>>> Technology & Service Solutions
>>>>>>> T +45 9730 5135
>>>>>>> M +45 2547 6050
>>>>>>> [hidden email]
>>>>>>> http://www.vestas.com
>>>>>>>
>>>>>>> Company reg. name: Vestas Wind Systems A/S
>>>>>>> This e-mail is subject to our e-mail disclaimer statement.
>>>>>>> Please refer to www.vestas.com/legal/notice
>>>>>>> If you have received this e-mail in error please contact the
>>> sender.
>>>>>>>
>>>>>>>
>>>>>>> -----Original Message-----
>>>>>>> From: R-sig-Geo [mailto:[hidden email]] On
>>> Behalf Of
>>>>> Agus Camacho
>>>>>>> Sent: 22. februar 2016 19:20
>>>>>>> To: [hidden email]
>>>>>>> Cc: r-sig-geo
>>>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
>>>>> reference system implicit in a .nc file
>>>>>>> Thanks Alex, but the locations still fall in the sea when i plot
>>> them
>>>>> using
>>>>>>> your recommended Solution. I looked at the sites you proposed and
>>> they
>>>>> have
>>>>>>> other values for lat_1, lat_0, etc..
>>>>>>>
>>>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>>>>>
>>>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>>>>>>> Dear all,
>>>>>>>>>
>>>>>>>>> Im trying to overlap these points:
>>>>>>>>>
>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>>>>>>> and a wrld_simpl object:
>>>>>>>>> library(maptools)
>>>>>>>>> data(wrld_simpl)
>>>>>>>>>
>>>>>>>>> Over this raster layer
>>>>>>>>>
>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>>>>>>> This rastr comes from a .nc file without a reference system. The
>>>>> author
>>>>>>>> of
>>>>>>>>> that .nc file gave me the following data about the .nc.
>>>>>>>>>
>>>>>>>>> The projection is *Lambert conformal conic* projection
>>>>>>>>> CEN_LAT = 38.0
>>>>>>>>> CEN_LON = -100.0
>>>>>>>>> TRUELAT1 = 25.
>>>>>>>>> TRUELAT2 = 45.
>>>>>>>>>
>>>>>>>>> However, despite i have gone through many sites in the internet,
>>> i
>>>>> cant
>>>>>>>>> figure it out:
>>>>>>>>>
>>>>>>>>> a) if that is all the data i need to set a reference system for
>>> my
>>>>> points
>>>>>>>>> and the wrld_simp object.
>>>>>>>>>
>>>>>>>>> b) how to change a typical CRS object with such data
>>>>>>>>>
>>>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>>>>>>
>>>>>>>>> Where do i enter the TRUELAT and CENLAT values?
>>>>>>>>> Are there any site that explains easily what the fields in the
>>> CRS
>>>>> mean
>>>>>>>> and
>>>>>>>>> how to change them?
>>>>>>>>>
>>>>>>>>> Thanks in advance.
>>>>>>>>>
>>>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>>>>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>>>>>>
>>>>>>>> I believe:
>>>>>>>> +lat_0  = CEN_LAT   Latitude of origin
>>>>>>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>>>>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>>>>>> +lon_0  = CEN_LON   Central meridian
>>>>>>>>
>>>>>>>> proj strings are defined by the proj4 libary. It's website listed
>>>> above
>>>>>>>> and the associated mailing lists or gis stackexchange would be the
>>>>>>>> places to get help on it.
>>>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>>>>>>
>>>>>>>> It often helps to browse similar projections on
>>>>>>>> http://spatialreference.org/
>>>>>>>> http://epsg.io/
>>>>>>>>
>>>>>>>> Enjoy,
>>>>>>>> Alex
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>> --
>>>>> Roger Bivand
>>>>> Department of Economics, Norwegian School of Economics,
>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>> e-mail: [hidden email]
>>>>> http://orcid.org/0000-0003-2392-6140
>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>> http://depsy.org/person/434412
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> --
>>>> Dr. Michael Sumner
>>>> Software and Database Engineer
>>>> Australian Antarctic Division
>>>> 203 Channel Highway
>>>> Kingston Tasmania 7050 Australia
>>>>
>>>>          [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>
>>> --
>>> Agustín Camacho Guerrero.
>>> Doutor em Zoologia.
>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>> Biociências, USP.
>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>
>>>          [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Thanks Chris, but:

> gdalinfo("geo_em.d01.nc")
NULL

In addition:

x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
                  of="GTiff",
                  ot="Byte",
                  output_Raster=TRUE,
                  verbose=TRUE)
> x
NULL

Also, when i try to install the netcdf lib:

> install.packages("netcdf")
Installing package into ‘C:/Users/Agus/Documents/R/win-library/3.2’
(as ‘lib’ is unspecified)
Warning in install.packages :
  package ‘netcdf’ is not available (for R version 3.2.3)


It seems to me that it might be necessary to know the original projection
of the locations, in order to reproject them to the projection of the .nc
file, right? or then, reproject the raster to the same projection of
wrld_simpl?

Also, when i try to reproject wrld_simpl to the proj specified by Dominik,
i have this error;

require(sp)
wrld_transf <- spTransform(wrld_simpl, CRS("+proj=lcc +lat_1=25.0
+lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere +a=6370000 +b=6370000
+units=m +no_defs"))
non finite transformation detected:
     [,1] [,2] [,3] [,4]
Error in .spTransform_Polygon(input[[i]], to_args = to_args, from_args =
from_args,  :
  failure in Polygons 145 Polygon 1 points




2016-02-23 14:22 GMT-07:00 Chris Reudenbach <[hidden email]>:

> Augustin
>
> just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are
> getting the information about the corner coordinates the subdatasets and so
> on. Together with Dominiks suggestion you can do something like this:
>
> library(gdalUtils)
> library(raster)
> Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES")
> wrffake<- "+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m
> +no_defs"
> x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
>                   of="GTiff",
>                   ot="Byte",
>                   output_Raster=TRUE,
>                   verbose=TRUE)
> wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619)
> extent(x) <- extent(wrfExt)
> projection(x) <- wrffake
> plot(x)
>
> Some remarks:
> (1) I just took the first pair of  coordinates as derived from gdalinfo("
> geo_em.d01.nc")
> you will find 4 different coordinate pairs (i did not proof which one is
> right
>
> The data is staggered (as outlined by Dominik) So some of the corner
> coordinates belongs to the staggered data and the others coordinates  to
> the unstaggered ones. You will find them marked
>
> If you have installed the netcdf libs you easily can use ncview
> geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more
> information of the header.
>
> Hope this helps
>
> cheers Chris
>
>
>
> Am 23.02.2016 um 21:11 schrieb Agus Camacho:
>
>> Thanks for that Dominik,
>>
>> Giving that projection to either the locations, the raster layer generated
>> from the .nc file, or both, still did not work. I keep having locations
>> that should be on land falling far on the sea. Might this be a problem
>> derived from using raster with a file whose original grid distances are
>> not
>> constant?
>>
>> Here is a link with the original file which has the original coordinate
>> data.
>>
>> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
>>
>>
>> 2016-02-23 12:07 GMT-07:00 Dominik Schneider <
>> [hidden email]
>>
>>> :
>>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>>>
>>> dealt with this.
>>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m
>>>
>>> +proj=lcc which is usually what wrf is run with.
>>> The tricky part is:
>>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> because every WRF run is different (the WRF Preprocessing System
>>> optimizes
>>> the projection for the domain).
>>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>>
>>> This gives:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> But, wrf users like to give out lat and  long so you need to assign it:
>>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> and then reproject the lat/long to lcc coordinates using this string:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> One word of caution, make sure you received the correct coordinates. Some
>>> variables are run cell center while some are run at cell edge. It looks
>>> like from your .nc file it was made by your collaborator so I assume they
>>> are right.
>>>
>>> That said, another word of caution, I found that the XLAT and XLONG
>>> variables from WRF output aren't very precise. There is a "geogrid" file
>>> from the preprocessing system that has the domain corners, resolution,
>>> nrow
>>> and ncol from which you can make a better grid using the native
>>> projection
>>> system (in my case it was a 4km grid). I suggest you try to get those.
>>>
>>> I hope this helps... I have to run but wanted to save people too much
>>> head
>>> scratching. I can get you running with more help tonight if you need.
>>> Dominik
>>>
>>>
>>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>>> wrote:
>>>
>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>>>> more giant crab this time.
>>>>
>>>> The creator of the .nc file also looked at this webpage:
>>>> http://www.pkrc.net/wrf-lambert.html
>>>> It seemed like the right proj4 string might be this one:
>>>>
>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>>      +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>>
>>>> However this projection also does not allow me to adequately plot the
>>>> locations on the raster.
>>>>
>>>> Here is the .nc file. it contains several layers.
>>>>
>>>>
>>>>
>>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>>
>>>>
>>>>
>>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>>
>>>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>>>>
>>>>> On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>>>>
>>>>>> I made an attempt at it too. Investigating the original data, I'm
>>>>>>>
>>>>>> not
>>>>
>>>>> sure that the projection information supplied is correct for the
>>>>>>>
>>>>>> data
>>>>
>>>>> linked. When I load up the data in a unprojected space, the
>>>>>>>
>>>>>> coordinates
>>>>
>>>>> don't look at all similar to any Lambert projected data I have, they
>>>>>>> actually look like Lat/Lon in some unprojected coordinate system,
>>>>>>> perhaps a different spheroid than expected.
>>>>>>>
>>>>>> Does anyone have a link to the original data? Is is possible that
>>>>>>
>>>>> this is
>>>>
>>>>> the General Oblique Transformation used by modellers - that is
>>>>>>
>>>>> something
>>>>
>>>>> that feels like longlat but is recentred and oblique? Example at the
>>>>>>
>>>>> very
>>>>
>>>>> end of my GEOSTAT talk last year (slides 81-83):
>>>>>>
>>>>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>>>>
>>>>>> Roger
>>>>>>
>>>>>>
>>>>>> For what it is worth, the General Oblique Transformation is not the
>>>>> only
>>>>> example - it's very common for modellers to have a mesh that has the
>>>>> "mostly-properties" of a projection, but is not actually describable
>>>>>
>>>> with
>>>>
>>>>> standard transform + affine parameters. The main cases that I've seen
>>>>>
>>>> are
>>>>
>>>>> polar stereographic, equal area or oblique Mercator. Often they really
>>>>>
>>>> are
>>>>
>>>>> simple transforms and you can reconstruct without loss, but it's not
>>>>> usually possible to tell without exploration. It's an interesting
>>>>> dis-connect to see code that builds a mesh with certain properties,
>>>>> then
>>>>> only stores longitudes and latitudes - when it could be done with
>>>>>
>>>> standard
>>>>
>>>>> tools and be stored and used much more efficiently.
>>>>>
>>>>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>>>> terminology conflated in this context too. )
>>>>>
>>>>> I'm also interested to explore the original data.
>>>>>
>>>>> Cheers, Mike.
>>>>>
>>>>>
>>>>>
>>>>> -Alex
>>>>>>>
>>>>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>>>>>
>>>>>>>> Hi
>>>>>>>>
>>>>>>>> I tried to make it work but I had to give up. I wanted to reproject
>>>>>>>>
>>>>>>> the
>>>>>
>>>>>> Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>>>>
>>>>>>> Perhaps someone can see what I did wrong. Here is what I did (data
>>>>>>>>
>>>>>>> in
>>>>
>>>>> R
>>>>>
>>>>>> binary format and figure in png format both attached):
>>>>>>
>>>>>>> library(raster)
>>>>>>>> library(maptools)
>>>>>>>> data(wrld_simpl)
>>>>>>>>
>>>>>>>> r <- raster("raster.grd")
>>>>>>>> projection(r)
>>>>>>>> ## > NA
>>>>>>>>
>>>>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>>>>>>>
>>>>>>> ",")
>>>>
>>>>> coordinates(uro) <- ~lon+lat
>>>>>>>>
>>>>>>>> ## Set projections for the 3 data sets
>>>>>>>>
>>>>>>>> ## Lamberth's confocal conic projection with given parameters
>>>>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>>>>>>>
>>>>>>> +lat_2=45.0
>>>>
>>>>> +ellps=WGS84"
>>>>>>
>>>>>>> projection(r)
>>>>>>>>
>>>>>>>> ## Assume that lon, lat are geographical coordinates (degrees
>>>>>>>>
>>>>>>> decimal)
>>>>
>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>>>>>>
>>>>>>>> ## wrld_simpl is in geographical coordinates
>>>>>>>> proj4string(wrld_simpl)
>>>>>>>>
>>>>>>>> ## Make figure in png format
>>>>>>>> ## Of course plotting data with 2 different projections will give
>>>>>>>> ## some distortions
>>>>>>>> pdf("uro.png")
>>>>>>>>
>>>>>>>> plot(r)
>>>>>>>> points(uro)
>>>>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>>>>>>>
>>>>>>> 'r'
>>>>
>>>>> dev.off()
>>>>>>>>
>>>>>>>> extent(r)
>>>>>>>> ## class       : Extent
>>>>>>>> ## xmin        : -131.4368
>>>>>>>> ## xmax        : -68.56323
>>>>>>>> ## ymin        : 12.35567
>>>>>>>> ## ymax        : 50.26619
>>>>>>>>
>>>>>>>> ## Reproject the raster to long-lat
>>>>>>>> ## This doesn't work (collapsed domain)
>>>>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>>>>>>>
>>>>>>> +ellps=WGS84 +towgs84=0,0,0")
>>>>>>
>>>>>>> ## Because
>>>>>>>>
>>>>>>>>> extent(rp)
>>>>>>>>>
>>>>>>>> ## class       : Extent
>>>>>>>> ## xmin        : -100.0015
>>>>>>>> ## xmax        : -99.68557
>>>>>>>> ## ymin        : 37.70658
>>>>>>>> ## ymax        : 38.00046
>>>>>>>>
>>>>>>>> ## Save data in R binary format
>>>>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>>>>>>
>>>>>>>>
>>>>>>>> Yours sincerely / Med venlig hilsen
>>>>>>>>
>>>>>>>> Frede Aakmann Tøgersen
>>>>>>>> Specialist, M.Sc., Ph.D.
>>>>>>>> Plant Performance & Modeling
>>>>>>>>
>>>>>>>> Technology & Service Solutions
>>>>>>>> T +45 9730 5135
>>>>>>>> M +45 2547 6050
>>>>>>>> [hidden email]
>>>>>>>> http://www.vestas.com
>>>>>>>>
>>>>>>>> Company reg. name: Vestas Wind Systems A/S
>>>>>>>> This e-mail is subject to our e-mail disclaimer statement.
>>>>>>>> Please refer to www.vestas.com/legal/notice
>>>>>>>> If you have received this e-mail in error please contact the
>>>>>>>>
>>>>>>> sender.
>>>>
>>>>>
>>>>>>>>
>>>>>>>> -----Original Message-----
>>>>>>>> From: R-sig-Geo [mailto:[hidden email]] On
>>>>>>>>
>>>>>>> Behalf Of
>>>>
>>>>> Agus Camacho
>>>>>>
>>>>>>> Sent: 22. februar 2016 19:20
>>>>>>>> To: [hidden email]
>>>>>>>> Cc: r-sig-geo
>>>>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
>>>>>>>>
>>>>>>> reference system implicit in a .nc file
>>>>>>
>>>>>>> Thanks Alex, but the locations still fall in the sea when i plot
>>>>>>>>
>>>>>>> them
>>>>
>>>>> using
>>>>>>
>>>>>>> your recommended Solution. I looked at the sites you proposed and
>>>>>>>>
>>>>>>> they
>>>>
>>>>> have
>>>>>>
>>>>>>> other values for lat_1, lat_0, etc..
>>>>>>>>
>>>>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>>>>>>
>>>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>>>>>>>
>>>>>>>>>> Dear all,
>>>>>>>>>>
>>>>>>>>>> Im trying to overlap these points:
>>>>>>>>>>
>>>>>>>>>>
>>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>>
>>>>> and a wrld_simpl object:
>>>>>>>>>> library(maptools)
>>>>>>>>>> data(wrld_simpl)
>>>>>>>>>>
>>>>>>>>>> Over this raster layer
>>>>>>>>>>
>>>>>>>>>>
>>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>>
>>>>> This rastr comes from a .nc file without a reference system. The
>>>>>>>>>>
>>>>>>>>> author
>>>>>>
>>>>>>> of
>>>>>>>>>
>>>>>>>>>> that .nc file gave me the following data about the .nc.
>>>>>>>>>>
>>>>>>>>>> The projection is *Lambert conformal conic* projection
>>>>>>>>>> CEN_LAT = 38.0
>>>>>>>>>> CEN_LON = -100.0
>>>>>>>>>> TRUELAT1 = 25.
>>>>>>>>>> TRUELAT2 = 45.
>>>>>>>>>>
>>>>>>>>>> However, despite i have gone through many sites in the internet,
>>>>>>>>>>
>>>>>>>>> i
>>>>
>>>>> cant
>>>>>>
>>>>>>> figure it out:
>>>>>>>>>>
>>>>>>>>>> a) if that is all the data i need to set a reference system for
>>>>>>>>>>
>>>>>>>>> my
>>>>
>>>>> points
>>>>>>
>>>>>>> and the wrld_simp object.
>>>>>>>>>>
>>>>>>>>>> b) how to change a typical CRS object with such data
>>>>>>>>>>
>>>>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>>>>>>>
>>>>>>>>>> Where do i enter the TRUELAT and CENLAT values?
>>>>>>>>>> Are there any site that explains easily what the fields in the
>>>>>>>>>>
>>>>>>>>> CRS
>>>>
>>>>> mean
>>>>>>
>>>>>>> and
>>>>>>>>>
>>>>>>>>>> how to change them?
>>>>>>>>>>
>>>>>>>>>> Thanks in advance.
>>>>>>>>>>
>>>>>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>>>>>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>>>>>>>
>>>>>>>>> I believe:
>>>>>>>>> +lat_0  = CEN_LAT   Latitude of origin
>>>>>>>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>>>>>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>>>>>>> +lon_0  = CEN_LON   Central meridian
>>>>>>>>>
>>>>>>>>> proj strings are defined by the proj4 libary. It's website listed
>>>>>>>>>
>>>>>>>> above
>>>>>
>>>>>> and the associated mailing lists or gis stackexchange would be the
>>>>>>>>> places to get help on it.
>>>>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>>>>>>>
>>>>>>>>> It often helps to browse similar projections on
>>>>>>>>> http://spatialreference.org/
>>>>>>>>> http://epsg.io/
>>>>>>>>>
>>>>>>>>> Enjoy,
>>>>>>>>> Alex
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>> R-sig-Geo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>
>>>>>> --
>>>>>> Roger Bivand
>>>>>> Department of Economics, Norwegian School of Economics,
>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>>> e-mail: [hidden email]
>>>>>> http://orcid.org/0000-0003-2392-6140
>>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>>> http://depsy.org/person/434412
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>
>>>>> --
>>>>> Dr. Michael Sumner
>>>>> Software and Database Engineer
>>>>> Australian Antarctic Division
>>>>> 203 Channel Highway
>>>>> Kingston Tasmania 7050 Australia
>>>>>
>>>>>          [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>>
>>>>
>>>> --
>>>> Agustín Camacho Guerrero.
>>>> Doutor em Zoologia.
>>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>>> Biociências, USP.
>>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>>
>>>>          [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>>
>>>
>>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Chris Reudenbach
In reply to this post by Chris Reudenbach
Agus

sorry I did miss the crucial part:

if you are doing as suggestest you have to define manually the
Lambertian ymin xmax ymin ymax using the header information of the nc file.

here an example for the unstaggered  data

library(gdalUtils)
library(raster)
library(proj4)
library(mapview)

r<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
                   of="GTiff",
                   ot="Byte",
                   output_Raster=TRUE,
                   verbose=TRUE)

finfo <- gdalinfo("geo_em.d01.nc")
##extract parameters
lat_1=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT1",
finfo))],"=")[[1]][2])
lat_2=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT2",
finfo))],"=")[[1]][2])
lat0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LAT",
finfo))],"=")[[1]][2])
lon0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LON",
finfo))],"=")[[1]][2])
dx=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DX",
finfo))],"=")[[1]][2])
dy=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DY",
finfo))],"=")[[1]][2])
y=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG",
finfo))],"=")[[1]][2])
x=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG",
finfo))],"=")[[1]][2])
x_0=0
y_0=0

## generate compliant prj4 string
projLcc=paste("+proj=lcc +lat_1=",lat_1," +lat_2=",lat_2,
            " +lat_0=",lat0," +lon_0=",lon0, " +x_0=",x_0,"
+y_0=",y_0,sep="")

## project centre coordinates
tr <- ptransform(cbind(lon0, lat0)/180*pi,'+proj=longlat +datum=WGS84
+no_defs',proj)

## calculate extent using the Lambertian units (m)
xmin=as.integer(tr[1,1]-(x/2*dx-dx))
xmax=as.integer(tr[1,1]+(x/2*dx-dx))
ymin=as.integer(tr[1,2]-(y/2*dy-dy))
ymax=as.integer(tr[1,2]+(y/2*dy-dy))
wrfLccExt<-extent(xmin,xmax,ymin,ymax)
extent(r) <- extent(wrfLccExt)
projection(r) <- projLcc
mapview(r)


cheers Chris


Am 23.02.2016 um 22:22 schrieb Chris Reudenbach:

> Augustin
>
> just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are
> getting the information about the corner coordinates the subdatasets
> and so on. Together with Dominiks suggestion you can do something like
> this:
>
> library(gdalUtils)
> library(raster)
> Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES")
> wrffake<- "+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m
> +no_defs"
> x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
>                   of="GTiff",
>                   ot="Byte",
>                   output_Raster=TRUE,
>                   verbose=TRUE)
> wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619)
> extent(x) <- extent(wrfExt)
> projection(x) <- wrffake
> plot(x)
>
> Some remarks:
> (1) I just took the first pair of  coordinates as derived from
> gdalinfo("geo_em.d01.nc")
> you will find 4 different coordinate pairs (i did not proof which one
> is right
>
> The data is staggered (as outlined by Dominik) So some of the corner
> coordinates belongs to the staggered data and the others coordinates  
> to the unstaggered ones. You will find them marked
>
> If you have installed the netcdf libs you easily can use ncview
> geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more
> information of the header.
>
> Hope this helps
>
> cheers Chris
>
>
>
> Am 23.02.2016 um 21:11 schrieb Agus Camacho:
>> Thanks for that Dominik,
>>
>> Giving that projection to either the locations, the raster layer
>> generated
>> from the .nc file, or both, still did not work. I keep having locations
>> that should be on land falling far on the sea. Might this be a problem
>> derived from using raster with a file whose original grid distances
>> are not
>> constant?
>>
>> Here is a link with the original file which has the original coordinate
>> data.
>>
>> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
>>
>>
>> 2016-02-23 12:07 GMT-07:00 Dominik Schneider
>> <[hidden email]
>>> :
>>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>>> dealt with this.
>>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m
>>>
>>> +proj=lcc which is usually what wrf is run with.
>>> The tricky part is:
>>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> because every WRF run is different (the WRF Preprocessing System
>>> optimizes
>>> the projection for the domain).
>>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>>
>>> This gives:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> But, wrf users like to give out lat and  long so you need to assign it:
>>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> and then reproject the lat/long to lcc coordinates using this string:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> One word of caution, make sure you received the correct coordinates.
>>> Some
>>> variables are run cell center while some are run at cell edge. It looks
>>> like from your .nc file it was made by your collaborator so I assume
>>> they
>>> are right.
>>>
>>> That said, another word of caution, I found that the XLAT and XLONG
>>> variables from WRF output aren't very precise. There is a "geogrid"
>>> file
>>> from the preprocessing system that has the domain corners,
>>> resolution, nrow
>>> and ncol from which you can make a better grid using the native
>>> projection
>>> system (in my case it was a 4km grid). I suggest you try to get those.
>>>
>>> I hope this helps... I have to run but wanted to save people too
>>> much head
>>> scratching. I can get you running with more help tonight if you need.
>>> Dominik
>>>
>>>
>>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>>> wrote:
>>>
>>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill
>>>> bring
>>>> more giant crab this time.
>>>>
>>>> The creator of the .nc file also looked at this webpage:
>>>> http://www.pkrc.net/wrf-lambert.html
>>>> It seemed like the right proj4 string might be this one:
>>>>
>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>>      +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>>
>>>> However this projection also does not allow me to adequately plot the
>>>> locations on the raster.
>>>>
>>>> Here is the .nc file. it contains several layers.
>>>>
>>>>
>>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0 
>>>>
>>>>
>>>>
>>>>
>>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>>
>>>>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]>
>>>>> wrote:
>>>>>
>>>>>> On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>>>>
>>>>>>> I made an attempt at it too. Investigating the original data, I'm
>>>> not
>>>>>>> sure that the projection information supplied is correct for the
>>>> data
>>>>>>> linked. When I load up the data in a unprojected space, the
>>>> coordinates
>>>>>>> don't look at all similar to any Lambert projected data I have,
>>>>>>> they
>>>>>>> actually look like Lat/Lon in some unprojected coordinate system,
>>>>>>> perhaps a different spheroid than expected.
>>>>>> Does anyone have a link to the original data? Is is possible that
>>>> this is
>>>>>> the General Oblique Transformation used by modellers - that is
>>>> something
>>>>>> that feels like longlat but is recentred and oblique? Example at the
>>>> very
>>>>>> end of my GEOSTAT talk last year (slides 81-83):
>>>>>>
>>>>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>>>>
>>>>>> Roger
>>>>>>
>>>>>>
>>>>> For what it is worth, the General Oblique Transformation is not
>>>>> the only
>>>>> example - it's very common for modellers to have a mesh that has the
>>>>> "mostly-properties" of a projection, but is not actually describable
>>>> with
>>>>> standard transform + affine parameters. The main cases that I've seen
>>>> are
>>>>> polar stereographic, equal area or oblique Mercator. Often they
>>>>> really
>>>> are
>>>>> simple transforms and you can reconstruct without loss, but it's not
>>>>> usually possible to tell without exploration. It's an interesting
>>>>> dis-connect to see code that builds a mesh with certain
>>>>> properties, then
>>>>> only stores longitudes and latitudes - when it could be done with
>>>> standard
>>>>> tools and be stored and used much more efficiently.
>>>>>
>>>>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>>>> terminology conflated in this context too. )
>>>>>
>>>>> I'm also interested to explore the original data.
>>>>>
>>>>> Cheers, Mike.
>>>>>
>>>>>
>>>>>
>>>>>>> -Alex
>>>>>>>
>>>>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>>>>>> Hi
>>>>>>>>
>>>>>>>> I tried to make it work but I had to give up. I wanted to
>>>>>>>> reproject
>>>>> the
>>>>>> Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>>>>>> Perhaps someone can see what I did wrong. Here is what I did (data
>>>> in
>>>>> R
>>>>>> binary format and figure in png format both attached):
>>>>>>>> library(raster)
>>>>>>>> library(maptools)
>>>>>>>> data(wrld_simpl)
>>>>>>>>
>>>>>>>> r <- raster("raster.grd")
>>>>>>>> projection(r)
>>>>>>>> ## > NA
>>>>>>>>
>>>>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>>> ",")
>>>>>>>> coordinates(uro) <- ~lon+lat
>>>>>>>>
>>>>>>>> ## Set projections for the 3 data sets
>>>>>>>>
>>>>>>>> ## Lamberth's confocal conic projection with given parameters
>>>>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>>> +lat_2=45.0
>>>>>> +ellps=WGS84"
>>>>>>>> projection(r)
>>>>>>>>
>>>>>>>> ## Assume that lon, lat are geographical coordinates (degrees
>>>> decimal)
>>>>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>>>>>>
>>>>>>>> ## wrld_simpl is in geographical coordinates
>>>>>>>> proj4string(wrld_simpl)
>>>>>>>>
>>>>>>>> ## Make figure in png format
>>>>>>>> ## Of course plotting data with 2 different projections will give
>>>>>>>> ## some distortions
>>>>>>>> pdf("uro.png")
>>>>>>>>
>>>>>>>> plot(r)
>>>>>>>> points(uro)
>>>>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>>> 'r'
>>>>>>>> dev.off()
>>>>>>>>
>>>>>>>> extent(r)
>>>>>>>> ## class       : Extent
>>>>>>>> ## xmin        : -131.4368
>>>>>>>> ## xmax        : -68.56323
>>>>>>>> ## ymin        : 12.35567
>>>>>>>> ## ymax        : 50.26619
>>>>>>>>
>>>>>>>> ## Reproject the raster to long-lat
>>>>>>>> ## This doesn't work (collapsed domain)
>>>>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>>>>> +ellps=WGS84 +towgs84=0,0,0")
>>>>>>>> ## Because
>>>>>>>>> extent(rp)
>>>>>>>> ## class       : Extent
>>>>>>>> ## xmin        : -100.0015
>>>>>>>> ## xmax        : -99.68557
>>>>>>>> ## ymin        : 37.70658
>>>>>>>> ## ymax        : 38.00046
>>>>>>>>
>>>>>>>> ## Save data in R binary format
>>>>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>>>>>>
>>>>>>>>
>>>>>>>> Yours sincerely / Med venlig hilsen
>>>>>>>>
>>>>>>>> Frede Aakmann Tøgersen
>>>>>>>> Specialist, M.Sc., Ph.D.
>>>>>>>> Plant Performance & Modeling
>>>>>>>>
>>>>>>>> Technology & Service Solutions
>>>>>>>> T +45 9730 5135
>>>>>>>> M +45 2547 6050
>>>>>>>> [hidden email]
>>>>>>>> http://www.vestas.com
>>>>>>>>
>>>>>>>> Company reg. name: Vestas Wind Systems A/S
>>>>>>>> This e-mail is subject to our e-mail disclaimer statement.
>>>>>>>> Please refer to www.vestas.com/legal/notice
>>>>>>>> If you have received this e-mail in error please contact the
>>>> sender.
>>>>>>>>
>>>>>>>>
>>>>>>>> -----Original Message-----
>>>>>>>> From: R-sig-Geo [mailto:[hidden email]] On
>>>> Behalf Of
>>>>>> Agus Camacho
>>>>>>>> Sent: 22. februar 2016 19:20
>>>>>>>> To: [hidden email]
>>>>>>>> Cc: r-sig-geo
>>>>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl
>>>>>>>> to a
>>>>>> reference system implicit in a .nc file
>>>>>>>> Thanks Alex, but the locations still fall in the sea when i plot
>>>> them
>>>>>> using
>>>>>>>> your recommended Solution. I looked at the sites you proposed and
>>>> they
>>>>>> have
>>>>>>>> other values for lat_1, lat_0, etc..
>>>>>>>>
>>>>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>>>>>>
>>>>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>>>>>>>> Dear all,
>>>>>>>>>>
>>>>>>>>>> Im trying to overlap these points:
>>>>>>>>>>
>>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0 
>>>>
>>>>>>>>>> and a wrld_simpl object:
>>>>>>>>>> library(maptools)
>>>>>>>>>> data(wrld_simpl)
>>>>>>>>>>
>>>>>>>>>> Over this raster layer
>>>>>>>>>>
>>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0 
>>>>
>>>>>>>>>> This rastr comes from a .nc file without a reference system. The
>>>>>> author
>>>>>>>>> of
>>>>>>>>>> that .nc file gave me the following data about the .nc.
>>>>>>>>>>
>>>>>>>>>> The projection is *Lambert conformal conic* projection
>>>>>>>>>> CEN_LAT = 38.0
>>>>>>>>>> CEN_LON = -100.0
>>>>>>>>>> TRUELAT1 = 25.
>>>>>>>>>> TRUELAT2 = 45.
>>>>>>>>>>
>>>>>>>>>> However, despite i have gone through many sites in the internet,
>>>> i
>>>>>> cant
>>>>>>>>>> figure it out:
>>>>>>>>>>
>>>>>>>>>> a) if that is all the data i need to set a reference system for
>>>> my
>>>>>> points
>>>>>>>>>> and the wrld_simp object.
>>>>>>>>>>
>>>>>>>>>> b) how to change a typical CRS object with such data
>>>>>>>>>>
>>>>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>>>>>>>
>>>>>>>>>> Where do i enter the TRUELAT and CENLAT values?
>>>>>>>>>> Are there any site that explains easily what the fields in the
>>>> CRS
>>>>>> mean
>>>>>>>>> and
>>>>>>>>>> how to change them?
>>>>>>>>>>
>>>>>>>>>> Thanks in advance.
>>>>>>>>>>
>>>>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>>>>>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>>>>>>>
>>>>>>>>> I believe:
>>>>>>>>> +lat_0  = CEN_LAT   Latitude of origin
>>>>>>>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>>>>>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>>>>>>> +lon_0  = CEN_LON   Central meridian
>>>>>>>>>
>>>>>>>>> proj strings are defined by the proj4 libary. It's website listed
>>>>> above
>>>>>>>>> and the associated mailing lists or gis stackexchange would be
>>>>>>>>> the
>>>>>>>>> places to get help on it.
>>>>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>>>>>>>
>>>>>>>>> It often helps to browse similar projections on
>>>>>>>>> http://spatialreference.org/
>>>>>>>>> http://epsg.io/
>>>>>>>>>
>>>>>>>>> Enjoy,
>>>>>>>>> Alex
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-sig-Geo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>> --
>>>>>> Roger Bivand
>>>>>> Department of Economics, Norwegian School of Economics,
>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>>> e-mail: [hidden email]
>>>>>> http://orcid.org/0000-0003-2392-6140
>>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>>> http://depsy.org/person/434412
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>> --
>>>>> Dr. Michael Sumner
>>>>> Software and Database Engineer
>>>>> Australian Antarctic Division
>>>>> 203 Channel Highway
>>>>> Kingston Tasmania 7050 Australia
>>>>>
>>>>>          [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>>
>>>> --
>>>> Agustín Camacho Guerrero.
>>>> Doutor em Zoologia.
>>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>>> Biociências, USP.
>>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>>
>>>>          [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Thanks again, but with gdalinfo  and gdal_translate giving NULL after
running to the end, i just cant do that..

2016-02-23 15:36 GMT-07:00 Chris Reudenbach <[hidden email]>:

> Agus
>
> sorry I did miss the crucial part:
>
> if you are doing as suggestest you have to define manually the Lambertian
> ymin xmax ymin ymax using the header information of the nc file.
>
> here an example for the unstaggered  data
>
> library(gdalUtils)
> library(raster)
> library(proj4)
> library(mapview)
>
> r<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
>                   of="GTiff",
>                   ot="Byte",
>                   output_Raster=TRUE,
>                   verbose=TRUE)
>
> finfo <- gdalinfo("geo_em.d01.nc")
> ##extract parameters
> lat_1=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT1",
> finfo))],"=")[[1]][2])
> lat_2=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#TRUELAT2",
> finfo))],"=")[[1]][2])
> lat0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LAT",
> finfo))],"=")[[1]][2])
> lon0=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#CEN_LON",
> finfo))],"=")[[1]][2])
> dx=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DX",
> finfo))],"=")[[1]][2])
> dy=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#DY",
> finfo))],"=")[[1]][2])
> y=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#SOUTH-NORTH_PATCH_END_UNSTAG",
> finfo))],"=")[[1]][2])
> x=as.numeric(strsplit(finfo[which(grepl("NC_GLOBAL#WEST-EAST_PATCH_END_UNSTAG",
> finfo))],"=")[[1]][2])
> x_0=0
> y_0=0
>
> ## generate compliant prj4 string
> projLcc=paste("+proj=lcc +lat_1=",lat_1," +lat_2=",lat_2,
>            " +lat_0=",lat0," +lon_0=",lon0, " +x_0=",x_0,"
> +y_0=",y_0,sep="")
>
> ## project centre coordinates
> tr <- ptransform(cbind(lon0, lat0)/180*pi,'+proj=longlat +datum=WGS84
> +no_defs',proj)
>
> ## calculate extent using the Lambertian units (m)
> xmin=as.integer(tr[1,1]-(x/2*dx-dx))
> xmax=as.integer(tr[1,1]+(x/2*dx-dx))
> ymin=as.integer(tr[1,2]-(y/2*dy-dy))
> ymax=as.integer(tr[1,2]+(y/2*dy-dy))
> wrfLccExt<-extent(xmin,xmax,ymin,ymax)
> extent(r) <- extent(wrfLccExt)
> projection(r) <- projLcc
> mapview(r)
>
>
> cheers Chris
>
>
>
> Am 23.02.2016 um 22:22 schrieb Chris Reudenbach:
>
>> Augustin
>>
>> just quick and dirty if you run gdalinfo("geo_em.d01.nc") your are
>> getting the information about the corner coordinates the subdatasets and so
>> on. Together with Dominiks suggestion you can do something like this:
>>
>> library(gdalUtils)
>> library(raster)
>> Sys.setenv(GDAL_NETCDF_BOTTOMUP="YES")
>> wrffake<- "+proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m
>> +no_defs"
>> x<-gdal_translate('NETCDF:"geo_em.d01.nc":LANDMASK', 'landmask.tif',
>>                   of="GTiff",
>>                   ot="Byte",
>>                   output_Raster=TRUE,
>>                   verbose=TRUE)
>> wrfExt<-extent(-151.29639,-48.703613,12.355667,50.26619)
>> extent(x) <- extent(wrfExt)
>> projection(x) <- wrffake
>> plot(x)
>>
>> Some remarks:
>> (1) I just took the first pair of  coordinates as derived from gdalinfo("
>> geo_em.d01.nc")
>> you will find 4 different coordinate pairs (i did not proof which one is
>> right
>>
>> The data is staggered (as outlined by Dominik) So some of the corner
>> coordinates belongs to the staggered data and the others coordinates  to
>> the unstaggered ones. You will find them marked
>>
>> If you have installed the netcdf libs you easily can use ncview
>> geo_em.d01.nc or ncdump -h geo_em.d01.nc to view the data or get more
>> information of the header.
>>
>> Hope this helps
>>
>> cheers Chris
>>
>>
>>
>> Am 23.02.2016 um 21:11 schrieb Agus Camacho:
>>
>>> Thanks for that Dominik,
>>>
>>> Giving that projection to either the locations, the raster layer
>>> generated
>>> from the .nc file, or both, still did not work. I keep having locations
>>> that should be on land falling far on the sea. Might this be a problem
>>> derived from using raster with a file whose original grid distances are
>>> not
>>> constant?
>>>
>>> Here is a link with the original file which has the original coordinate
>>> data.
>>>
>>> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
>>>
>>>
>>> 2016-02-23 12:07 GMT-07:00 Dominik Schneider <
>>> [hidden email]
>>>
>>>> :
>>>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>>>> dealt with this.
>>>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>>>> +a=6370000 +b=6370000 +units=m
>>>>
>>>> +proj=lcc which is usually what wrf is run with.
>>>> The tricky part is:
>>>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>>> because every WRF run is different (the WRF Preprocessing System
>>>> optimizes
>>>> the projection for the domain).
>>>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>>>
>>>> This gives:
>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>>> +ellps=sphere
>>>> +a=6370000 +b=6370000 +units=m +no_defs
>>>>
>>>> But, wrf users like to give out lat and  long so you need to assign it:
>>>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>>
>>>> and then reproject the lat/long to lcc coordinates using this string:
>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>>> +ellps=sphere
>>>> +a=6370000 +b=6370000 +units=m +no_defs
>>>>
>>>> One word of caution, make sure you received the correct coordinates.
>>>> Some
>>>> variables are run cell center while some are run at cell edge. It looks
>>>> like from your .nc file it was made by your collaborator so I assume
>>>> they
>>>> are right.
>>>>
>>>> That said, another word of caution, I found that the XLAT and XLONG
>>>> variables from WRF output aren't very precise. There is a "geogrid" file
>>>> from the preprocessing system that has the domain corners, resolution,
>>>> nrow
>>>> and ncol from which you can make a better grid using the native
>>>> projection
>>>> system (in my case it was a 4km grid). I suggest you try to get those.
>>>>
>>>> I hope this helps... I have to run but wanted to save people too much
>>>> head
>>>> scratching. I can get you running with more help tonight if you need.
>>>> Dominik
>>>>
>>>>
>>>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>>>> wrote:
>>>>
>>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill
>>>>> bring
>>>>> more giant crab this time.
>>>>>
>>>>> The creator of the .nc file also looked at this webpage:
>>>>> http://www.pkrc.net/wrf-lambert.html
>>>>> It seemed like the right proj4 string might be this one:
>>>>>
>>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>>>      +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>>>
>>>>> However this projection also does not allow me to adequately plot the
>>>>> locations on the raster.
>>>>>
>>>>> Here is the .nc file. it contains several layers.
>>>>>
>>>>>
>>>>>
>>>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>>>
>>>>>
>>>>>
>>>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>>>
>>>>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>>>>>
>>>>>> On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>>>>>
>>>>>>> I made an attempt at it too. Investigating the original data, I'm
>>>>>>>>
>>>>>>> not
>>>>>
>>>>>> sure that the projection information supplied is correct for the
>>>>>>>>
>>>>>>> data
>>>>>
>>>>>> linked. When I load up the data in a unprojected space, the
>>>>>>>>
>>>>>>> coordinates
>>>>>
>>>>>> don't look at all similar to any Lambert projected data I have, they
>>>>>>>> actually look like Lat/Lon in some unprojected coordinate system,
>>>>>>>> perhaps a different spheroid than expected.
>>>>>>>>
>>>>>>> Does anyone have a link to the original data? Is is possible that
>>>>>>>
>>>>>> this is
>>>>>
>>>>>> the General Oblique Transformation used by modellers - that is
>>>>>>>
>>>>>> something
>>>>>
>>>>>> that feels like longlat but is recentred and oblique? Example at the
>>>>>>>
>>>>>> very
>>>>>
>>>>>> end of my GEOSTAT talk last year (slides 81-83):
>>>>>>>
>>>>>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>>>>>
>>>>>>> Roger
>>>>>>>
>>>>>>>
>>>>>>> For what it is worth, the General Oblique Transformation is not the
>>>>>> only
>>>>>> example - it's very common for modellers to have a mesh that has the
>>>>>> "mostly-properties" of a projection, but is not actually describable
>>>>>>
>>>>> with
>>>>>
>>>>>> standard transform + affine parameters. The main cases that I've seen
>>>>>>
>>>>> are
>>>>>
>>>>>> polar stereographic, equal area or oblique Mercator. Often they really
>>>>>>
>>>>> are
>>>>>
>>>>>> simple transforms and you can reconstruct without loss, but it's not
>>>>>> usually possible to tell without exploration. It's an interesting
>>>>>> dis-connect to see code that builds a mesh with certain properties,
>>>>>> then
>>>>>> only stores longitudes and latitudes - when it could be done with
>>>>>>
>>>>> standard
>>>>>
>>>>>> tools and be stored and used much more efficiently.
>>>>>>
>>>>>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>>>>> terminology conflated in this context too. )
>>>>>>
>>>>>> I'm also interested to explore the original data.
>>>>>>
>>>>>> Cheers, Mike.
>>>>>>
>>>>>>
>>>>>>
>>>>>> -Alex
>>>>>>>>
>>>>>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>>>>>>
>>>>>>>>> Hi
>>>>>>>>>
>>>>>>>>> I tried to make it work but I had to give up. I wanted to reproject
>>>>>>>>>
>>>>>>>> the
>>>>>>
>>>>>>> Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>>>>>
>>>>>>>> Perhaps someone can see what I did wrong. Here is what I did (data
>>>>>>>>>
>>>>>>>> in
>>>>>
>>>>>> R
>>>>>>
>>>>>>> binary format and figure in png format both attached):
>>>>>>>
>>>>>>>> library(raster)
>>>>>>>>> library(maptools)
>>>>>>>>> data(wrld_simpl)
>>>>>>>>>
>>>>>>>>> r <- raster("raster.grd")
>>>>>>>>> projection(r)
>>>>>>>>> ## > NA
>>>>>>>>>
>>>>>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>>>>>>>>
>>>>>>>> ",")
>>>>>
>>>>>> coordinates(uro) <- ~lon+lat
>>>>>>>>>
>>>>>>>>> ## Set projections for the 3 data sets
>>>>>>>>>
>>>>>>>>> ## Lamberth's confocal conic projection with given parameters
>>>>>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>>>>>>>>
>>>>>>>> +lat_2=45.0
>>>>>
>>>>>> +ellps=WGS84"
>>>>>>>
>>>>>>>> projection(r)
>>>>>>>>>
>>>>>>>>> ## Assume that lon, lat are geographical coordinates (degrees
>>>>>>>>>
>>>>>>>> decimal)
>>>>>
>>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>>>>>>>
>>>>>>>>> ## wrld_simpl is in geographical coordinates
>>>>>>>>> proj4string(wrld_simpl)
>>>>>>>>>
>>>>>>>>> ## Make figure in png format
>>>>>>>>> ## Of course plotting data with 2 different projections will give
>>>>>>>>> ## some distortions
>>>>>>>>> pdf("uro.png")
>>>>>>>>>
>>>>>>>>> plot(r)
>>>>>>>>> points(uro)
>>>>>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>>>>>>>>
>>>>>>>> 'r'
>>>>>
>>>>>> dev.off()
>>>>>>>>>
>>>>>>>>> extent(r)
>>>>>>>>> ## class       : Extent
>>>>>>>>> ## xmin        : -131.4368
>>>>>>>>> ## xmax        : -68.56323
>>>>>>>>> ## ymin        : 12.35567
>>>>>>>>> ## ymax        : 50.26619
>>>>>>>>>
>>>>>>>>> ## Reproject the raster to long-lat
>>>>>>>>> ## This doesn't work (collapsed domain)
>>>>>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>>>>>>>>
>>>>>>>> +ellps=WGS84 +towgs84=0,0,0")
>>>>>>>
>>>>>>>> ## Because
>>>>>>>>>
>>>>>>>>>> extent(rp)
>>>>>>>>>>
>>>>>>>>> ## class       : Extent
>>>>>>>>> ## xmin        : -100.0015
>>>>>>>>> ## xmax        : -99.68557
>>>>>>>>> ## ymin        : 37.70658
>>>>>>>>> ## ymax        : 38.00046
>>>>>>>>>
>>>>>>>>> ## Save data in R binary format
>>>>>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> Yours sincerely / Med venlig hilsen
>>>>>>>>>
>>>>>>>>> Frede Aakmann Tøgersen
>>>>>>>>> Specialist, M.Sc., Ph.D.
>>>>>>>>> Plant Performance & Modeling
>>>>>>>>>
>>>>>>>>> Technology & Service Solutions
>>>>>>>>> T +45 9730 5135
>>>>>>>>> M +45 2547 6050
>>>>>>>>> [hidden email]
>>>>>>>>> http://www.vestas.com
>>>>>>>>>
>>>>>>>>> Company reg. name: Vestas Wind Systems A/S
>>>>>>>>> This e-mail is subject to our e-mail disclaimer statement.
>>>>>>>>> Please refer to www.vestas.com/legal/notice
>>>>>>>>> If you have received this e-mail in error please contact the
>>>>>>>>>
>>>>>>>> sender.
>>>>>
>>>>>>
>>>>>>>>>
>>>>>>>>> -----Original Message-----
>>>>>>>>> From: R-sig-Geo [mailto:[hidden email]] On
>>>>>>>>>
>>>>>>>> Behalf Of
>>>>>
>>>>>> Agus Camacho
>>>>>>>
>>>>>>>> Sent: 22. februar 2016 19:20
>>>>>>>>> To: [hidden email]
>>>>>>>>> Cc: r-sig-geo
>>>>>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
>>>>>>>>>
>>>>>>>> reference system implicit in a .nc file
>>>>>>>
>>>>>>>> Thanks Alex, but the locations still fall in the sea when i plot
>>>>>>>>>
>>>>>>>> them
>>>>>
>>>>>> using
>>>>>>>
>>>>>>>> your recommended Solution. I looked at the sites you proposed and
>>>>>>>>>
>>>>>>>> they
>>>>>
>>>>>> have
>>>>>>>
>>>>>>>> other values for lat_1, lat_0, etc..
>>>>>>>>>
>>>>>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>>>>>>>
>>>>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>>>>>>>>
>>>>>>>>>>> Dear all,
>>>>>>>>>>>
>>>>>>>>>>> Im trying to overlap these points:
>>>>>>>>>>>
>>>>>>>>>>>
>>>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>>>
>>>>>> and a wrld_simpl object:
>>>>>>>>>>> library(maptools)
>>>>>>>>>>> data(wrld_simpl)
>>>>>>>>>>>
>>>>>>>>>>> Over this raster layer
>>>>>>>>>>>
>>>>>>>>>>>
>>>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>>>
>>>>>> This rastr comes from a .nc file without a reference system. The
>>>>>>>>>>>
>>>>>>>>>> author
>>>>>>>
>>>>>>>> of
>>>>>>>>>>
>>>>>>>>>>> that .nc file gave me the following data about the .nc.
>>>>>>>>>>>
>>>>>>>>>>> The projection is *Lambert conformal conic* projection
>>>>>>>>>>> CEN_LAT = 38.0
>>>>>>>>>>> CEN_LON = -100.0
>>>>>>>>>>> TRUELAT1 = 25.
>>>>>>>>>>> TRUELAT2 = 45.
>>>>>>>>>>>
>>>>>>>>>>> However, despite i have gone through many sites in the internet,
>>>>>>>>>>>
>>>>>>>>>> i
>>>>>
>>>>>> cant
>>>>>>>
>>>>>>>> figure it out:
>>>>>>>>>>>
>>>>>>>>>>> a) if that is all the data i need to set a reference system for
>>>>>>>>>>>
>>>>>>>>>> my
>>>>>
>>>>>> points
>>>>>>>
>>>>>>>> and the wrld_simp object.
>>>>>>>>>>>
>>>>>>>>>>> b) how to change a typical CRS object with such data
>>>>>>>>>>>
>>>>>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>>>>>>>>
>>>>>>>>>>> Where do i enter the TRUELAT and CENLAT values?
>>>>>>>>>>> Are there any site that explains easily what the fields in the
>>>>>>>>>>>
>>>>>>>>>> CRS
>>>>>
>>>>>> mean
>>>>>>>
>>>>>>>> and
>>>>>>>>>>
>>>>>>>>>>> how to change them?
>>>>>>>>>>>
>>>>>>>>>>> Thanks in advance.
>>>>>>>>>>>
>>>>>>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>>>>>>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>>>>>>>>
>>>>>>>>>> I believe:
>>>>>>>>>> +lat_0  = CEN_LAT   Latitude of origin
>>>>>>>>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>>>>>>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>>>>>>>> +lon_0  = CEN_LON   Central meridian
>>>>>>>>>>
>>>>>>>>>> proj strings are defined by the proj4 libary. It's website listed
>>>>>>>>>>
>>>>>>>>> above
>>>>>>
>>>>>>> and the associated mailing lists or gis stackexchange would be the
>>>>>>>>>> places to get help on it.
>>>>>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>>>>>>>>
>>>>>>>>>> It often helps to browse similar projections on
>>>>>>>>>> http://spatialreference.org/
>>>>>>>>>> http://epsg.io/
>>>>>>>>>>
>>>>>>>>>> Enjoy,
>>>>>>>>>> Alex
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>> _______________________________________________
>>>>>>>> R-sig-Geo mailing list
>>>>>>>> [hidden email]
>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>
>>>>>>> --
>>>>>>> Roger Bivand
>>>>>>> Department of Economics, Norwegian School of Economics,
>>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>>>> e-mail: [hidden email]
>>>>>>> http://orcid.org/0000-0003-2392-6140
>>>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>>>> http://depsy.org/person/434412
>>>>>>> _______________________________________________
>>>>>>> R-sig-Geo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>
>>>>>> --
>>>>>> Dr. Michael Sumner
>>>>>> Software and Database Engineer
>>>>>> Australian Antarctic Division
>>>>>> 203 Channel Highway
>>>>>> Kingston Tasmania 7050 Australia
>>>>>>
>>>>>>          [[alternative HTML version deleted]]
>>>>>>
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> Agustín Camacho Guerrero.
>>>>> Doutor em Zoologia.
>>>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>>>> Biociências, USP.
>>>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>>>
>>>>>          [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>>
>>>>
>>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Michael Sumner-2
In reply to this post by dschneiderch
On Wed, 24 Feb 2016 at 06:07 Dominik Schneider <
[hidden email]> wrote:

> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
> dealt with this.
> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
> +a=6370000 +b=6370000 +units=m
>
> +proj=lcc which is usually what wrf is run with.
> The tricky part is:
> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
> because every WRF run is different (the WRF Preprocessing System optimizes
> the projection for the domain).
> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>
> This gives:
> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
> +a=6370000 +b=6370000 +units=m +no_defs
>
>
This looks right to me:

library(raster)
library(rgdal)
prj <- "+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
+ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs"

r <- raster("results_us_future_output_none_0.nc", varname = "dctmx")
## use  print(r) to see the nc header dump

## lons/lats from NetCDF
lon <- raster("results_us_future_output_none_0.nc", varname = "lon")
lat <- raster("results_us_future_output_none_0.nc", varname = "lat")

## the true projected coordinates (??)
xy <- project(cbind(values(lon), values(lat)), prj)

## looks right
plot(xy, pch = ".")

## this fails, so it's not exactly right
##grd <- rasterFromXYZ(cbind(xy, 0)

## could use sp::points2grid tools to control tolerances, or just fudge it

ex <- extent(xy)
resolution <- c(round(max(diff(sort(unique(xy[,1]))))),
round(max(diff(sort(unique(xy[,2]))))))
ex <- ex + c(-1, 1, -1, 1) * rep(resolution / 2, 2)

## finally
mp <- setExtent(r, ex)
projection(mp) <- prj

library(maptools)
data(wrld_simpl)
namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
"Canada", "Mexico")), prj)

plot(mp)
plot(namer, add = TRUE)


("Why does this crucial metadata get thrown away?", he sighs . . .)

HTH




> But, wrf users like to give out lat and  long so you need to assign it:
> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>
> and then reproject the lat/long to lcc coordinates using this string:
> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0 +ellps=sphere
> +a=6370000 +b=6370000 +units=m +no_defs
>
> One word of caution, make sure you received the correct coordinates. Some
> variables are run cell center while some are run at cell edge. It looks
> like from your .nc file it was made by your collaborator so I assume they
> are right.
>
> That said, another word of caution, I found that the XLAT and XLONG
> variables from WRF output aren't very precise. There is a "geogrid" file
> from the preprocessing system that has the domain corners, resolution, nrow
> and ncol from which you can make a better grid using the native projection
> system (in my case it was a 4km grid). I suggest you try to get those.
>
> I hope this helps... I have to run but wanted to save people too much head
> scratching. I can get you running with more help tonight if you need.
> Dominik
>
>
> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
> wrote:
>
>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>> more giant crab this time.
>>
>> The creator of the .nc file also looked at this webpage:
>> http://www.pkrc.net/wrf-lambert.html
>> It seemed like the right proj4 string might be this one:
>>
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>     +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>
>> However this projection also does not allow me to adequately plot the
>> locations on the raster.
>>
>> Here is the .nc file. it contains several layers.
>>
>>
>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>
>>
>>
>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>
>> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>> >
>> > > On Tue, 23 Feb 2016, Alex Mandel wrote:
>> > >
>> > > > I made an attempt at it too. Investigating the original data, I'm
>> not
>> > > > sure that the projection information supplied is correct for the
>> data
>> > > > linked. When I load up the data in a unprojected space, the
>> coordinates
>> > > > don't look at all similar to any Lambert projected data I have, they
>> > > > actually look like Lat/Lon in some unprojected coordinate system,
>> > > > perhaps a different spheroid than expected.
>> > >
>> > > Does anyone have a link to the original data? Is is possible that
>> this is
>> > > the General Oblique Transformation used by modellers - that is
>> something
>> > > that feels like longlat but is recentred and oblique? Example at the
>> very
>> > > end of my GEOSTAT talk last year (slides 81-83):
>> > >
>> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf
>> > >
>> > > Roger
>> > >
>> > >
>> > For what it is worth, the General Oblique Transformation is not the only
>> > example - it's very common for modellers to have a mesh that has the
>> > "mostly-properties" of a projection, but is not actually describable
>> with
>> > standard transform + affine parameters. The main cases that I've seen
>> are
>> > polar stereographic, equal area or oblique Mercator. Often they really
>> are
>> > simple transforms and you can reconstruct without loss, but it's not
>> > usually possible to tell without exploration. It's an interesting
>> > dis-connect to see code that builds a mesh with certain properties, then
>> > only stores longitudes and latitudes - when it could be done with
>> standard
>> > tools and be stored and used much more efficiently.
>> >
>> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>> > terminology conflated in this context too. )
>> >
>> > I'm also interested to explore the original data.
>> >
>> > Cheers, Mike.
>> >
>> >
>> >
>> > > >
>> > > > -Alex
>> > > >
>> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>> > > >> Hi
>> > > >>
>> > > >> I tried to make it work but I had to give up. I wanted to reproject
>> > the
>> > > Lamberth conformal conic coordinates to long-lat but it didn't work.
>> > > >>
>> > > >> Perhaps someone can see what I did wrong. Here is what I did (data
>> in
>> > R
>> > > binary format and figure in png format both attached):
>> > > >>
>> > > >> library(raster)
>> > > >> library(maptools)
>> > > >> data(wrld_simpl)
>> > > >>
>> > > >> r <- raster("raster.grd")
>> > > >> projection(r)
>> > > >> ## > NA
>> > > >>
>> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>> ",")
>> > > >> coordinates(uro) <- ~lon+lat
>> > > >>
>> > > >> ## Set projections for the 3 data sets
>> > > >>
>> > > >> ## Lamberth's confocal conic projection with given parameters
>> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>> +lat_2=45.0
>> > > +ellps=WGS84"
>> > > >> projection(r)
>> > > >>
>> > > >> ## Assume that lon, lat are geographical coordinates (degrees
>> decimal)
>> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>> > > >>
>> > > >> ## wrld_simpl is in geographical coordinates
>> > > >> proj4string(wrld_simpl)
>> > > >>
>> > > >> ## Make figure in png format
>> > > >> ## Of course plotting data with 2 different projections will give
>> > > >> ## some distortions
>> > > >> pdf("uro.png")
>> > > >>
>> > > >> plot(r)
>> > > >> points(uro)
>> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>> 'r'
>> > > >>
>> > > >> dev.off()
>> > > >>
>> > > >> extent(r)
>> > > >> ## class       : Extent
>> > > >> ## xmin        : -131.4368
>> > > >> ## xmax        : -68.56323
>> > > >> ## ymin        : 12.35567
>> > > >> ## ymax        : 50.26619
>> > > >>
>> > > >> ## Reproject the raster to long-lat
>> > > >> ## This doesn't work (collapsed domain)
>> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>> > > +ellps=WGS84 +towgs84=0,0,0")
>> > > >>
>> > > >> ## Because
>> > > >>> extent(rp)
>> > > >> ## class       : Extent
>> > > >> ## xmin        : -100.0015
>> > > >> ## xmax        : -99.68557
>> > > >> ## ymin        : 37.70658
>> > > >> ## ymax        : 38.00046
>> > > >>
>> > > >> ## Save data in R binary format
>> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>> > > >>
>> > > >>
>> > > >> Yours sincerely / Med venlig hilsen
>> > > >>
>> > > >> Frede Aakmann Tøgersen
>> > > >> Specialist, M.Sc., Ph.D.
>> > > >> Plant Performance & Modeling
>> > > >>
>> > > >> Technology & Service Solutions
>> > > >> T +45 9730 5135
>> > > >> M +45 2547 6050
>> > > >> [hidden email]
>> > > >> http://www.vestas.com
>> > > >>
>> > > >> Company reg. name: Vestas Wind Systems A/S
>> > > >> This e-mail is subject to our e-mail disclaimer statement.
>> > > >> Please refer to www.vestas.com/legal/notice
>> > > >> If you have received this e-mail in error please contact the
>> sender.
>> > > >>
>> > > >>
>> > > >>
>> > > >> -----Original Message-----
>> > > >> From: R-sig-Geo [mailto:[hidden email]] On
>> Behalf Of
>> > > Agus Camacho
>> > > >> Sent: 22. februar 2016 19:20
>> > > >> To: [hidden email]
>> > > >> Cc: r-sig-geo
>> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to a
>> > > reference system implicit in a .nc file
>> > > >>
>> > > >> Thanks Alex, but the locations still fall in the sea when i plot
>> them
>> > > using
>> > > >> your recommended Solution. I looked at the sites you proposed and
>> they
>> > > have
>> > > >> other values for lat_1, lat_0, etc..
>> > > >>
>> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>> > > >>
>> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>> > > >>>> Dear all,
>> > > >>>>
>> > > >>>> Im trying to overlap these points:
>> > > >>>>
>> > > >>>
>> > >
>> >
>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>> > > >>>>
>> > > >>>> and a wrld_simpl object:
>> > > >>>> library(maptools)
>> > > >>>> data(wrld_simpl)
>> > > >>>>
>> > > >>>> Over this raster layer
>> > > >>>>
>> > > >>>
>> > >
>> >
>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>> > > >>>>
>> > > >>>> This rastr comes from a .nc file without a reference system. The
>> > > author
>> > > >>> of
>> > > >>>> that .nc file gave me the following data about the .nc.
>> > > >>>>
>> > > >>>> The projection is *Lambert conformal conic* projection
>> > > >>>> CEN_LAT = 38.0
>> > > >>>> CEN_LON = -100.0
>> > > >>>> TRUELAT1 = 25.
>> > > >>>> TRUELAT2 = 45.
>> > > >>>>
>> > > >>>> However, despite i have gone through many sites in the internet,
>> i
>> > > cant
>> > > >>>> figure it out:
>> > > >>>>
>> > > >>>> a) if that is all the data i need to set a reference system for
>> my
>> > > points
>> > > >>>> and the wrld_simp object.
>> > > >>>>
>> > > >>>> b) how to change a typical CRS object with such data
>> > > >>>>
>> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>> > > >>>>
>> > > >>>> Where do i enter the TRUELAT and CENLAT values?
>> > > >>>> Are there any site that explains easily what the fields in the
>> CRS
>> > > mean
>> > > >>> and
>> > > >>>> how to change them?
>> > > >>>>
>> > > >>>> Thanks in advance.
>> > > >>>>
>> > > >>>
>> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
>> > > >>> https://trac.osgeo.org/proj/wiki/GenParms
>> > > >>>
>> > > >>> I believe:
>> > > >>> +lat_0  = CEN_LAT   Latitude of origin
>> > > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>> > > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>> > > >>> +lon_0  = CEN_LON   Central meridian
>> > > >>>
>> > > >>> proj strings are defined by the proj4 libary. It's website listed
>> > above
>> > > >>> and the associated mailing lists or gis stackexchange would be the
>> > > >>> places to get help on it.
>> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
>> > > >>>
>> > > >>> It often helps to browse similar projections on
>> > > >>> http://spatialreference.org/
>> > > >>> http://epsg.io/
>> > > >>>
>> > > >>> Enjoy,
>> > > >>> Alex
>> > > >>>
>> > > >>>
>> > > >>
>> > > >>
>> > > >
>> > > > _______________________________________________
>> > > > R-sig-Geo mailing list
>> > > > [hidden email]
>> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> > >
>> > > --
>> > > Roger Bivand
>> > > Department of Economics, Norwegian School of Economics,
>> > > Helleveien 30, N-5045 Bergen, Norway.
>> > > voice: +47 55 95 93 55; fax +47 55 95 91 00
>> > > e-mail: [hidden email]
>> > > http://orcid.org/0000-0003-2392-6140
>> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>> > > http://depsy.org/person/434412
>> > > _______________________________________________
>> > > R-sig-Geo mailing list
>> > > [hidden email]
>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>> > --
>> > Dr. Michael Sumner
>> > Software and Database Engineer
>> > Australian Antarctic Division
>> > 203 Channel Highway
>> > Kingston Tasmania 7050 Australia
>> >
>> >         [[alternative HTML version deleted]]
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > [hidden email]
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>>
>>
>> --
>> Agustín Camacho Guerrero.
>> Doutor em Zoologia.
>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>> Biociências, USP.
>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>> São Paulo - SP, CEP: 05508-090, Brasil.
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
Thanks Michael, that was a good job and despite the map looks nice now, the
locations cant still be plotted adequately, either as raw lon lat or as
spatial points with different projections.


x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
x=x[,1:3]
colnames(x)
coordinates(x)=~lon+lat
projection(x)="+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
x=spTransform(x, CRS("+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
+lon_0=-100.0 +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs"))
points(x)

2016-02-23 16:24 GMT-07:00 Michael Sumner <[hidden email]>:

>
>
> On Wed, 24 Feb 2016 at 06:07 Dominik Schneider <
> [hidden email]> wrote:
>
>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>> dealt with this.
>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>> +a=6370000 +b=6370000 +units=m
>>
>> +proj=lcc which is usually what wrf is run with.
>> The tricky part is:
>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> because every WRF run is different (the WRF Preprocessing System
>> optimizes the projection for the domain).
>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>
>> This gives:
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>>
> This looks right to me:
>
> library(raster)
> library(rgdal)
> prj <- "+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs"
>
> r <- raster("results_us_future_output_none_0.nc", varname = "dctmx")
> ## use  print(r) to see the nc header dump
>
> ## lons/lats from NetCDF
> lon <- raster("results_us_future_output_none_0.nc", varname = "lon")
> lat <- raster("results_us_future_output_none_0.nc", varname = "lat")
>
> ## the true projected coordinates (??)
> xy <- project(cbind(values(lon), values(lat)), prj)
>
> ## looks right
> plot(xy, pch = ".")
>
> ## this fails, so it's not exactly right
> ##grd <- rasterFromXYZ(cbind(xy, 0)
>
> ## could use sp::points2grid tools to control tolerances, or just fudge it
>
> ex <- extent(xy)
> resolution <- c(round(max(diff(sort(unique(xy[,1]))))),
> round(max(diff(sort(unique(xy[,2]))))))
> ex <- ex + c(-1, 1, -1, 1) * rep(resolution / 2, 2)
>
> ## finally
> mp <- setExtent(r, ex)
> projection(mp) <- prj
>
> library(maptools)
> data(wrld_simpl)
> namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
> "Canada", "Mexico")), prj)
>
> plot(mp)
> plot(namer, add = TRUE)
>
>
> ("Why does this crucial metadata get thrown away?", he sighs . . .)
>
> HTH
>
>
>
>
>> But, wrf users like to give out lat and  long so you need to assign it:
>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>> and then reproject the lat/long to lcc coordinates using this string:
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>> One word of caution, make sure you received the correct coordinates. Some
>> variables are run cell center while some are run at cell edge. It looks
>> like from your .nc file it was made by your collaborator so I assume they
>> are right.
>>
>> That said, another word of caution, I found that the XLAT and XLONG
>> variables from WRF output aren't very precise. There is a "geogrid" file
>> from the preprocessing system that has the domain corners, resolution, nrow
>> and ncol from which you can make a better grid using the native projection
>> system (in my case it was a 4km grid). I suggest you try to get those.
>>
>> I hope this helps... I have to run but wanted to save people too much
>> head scratching. I can get you running with more help tonight if you need.
>> Dominik
>>
>>
>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>> wrote:
>>
>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>>> more giant crab this time.
>>>
>>> The creator of the .nc file also looked at this webpage:
>>> http://www.pkrc.net/wrf-lambert.html
>>> It seemed like the right proj4 string might be this one:
>>>
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>     +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>
>>> However this projection also does not allow me to adequately plot the
>>> locations on the raster.
>>>
>>> Here is the .nc file. it contains several layers.
>>>
>>>
>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>
>>>
>>>
>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>
>>> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>> >
>>> > > On Tue, 23 Feb 2016, Alex Mandel wrote:
>>> > >
>>> > > > I made an attempt at it too. Investigating the original data, I'm
>>> not
>>> > > > sure that the projection information supplied is correct for the
>>> data
>>> > > > linked. When I load up the data in a unprojected space, the
>>> coordinates
>>> > > > don't look at all similar to any Lambert projected data I have,
>>> they
>>> > > > actually look like Lat/Lon in some unprojected coordinate system,
>>> > > > perhaps a different spheroid than expected.
>>> > >
>>> > > Does anyone have a link to the original data? Is is possible that
>>> this is
>>> > > the General Oblique Transformation used by modellers - that is
>>> something
>>> > > that feels like longlat but is recentred and oblique? Example at the
>>> very
>>> > > end of my GEOSTAT talk last year (slides 81-83):
>>> > >
>>> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>> > >
>>> > > Roger
>>> > >
>>> > >
>>> > For what it is worth, the General Oblique Transformation is not the
>>> only
>>> > example - it's very common for modellers to have a mesh that has the
>>> > "mostly-properties" of a projection, but is not actually describable
>>> with
>>> > standard transform + affine parameters. The main cases that I've seen
>>> are
>>> > polar stereographic, equal area or oblique Mercator. Often they really
>>> are
>>> > simple transforms and you can reconstruct without loss, but it's not
>>> > usually possible to tell without exploration. It's an interesting
>>> > dis-connect to see code that builds a mesh with certain properties,
>>> then
>>> > only stores longitudes and latitudes - when it could be done with
>>> standard
>>> > tools and be stored and used much more efficiently.
>>> >
>>> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>> > terminology conflated in this context too. )
>>> >
>>> > I'm also interested to explore the original data.
>>> >
>>> > Cheers, Mike.
>>> >
>>> >
>>> >
>>> > > >
>>> > > > -Alex
>>> > > >
>>> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>> > > >> Hi
>>> > > >>
>>> > > >> I tried to make it work but I had to give up. I wanted to
>>> reproject
>>> > the
>>> > > Lamberth conformal conic coordinates to long-lat but it didn't work.
>>> > > >>
>>> > > >> Perhaps someone can see what I did wrong. Here is what I did
>>> (data in
>>> > R
>>> > > binary format and figure in png format both attached):
>>> > > >>
>>> > > >> library(raster)
>>> > > >> library(maptools)
>>> > > >> data(wrld_simpl)
>>> > > >>
>>> > > >> r <- raster("raster.grd")
>>> > > >> projection(r)
>>> > > >> ## > NA
>>> > > >>
>>> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>> ",")
>>> > > >> coordinates(uro) <- ~lon+lat
>>> > > >>
>>> > > >> ## Set projections for the 3 data sets
>>> > > >>
>>> > > >> ## Lamberth's confocal conic projection with given parameters
>>> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>> +lat_2=45.0
>>> > > +ellps=WGS84"
>>> > > >> projection(r)
>>> > > >>
>>> > > >> ## Assume that lon, lat are geographical coordinates (degrees
>>> decimal)
>>> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>> > > >>
>>> > > >> ## wrld_simpl is in geographical coordinates
>>> > > >> proj4string(wrld_simpl)
>>> > > >>
>>> > > >> ## Make figure in png format
>>> > > >> ## Of course plotting data with 2 different projections will give
>>> > > >> ## some distortions
>>> > > >> pdf("uro.png")
>>> > > >>
>>> > > >> plot(r)
>>> > > >> points(uro)
>>> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>> 'r'
>>> > > >>
>>> > > >> dev.off()
>>> > > >>
>>> > > >> extent(r)
>>> > > >> ## class       : Extent
>>> > > >> ## xmin        : -131.4368
>>> > > >> ## xmax        : -68.56323
>>> > > >> ## ymin        : 12.35567
>>> > > >> ## ymax        : 50.26619
>>> > > >>
>>> > > >> ## Reproject the raster to long-lat
>>> > > >> ## This doesn't work (collapsed domain)
>>> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>> > > +ellps=WGS84 +towgs84=0,0,0")
>>> > > >>
>>> > > >> ## Because
>>> > > >>> extent(rp)
>>> > > >> ## class       : Extent
>>> > > >> ## xmin        : -100.0015
>>> > > >> ## xmax        : -99.68557
>>> > > >> ## ymin        : 37.70658
>>> > > >> ## ymax        : 38.00046
>>> > > >>
>>> > > >> ## Save data in R binary format
>>> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>> > > >>
>>> > > >>
>>> > > >> Yours sincerely / Med venlig hilsen
>>> > > >>
>>> > > >> Frede Aakmann Tøgersen
>>> > > >> Specialist, M.Sc., Ph.D.
>>> > > >> Plant Performance & Modeling
>>> > > >>
>>> > > >> Technology & Service Solutions
>>> > > >> T +45 9730 5135
>>> > > >> M +45 2547 6050
>>> > > >> [hidden email]
>>> > > >> http://www.vestas.com
>>> > > >>
>>> > > >> Company reg. name: Vestas Wind Systems A/S
>>> > > >> This e-mail is subject to our e-mail disclaimer statement.
>>> > > >> Please refer to www.vestas.com/legal/notice
>>> > > >> If you have received this e-mail in error please contact the
>>> sender.
>>> > > >>
>>> > > >>
>>> > > >>
>>> > > >> -----Original Message-----
>>> > > >> From: R-sig-Geo [mailto:[hidden email]] On
>>> Behalf Of
>>> > > Agus Camacho
>>> > > >> Sent: 22. februar 2016 19:20
>>> > > >> To: [hidden email]
>>> > > >> Cc: r-sig-geo
>>> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to
>>> a
>>> > > reference system implicit in a .nc file
>>> > > >>
>>> > > >> Thanks Alex, but the locations still fall in the sea when i plot
>>> them
>>> > > using
>>> > > >> your recommended Solution. I looked at the sites you proposed and
>>> they
>>> > > have
>>> > > >> other values for lat_1, lat_0, etc..
>>> > > >>
>>> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>> > > >>
>>> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>> > > >>>> Dear all,
>>> > > >>>>
>>> > > >>>> Im trying to overlap these points:
>>> > > >>>>
>>> > > >>>
>>> > >
>>> >
>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>> > > >>>>
>>> > > >>>> and a wrld_simpl object:
>>> > > >>>> library(maptools)
>>> > > >>>> data(wrld_simpl)
>>> > > >>>>
>>> > > >>>> Over this raster layer
>>> > > >>>>
>>> > > >>>
>>> > >
>>> >
>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>> > > >>>>
>>> > > >>>> This rastr comes from a .nc file without a reference system. The
>>> > > author
>>> > > >>> of
>>> > > >>>> that .nc file gave me the following data about the .nc.
>>> > > >>>>
>>> > > >>>> The projection is *Lambert conformal conic* projection
>>> > > >>>> CEN_LAT = 38.0
>>> > > >>>> CEN_LON = -100.0
>>> > > >>>> TRUELAT1 = 25.
>>> > > >>>> TRUELAT2 = 45.
>>> > > >>>>
>>> > > >>>> However, despite i have gone through many sites in the
>>> internet, i
>>> > > cant
>>> > > >>>> figure it out:
>>> > > >>>>
>>> > > >>>> a) if that is all the data i need to set a reference system for
>>> my
>>> > > points
>>> > > >>>> and the wrld_simp object.
>>> > > >>>>
>>> > > >>>> b) how to change a typical CRS object with such data
>>> > > >>>>
>>> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>> > > >>>>
>>> > > >>>> Where do i enter the TRUELAT and CENLAT values?
>>> > > >>>> Are there any site that explains easily what the fields in the
>>> CRS
>>> > > mean
>>> > > >>> and
>>> > > >>>> how to change them?
>>> > > >>>>
>>> > > >>>> Thanks in advance.
>>> > > >>>>
>>> > > >>>
>>> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>> > > >>> https://trac.osgeo.org/proj/wiki/GenParms
>>> > > >>>
>>> > > >>> I believe:
>>> > > >>> +lat_0  = CEN_LAT   Latitude of origin
>>> > > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>> > > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>> > > >>> +lon_0  = CEN_LON   Central meridian
>>> > > >>>
>>> > > >>> proj strings are defined by the proj4 libary. It's website listed
>>> > above
>>> > > >>> and the associated mailing lists or gis stackexchange would be
>>> the
>>> > > >>> places to get help on it.
>>> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>> > > >>>
>>> > > >>> It often helps to browse similar projections on
>>> > > >>> http://spatialreference.org/
>>> > > >>> http://epsg.io/
>>> > > >>>
>>> > > >>> Enjoy,
>>> > > >>> Alex
>>> > > >>>
>>> > > >>>
>>> > > >>
>>> > > >>
>>> > > >
>>> > > > _______________________________________________
>>> > > > R-sig-Geo mailing list
>>> > > > [hidden email]
>>> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> > >
>>> > > --
>>> > > Roger Bivand
>>> > > Department of Economics, Norwegian School of Economics,
>>> > > Helleveien 30, N-5045 Bergen, Norway.
>>> > > voice: +47 55 95 93 55; fax +47 55 95 91 00
>>> > > e-mail: [hidden email]
>>> > > http://orcid.org/0000-0003-2392-6140
>>> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>> > > http://depsy.org/person/434412
>>> > > _______________________________________________
>>> > > R-sig-Geo mailing list
>>> > > [hidden email]
>>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> >
>>> > --
>>> > Dr. Michael Sumner
>>> > Software and Database Engineer
>>> > Australian Antarctic Division
>>> > 203 Channel Highway
>>> > Kingston Tasmania 7050 Australia
>>> >
>>> >         [[alternative HTML version deleted]]
>>> >
>>> > _______________________________________________
>>> > R-sig-Geo mailing list
>>> > [hidden email]
>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> >
>>>
>>>
>>>
>>> --
>>> Agustín Camacho Guerrero.
>>> Doutor em Zoologia.
>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>> Biociências, USP.
>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
>


--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

dschneiderch
just glancing through the email chain. looks like you were able to get
things plotted....but you are still having an issue with some points?

here is a link for information regarding the grid staggering used in WRF.
http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm

> Some remarks:
> (1) I just took the first pair of  coordinates as derived from gdalinfo("
> geo_em.d01.nc")
> you will find 4 different coordinate pairs (i did not proof which one is
> right
> The data is staggered (as outlined by Dominik) So some of the corner
> coordinates belongs to the staggered data and the others coordinates  to
> the unstaggered ones.



On Tue, Feb 23, 2016 at 4:30 PM, Agus Camacho <[hidden email]>
wrote:

> Thanks Michael, that was a good job and despite the map looks nice now,
> the locations cant still be plotted adequately, either as raw lon lat or as
> spatial points with different projections.
>
>
> x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
> Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
> x=x[,1:3]
> colnames(x)
> coordinates(x)=~lon+lat
> projection(x)="+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
> x=spTransform(x, CRS("+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
> +lon_0=-100.0 +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs"))
> points(x)
>
> 2016-02-23 16:24 GMT-07:00 Michael Sumner <[hidden email]>:
>
>>
>>
>> On Wed, 24 Feb 2016 at 06:07 Dominik Schneider <
>> [hidden email]> wrote:
>>
>>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>>> dealt with this.
>>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m
>>>
>>> +proj=lcc which is usually what wrf is run with.
>>> The tricky part is:
>>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> because every WRF run is different (the WRF Preprocessing System
>>> optimizes the projection for the domain).
>>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>>
>>> This gives:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>>
>> This looks right to me:
>>
>> library(raster)
>> library(rgdal)
>> prj <- "+proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs"
>>
>> r <- raster("results_us_future_output_none_0.nc", varname = "dctmx")
>> ## use  print(r) to see the nc header dump
>>
>> ## lons/lats from NetCDF
>> lon <- raster("results_us_future_output_none_0.nc", varname = "lon")
>> lat <- raster("results_us_future_output_none_0.nc", varname = "lat")
>>
>> ## the true projected coordinates (??)
>> xy <- project(cbind(values(lon), values(lat)), prj)
>>
>> ## looks right
>> plot(xy, pch = ".")
>>
>> ## this fails, so it's not exactly right
>> ##grd <- rasterFromXYZ(cbind(xy, 0)
>>
>> ## could use sp::points2grid tools to control tolerances, or just fudge it
>>
>> ex <- extent(xy)
>> resolution <- c(round(max(diff(sort(unique(xy[,1]))))),
>> round(max(diff(sort(unique(xy[,2]))))))
>> ex <- ex + c(-1, 1, -1, 1) * rep(resolution / 2, 2)
>>
>> ## finally
>> mp <- setExtent(r, ex)
>> projection(mp) <- prj
>>
>> library(maptools)
>> data(wrld_simpl)
>> namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
>> "Canada", "Mexico")), prj)
>>
>> plot(mp)
>> plot(namer, add = TRUE)
>>
>>
>> ("Why does this crucial metadata get thrown away?", he sighs . . .)
>>
>> HTH
>>
>>
>>
>>
>>> But, wrf users like to give out lat and  long so you need to assign it:
>>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> and then reproject the lat/long to lcc coordinates using this string:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> One word of caution, make sure you received the correct coordinates.
>>> Some variables are run cell center while some are run at cell edge. It
>>> looks like from your .nc file it was made by your collaborator so I assume
>>> they are right.
>>>
>>> That said, another word of caution, I found that the XLAT and XLONG
>>> variables from WRF output aren't very precise. There is a "geogrid" file
>>> from the preprocessing system that has the domain corners, resolution, nrow
>>> and ncol from which you can make a better grid using the native projection
>>> system (in my case it was a 4km grid). I suggest you try to get those.
>>>
>>> I hope this helps... I have to run but wanted to save people too much
>>> head scratching. I can get you running with more help tonight if you need.
>>> Dominik
>>>
>>>
>>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>>> wrote:
>>>
>>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill
>>>> bring
>>>> more giant crab this time.
>>>>
>>>> The creator of the .nc file also looked at this webpage:
>>>> http://www.pkrc.net/wrf-lambert.html
>>>> It seemed like the right proj4 string might be this one:
>>>>
>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>>     +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>>
>>>> However this projection also does not allow me to adequately plot the
>>>> locations on the raster.
>>>>
>>>> Here is the .nc file. it contains several layers.
>>>>
>>>>
>>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>>
>>>>
>>>>
>>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>>
>>>> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]>
>>>> wrote:
>>>> >
>>>> > > On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>> > >
>>>> > > > I made an attempt at it too. Investigating the original data, I'm
>>>> not
>>>> > > > sure that the projection information supplied is correct for the
>>>> data
>>>> > > > linked. When I load up the data in a unprojected space, the
>>>> coordinates
>>>> > > > don't look at all similar to any Lambert projected data I have,
>>>> they
>>>> > > > actually look like Lat/Lon in some unprojected coordinate system,
>>>> > > > perhaps a different spheroid than expected.
>>>> > >
>>>> > > Does anyone have a link to the original data? Is is possible that
>>>> this is
>>>> > > the General Oblique Transformation used by modellers - that is
>>>> something
>>>> > > that feels like longlat but is recentred and oblique? Example at
>>>> the very
>>>> > > end of my GEOSTAT talk last year (slides 81-83):
>>>> > >
>>>> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>> > >
>>>> > > Roger
>>>> > >
>>>> > >
>>>> > For what it is worth, the General Oblique Transformation is not the
>>>> only
>>>> > example - it's very common for modellers to have a mesh that has the
>>>> > "mostly-properties" of a projection, but is not actually describable
>>>> with
>>>> > standard transform + affine parameters. The main cases that I've seen
>>>> are
>>>> > polar stereographic, equal area or oblique Mercator. Often they
>>>> really are
>>>> > simple transforms and you can reconstruct without loss, but it's not
>>>> > usually possible to tell without exploration. It's an interesting
>>>> > dis-connect to see code that builds a mesh with certain properties,
>>>> then
>>>> > only stores longitudes and latitudes - when it could be done with
>>>> standard
>>>> > tools and be stored and used much more efficiently.
>>>> >
>>>> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>>> > terminology conflated in this context too. )
>>>> >
>>>> > I'm also interested to explore the original data.
>>>> >
>>>> > Cheers, Mike.
>>>> >
>>>> >
>>>> >
>>>> > > >
>>>> > > > -Alex
>>>> > > >
>>>> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>> > > >> Hi
>>>> > > >>
>>>> > > >> I tried to make it work but I had to give up. I wanted to
>>>> reproject
>>>> > the
>>>> > > Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>> > > >>
>>>> > > >> Perhaps someone can see what I did wrong. Here is what I did
>>>> (data in
>>>> > R
>>>> > > binary format and figure in png format both attached):
>>>> > > >>
>>>> > > >> library(raster)
>>>> > > >> library(maptools)
>>>> > > >> data(wrld_simpl)
>>>> > > >>
>>>> > > >> r <- raster("raster.grd")
>>>> > > >> projection(r)
>>>> > > >> ## > NA
>>>> > > >>
>>>> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>>> ",")
>>>> > > >> coordinates(uro) <- ~lon+lat
>>>> > > >>
>>>> > > >> ## Set projections for the 3 data sets
>>>> > > >>
>>>> > > >> ## Lamberth's confocal conic projection with given parameters
>>>> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>>> +lat_2=45.0
>>>> > > +ellps=WGS84"
>>>> > > >> projection(r)
>>>> > > >>
>>>> > > >> ## Assume that lon, lat are geographical coordinates (degrees
>>>> decimal)
>>>> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>> > > >>
>>>> > > >> ## wrld_simpl is in geographical coordinates
>>>> > > >> proj4string(wrld_simpl)
>>>> > > >>
>>>> > > >> ## Make figure in png format
>>>> > > >> ## Of course plotting data with 2 different projections will give
>>>> > > >> ## some distortions
>>>> > > >> pdf("uro.png")
>>>> > > >>
>>>> > > >> plot(r)
>>>> > > >> points(uro)
>>>> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent
>>>> of 'r'
>>>> > > >>
>>>> > > >> dev.off()
>>>> > > >>
>>>> > > >> extent(r)
>>>> > > >> ## class       : Extent
>>>> > > >> ## xmin        : -131.4368
>>>> > > >> ## xmax        : -68.56323
>>>> > > >> ## ymin        : 12.35567
>>>> > > >> ## ymax        : 50.26619
>>>> > > >>
>>>> > > >> ## Reproject the raster to long-lat
>>>> > > >> ## This doesn't work (collapsed domain)
>>>> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>>> > > +ellps=WGS84 +towgs84=0,0,0")
>>>> > > >>
>>>> > > >> ## Because
>>>> > > >>> extent(rp)
>>>> > > >> ## class       : Extent
>>>> > > >> ## xmin        : -100.0015
>>>> > > >> ## xmax        : -99.68557
>>>> > > >> ## ymin        : 37.70658
>>>> > > >> ## ymax        : 38.00046
>>>> > > >>
>>>> > > >> ## Save data in R binary format
>>>> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>> > > >>
>>>> > > >>
>>>> > > >> Yours sincerely / Med venlig hilsen
>>>> > > >>
>>>> > > >> Frede Aakmann Tøgersen
>>>> > > >> Specialist, M.Sc., Ph.D.
>>>> > > >> Plant Performance & Modeling
>>>> > > >>
>>>> > > >> Technology & Service Solutions
>>>> > > >> T +45 9730 5135
>>>> > > >> M +45 2547 6050
>>>> > > >> [hidden email]
>>>> > > >> http://www.vestas.com
>>>> > > >>
>>>> > > >> Company reg. name: Vestas Wind Systems A/S
>>>> > > >> This e-mail is subject to our e-mail disclaimer statement.
>>>> > > >> Please refer to www.vestas.com/legal/notice
>>>> > > >> If you have received this e-mail in error please contact the
>>>> sender.
>>>> > > >>
>>>> > > >>
>>>> > > >>
>>>> > > >> -----Original Message-----
>>>> > > >> From: R-sig-Geo [mailto:[hidden email]] On
>>>> Behalf Of
>>>> > > Agus Camacho
>>>> > > >> Sent: 22. februar 2016 19:20
>>>> > > >> To: [hidden email]
>>>> > > >> Cc: r-sig-geo
>>>> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl
>>>> to a
>>>> > > reference system implicit in a .nc file
>>>> > > >>
>>>> > > >> Thanks Alex, but the locations still fall in the sea when i plot
>>>> them
>>>> > > using
>>>> > > >> your recommended Solution. I looked at the sites you proposed
>>>> and they
>>>> > > have
>>>> > > >> other values for lat_1, lat_0, etc..
>>>> > > >>
>>>> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>> > > >>
>>>> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>> > > >>>> Dear all,
>>>> > > >>>>
>>>> > > >>>> Im trying to overlap these points:
>>>> > > >>>>
>>>> > > >>>
>>>> > >
>>>> >
>>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>> > > >>>>
>>>> > > >>>> and a wrld_simpl object:
>>>> > > >>>> library(maptools)
>>>> > > >>>> data(wrld_simpl)
>>>> > > >>>>
>>>> > > >>>> Over this raster layer
>>>> > > >>>>
>>>> > > >>>
>>>> > >
>>>> >
>>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>> > > >>>>
>>>> > > >>>> This rastr comes from a .nc file without a reference system.
>>>> The
>>>> > > author
>>>> > > >>> of
>>>> > > >>>> that .nc file gave me the following data about the .nc.
>>>> > > >>>>
>>>> > > >>>> The projection is *Lambert conformal conic* projection
>>>> > > >>>> CEN_LAT = 38.0
>>>> > > >>>> CEN_LON = -100.0
>>>> > > >>>> TRUELAT1 = 25.
>>>> > > >>>> TRUELAT2 = 45.
>>>> > > >>>>
>>>> > > >>>> However, despite i have gone through many sites in the
>>>> internet, i
>>>> > > cant
>>>> > > >>>> figure it out:
>>>> > > >>>>
>>>> > > >>>> a) if that is all the data i need to set a reference system
>>>> for my
>>>> > > points
>>>> > > >>>> and the wrld_simp object.
>>>> > > >>>>
>>>> > > >>>> b) how to change a typical CRS object with such data
>>>> > > >>>>
>>>> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>> > > >>>>
>>>> > > >>>> Where do i enter the TRUELAT and CENLAT values?
>>>> > > >>>> Are there any site that explains easily what the fields in the
>>>> CRS
>>>> > > mean
>>>> > > >>> and
>>>> > > >>>> how to change them?
>>>> > > >>>>
>>>> > > >>>> Thanks in advance.
>>>> > > >>>>
>>>> > > >>>
>>>> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>> > > >>> https://trac.osgeo.org/proj/wiki/GenParms
>>>> > > >>>
>>>> > > >>> I believe:
>>>> > > >>> +lat_0  = CEN_LAT   Latitude of origin
>>>> > > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>> > > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>> > > >>> +lon_0  = CEN_LON   Central meridian
>>>> > > >>>
>>>> > > >>> proj strings are defined by the proj4 libary. It's website
>>>> listed
>>>> > above
>>>> > > >>> and the associated mailing lists or gis stackexchange would be
>>>> the
>>>> > > >>> places to get help on it.
>>>> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>> > > >>>
>>>> > > >>> It often helps to browse similar projections on
>>>> > > >>> http://spatialreference.org/
>>>> > > >>> http://epsg.io/
>>>> > > >>>
>>>> > > >>> Enjoy,
>>>> > > >>> Alex
>>>> > > >>>
>>>> > > >>>
>>>> > > >>
>>>> > > >>
>>>> > > >
>>>> > > > _______________________________________________
>>>> > > > R-sig-Geo mailing list
>>>> > > > [hidden email]
>>>> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> > >
>>>> > > --
>>>> > > Roger Bivand
>>>> > > Department of Economics, Norwegian School of Economics,
>>>> > > Helleveien 30, N-5045 Bergen, Norway.
>>>> > > voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>> > > e-mail: [hidden email]
>>>> > > http://orcid.org/0000-0003-2392-6140
>>>> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>> > > http://depsy.org/person/434412
>>>> > > _______________________________________________
>>>> > > R-sig-Geo mailing list
>>>> > > [hidden email]
>>>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> >
>>>> > --
>>>> > Dr. Michael Sumner
>>>> > Software and Database Engineer
>>>> > Australian Antarctic Division
>>>> > 203 Channel Highway
>>>> > Kingston Tasmania 7050 Australia
>>>> >
>>>> >         [[alternative HTML version deleted]]
>>>> >
>>>> > _______________________________________________
>>>> > R-sig-Geo mailing list
>>>> > [hidden email]
>>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>> >
>>>>
>>>>
>>>>
>>>> --
>>>> Agustín Camacho Guerrero.
>>>> Doutor em Zoologia.
>>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>>> Biociências, USP.
>>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>> --
>> Dr. Michael Sumner
>> Software and Database Engineer
>> Australian Antarctic Division
>> 203 Channel Highway
>> Kingston Tasmania 7050 Australia
>>
>>
>
>
> --
> Agustín Camacho Guerrero.
> Doutor em Zoologia.
> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
> Biociências, USP.
> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
> São Paulo - SP, CEP: 05508-090, Brasil.
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

agus
In reply to this post by agus
With the help of everybody i finally got this to work, so here is a script
that does the job of reprojecting both, a raster layer obtained from a .nc
and some locations in order to overplot them using plot.raster or mapview.
Im using a combination of the advices of Dominik, Michael and Chris.


require(ncdf4)
require(raster)


setwd("C:/~")


r=raster("geo_em.d01.nc",
              varname="dctmx")# days of ctmax events

# Set extent and projections of rasters for plotting
# chris gave me the orig data fom the nc file because i could not install
gdal
xmin= -3545999
xmax= 3546000
ymin = -2286000
ymax=2286000
pr<- "+proj=lcc +lat_1=25 +lat_2=45 +lat_0=38.000008 +lon_0=-100 +x_0=0
+y_0=0"

wrfLccExt<-extent(xmin,xmax,ymin,ymax)

extent(r) <- extent(wrfLccExt)
projection(r) <- pr

# get and prepare  urosaurus locations

x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
x=x[,1:3]
colnames(x)
coordinates(x)=~lon+lat
proj4string(x)<- CRS("+proj=longlat +datum=WGS84")
x=spTransform(x, pr)


# get prepare and plot wrld_simpl
require(maptools)
require(sp)
data(wrld_simpl)
namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
"Canada", "Mexico")), prj)

#plot with raster
plot(r, cex.main=.7,legend=F)
points(x)
plot(namer, add = TRUE)

# plot with mapview (cool!)
m=mapview(r)
u=mapview(x)
m+u


Thanks to all!
Agus

2016-02-23 13:11 GMT-07:00 Agus Camacho <[hidden email]>:

> Thanks for that Dominik,
>
> Giving that projection to either the locations, the raster layer generated
> from the .nc file, or both, still did not work. I keep having locations
> that should be on land falling far on the sea. Might this be a problem
> derived from using raster with a file whose original grid distances are not
> constant?
>
> Here is a link with the original file which has the original coordinate
> data.
>
> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
>
>
> 2016-02-23 12:07 GMT-07:00 Dominik Schneider <
> [hidden email]>:
>
>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>> dealt with this.
>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>> +a=6370000 +b=6370000 +units=m
>>
>> +proj=lcc which is usually what wrf is run with.
>> The tricky part is:
>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> because every WRF run is different (the WRF Preprocessing System
>> optimizes the projection for the domain).
>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>
>> This gives:
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>> But, wrf users like to give out lat and  long so you need to assign it:
>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>> and then reproject the lat/long to lcc coordinates using this string:
>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>
>> One word of caution, make sure you received the correct coordinates. Some
>> variables are run cell center while some are run at cell edge. It looks
>> like from your .nc file it was made by your collaborator so I assume they
>> are right.
>>
>> That said, another word of caution, I found that the XLAT and XLONG
>> variables from WRF output aren't very precise. There is a "geogrid" file
>> from the preprocessing system that has the domain corners, resolution, nrow
>> and ncol from which you can make a better grid using the native projection
>> system (in my case it was a 4km grid). I suggest you try to get those.
>>
>> I hope this helps... I have to run but wanted to save people too much
>> head scratching. I can get you running with more help tonight if you need.
>> Dominik
>>
>>
>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>> wrote:
>>
>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>>> more giant crab this time.
>>>
>>> The creator of the .nc file also looked at this webpage:
>>> http://www.pkrc.net/wrf-lambert.html
>>> It seemed like the right proj4 string might be this one:
>>>
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>     +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>
>>> However this projection also does not allow me to adequately plot the
>>> locations on the raster.
>>>
>>> Here is the .nc file. it contains several layers.
>>>
>>>
>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>
>>>
>>>
>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>
>>> > On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>> >
>>> > > On Tue, 23 Feb 2016, Alex Mandel wrote:
>>> > >
>>> > > > I made an attempt at it too. Investigating the original data, I'm
>>> not
>>> > > > sure that the projection information supplied is correct for the
>>> data
>>> > > > linked. When I load up the data in a unprojected space, the
>>> coordinates
>>> > > > don't look at all similar to any Lambert projected data I have,
>>> they
>>> > > > actually look like Lat/Lon in some unprojected coordinate system,
>>> > > > perhaps a different spheroid than expected.
>>> > >
>>> > > Does anyone have a link to the original data? Is is possible that
>>> this is
>>> > > the General Oblique Transformation used by modellers - that is
>>> something
>>> > > that feels like longlat but is recentred and oblique? Example at the
>>> very
>>> > > end of my GEOSTAT talk last year (slides 81-83):
>>> > >
>>> > > http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>> > >
>>> > > Roger
>>> > >
>>> > >
>>> > For what it is worth, the General Oblique Transformation is not the
>>> only
>>> > example - it's very common for modellers to have a mesh that has the
>>> > "mostly-properties" of a projection, but is not actually describable
>>> with
>>> > standard transform + affine parameters. The main cases that I've seen
>>> are
>>> > polar stereographic, equal area or oblique Mercator. Often they really
>>> are
>>> > simple transforms and you can reconstruct without loss, but it's not
>>> > usually possible to tell without exploration. It's an interesting
>>> > dis-connect to see code that builds a mesh with certain properties,
>>> then
>>> > only stores longitudes and latitudes - when it could be done with
>>> standard
>>> > tools and be stored and used much more efficiently.
>>> >
>>> > (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>> > terminology conflated in this context too. )
>>> >
>>> > I'm also interested to explore the original data.
>>> >
>>> > Cheers, Mike.
>>> >
>>> >
>>> >
>>> > > >
>>> > > > -Alex
>>> > > >
>>> > > > On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>> > > >> Hi
>>> > > >>
>>> > > >> I tried to make it work but I had to give up. I wanted to
>>> reproject
>>> > the
>>> > > Lamberth conformal conic coordinates to long-lat but it didn't work.
>>> > > >>
>>> > > >> Perhaps someone can see what I did wrong. Here is what I did
>>> (data in
>>> > R
>>> > > binary format and figure in png format both attached):
>>> > > >>
>>> > > >> library(raster)
>>> > > >> library(maptools)
>>> > > >> data(wrld_simpl)
>>> > > >>
>>> > > >> r <- raster("raster.grd")
>>> > > >> projection(r)
>>> > > >> ## > NA
>>> > > >>
>>> > > >> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>> ",")
>>> > > >> coordinates(uro) <- ~lon+lat
>>> > > >>
>>> > > >> ## Set projections for the 3 data sets
>>> > > >>
>>> > > >> ## Lamberth's confocal conic projection with given parameters
>>> > > >> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>> +lat_2=45.0
>>> > > +ellps=WGS84"
>>> > > >> projection(r)
>>> > > >>
>>> > > >> ## Assume that lon, lat are geographical coordinates (degrees
>>> decimal)
>>> > > >> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>> > > >>
>>> > > >> ## wrld_simpl is in geographical coordinates
>>> > > >> proj4string(wrld_simpl)
>>> > > >>
>>> > > >> ## Make figure in png format
>>> > > >> ## Of course plotting data with 2 different projections will give
>>> > > >> ## some distortions
>>> > > >> pdf("uro.png")
>>> > > >>
>>> > > >> plot(r)
>>> > > >> points(uro)
>>> > > >> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>> 'r'
>>> > > >>
>>> > > >> dev.off()
>>> > > >>
>>> > > >> extent(r)
>>> > > >> ## class       : Extent
>>> > > >> ## xmin        : -131.4368
>>> > > >> ## xmax        : -68.56323
>>> > > >> ## ymin        : 12.35567
>>> > > >> ## ymax        : 50.26619
>>> > > >>
>>> > > >> ## Reproject the raster to long-lat
>>> > > >> ## This doesn't work (collapsed domain)
>>> > > >> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>> > > +ellps=WGS84 +towgs84=0,0,0")
>>> > > >>
>>> > > >> ## Because
>>> > > >>> extent(rp)
>>> > > >> ## class       : Extent
>>> > > >> ## xmin        : -100.0015
>>> > > >> ## xmax        : -99.68557
>>> > > >> ## ymin        : 37.70658
>>> > > >> ## ymax        : 38.00046
>>> > > >>
>>> > > >> ## Save data in R binary format
>>> > > >> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>> > > >>
>>> > > >>
>>> > > >> Yours sincerely / Med venlig hilsen
>>> > > >>
>>> > > >> Frede Aakmann Tøgersen
>>> > > >> Specialist, M.Sc., Ph.D.
>>> > > >> Plant Performance & Modeling
>>> > > >>
>>> > > >> Technology & Service Solutions
>>> > > >> T +45 9730 5135
>>> > > >> M +45 2547 6050
>>> > > >> [hidden email]
>>> > > >> http://www.vestas.com
>>> > > >>
>>> > > >> Company reg. name: Vestas Wind Systems A/S
>>> > > >> This e-mail is subject to our e-mail disclaimer statement.
>>> > > >> Please refer to www.vestas.com/legal/notice
>>> > > >> If you have received this e-mail in error please contact the
>>> sender.
>>> > > >>
>>> > > >>
>>> > > >>
>>> > > >> -----Original Message-----
>>> > > >> From: R-sig-Geo [mailto:[hidden email]] On
>>> Behalf Of
>>> > > Agus Camacho
>>> > > >> Sent: 22. februar 2016 19:20
>>> > > >> To: [hidden email]
>>> > > >> Cc: r-sig-geo
>>> > > >> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to
>>> a
>>> > > reference system implicit in a .nc file
>>> > > >>
>>> > > >> Thanks Alex, but the locations still fall in the sea when i plot
>>> them
>>> > > using
>>> > > >> your recommended Solution. I looked at the sites you proposed and
>>> they
>>> > > have
>>> > > >> other values for lat_1, lat_0, etc..
>>> > > >>
>>> > > >> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>> > > >>
>>> > > >>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>> > > >>>> Dear all,
>>> > > >>>>
>>> > > >>>> Im trying to overlap these points:
>>> > > >>>>
>>> > > >>>
>>> > >
>>> >
>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>> > > >>>>
>>> > > >>>> and a wrld_simpl object:
>>> > > >>>> library(maptools)
>>> > > >>>> data(wrld_simpl)
>>> > > >>>>
>>> > > >>>> Over this raster layer
>>> > > >>>>
>>> > > >>>
>>> > >
>>> >
>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>> > > >>>>
>>> > > >>>> This rastr comes from a .nc file without a reference system. The
>>> > > author
>>> > > >>> of
>>> > > >>>> that .nc file gave me the following data about the .nc.
>>> > > >>>>
>>> > > >>>> The projection is *Lambert conformal conic* projection
>>> > > >>>> CEN_LAT = 38.0
>>> > > >>>> CEN_LON = -100.0
>>> > > >>>> TRUELAT1 = 25.
>>> > > >>>> TRUELAT2 = 45.
>>> > > >>>>
>>> > > >>>> However, despite i have gone through many sites in the
>>> internet, i
>>> > > cant
>>> > > >>>> figure it out:
>>> > > >>>>
>>> > > >>>> a) if that is all the data i need to set a reference system for
>>> my
>>> > > points
>>> > > >>>> and the wrld_simp object.
>>> > > >>>>
>>> > > >>>> b) how to change a typical CRS object with such data
>>> > > >>>>
>>> > > >>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>> > > >>>>
>>> > > >>>> Where do i enter the TRUELAT and CENLAT values?
>>> > > >>>> Are there any site that explains easily what the fields in the
>>> CRS
>>> > > mean
>>> > > >>> and
>>> > > >>>> how to change them?
>>> > > >>>>
>>> > > >>>> Thanks in advance.
>>> > > >>>>
>>> > > >>>
>>> > > >>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>> > > >>> https://trac.osgeo.org/proj/wiki/GenParms
>>> > > >>>
>>> > > >>> I believe:
>>> > > >>> +lat_0  = CEN_LAT   Latitude of origin
>>> > > >>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>> > > >>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>> > > >>> +lon_0  = CEN_LON   Central meridian
>>> > > >>>
>>> > > >>> proj strings are defined by the proj4 libary. It's website listed
>>> > above
>>> > > >>> and the associated mailing lists or gis stackexchange would be
>>> the
>>> > > >>> places to get help on it.
>>> > > >>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>> > > >>>
>>> > > >>> It often helps to browse similar projections on
>>> > > >>> http://spatialreference.org/
>>> > > >>> http://epsg.io/
>>> > > >>>
>>> > > >>> Enjoy,
>>> > > >>> Alex
>>> > > >>>
>>> > > >>>
>>> > > >>
>>> > > >>
>>> > > >
>>> > > > _______________________________________________
>>> > > > R-sig-Geo mailing list
>>> > > > [hidden email]
>>> > > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> > >
>>> > > --
>>> > > Roger Bivand
>>> > > Department of Economics, Norwegian School of Economics,
>>> > > Helleveien 30, N-5045 Bergen, Norway.
>>> > > voice: +47 55 95 93 55; fax +47 55 95 91 00
>>> > > e-mail: [hidden email]
>>> > > http://orcid.org/0000-0003-2392-6140
>>> > > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>> > > http://depsy.org/person/434412
>>> > > _______________________________________________
>>> > > R-sig-Geo mailing list
>>> > > [hidden email]
>>> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> >
>>> > --
>>> > Dr. Michael Sumner
>>> > Software and Database Engineer
>>> > Australian Antarctic Division
>>> > 203 Channel Highway
>>> > Kingston Tasmania 7050 Australia
>>> >
>>> >         [[alternative HTML version deleted]]
>>> >
>>> > _______________________________________________
>>> > R-sig-Geo mailing list
>>> > [hidden email]
>>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>> >
>>>
>>>
>>>
>>> --
>>> Agustín Camacho Guerrero.
>>> Doutor em Zoologia.
>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>> Biociências, USP.
>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>>
>
>
> --
> Agustín Camacho Guerrero.
> Doutor em Zoologia.
> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
> Biociências, USP.
> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
> São Paulo - SP, CEP: 05508-090, Brasil.
>



--
Agustín Camacho Guerrero.
Doutor em Zoologia.
Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
Biociências, USP.
Rua do Matão, trav. 14, nº 321, Cidade Universitária,
São Paulo - SP, CEP: 05508-090, Brasil.

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: adapting spatial points and wrld_smpl to a reference system implicit in a .nc file

Chris Reudenbach
Agus,

sorry for the addon but I think I have to provide a correction of the
corner coordinates (e.g. the extent values):

In the example that you have posted below I did calculate the extent
using the domain center coordinate and the WRF grid resolution in meter
and the number of rows and cols.

Since Dominik provides a link to the file description of the netcdf file
I think it is more accurate to reproject the  corner coordinates  as
given by the netcdf header variables (NC_GLOBAL#corner_lats,
NC_GLOBAL#corner_lons). Assuming that your variable "dctmx" (which I can
not identfy it in the nc file)  is of type "Mass" (staggered = M) the
correct corner coordinates are stored as the first 4 entrys of the dump
snippet below:

lats<-"
NC_GLOBAL#corner_lats={12.355667,50.26619,50.26619,12.355667,12.308136,50.18787,50.18787,12.308136,12.210785,50.403816,50.403816,12.210785,12.163345,50.325382,50.325382,12.163345}"

lons<-"
NC_GLOBAL#corner_lons={-131.43678,-151.29639,-48.703613,-68.563232,-131.5851,-151.51157,-48.488434,-68.414917,-131.38828,-151.41891,-48.581085,-68.611725,-131.53641,-151.63441,-48.36557,-68.463593}"

after cleaning and converting the strings you may calculate the corner
coordinates:


library(proj4)
## project mass corner coordinates to lambertian
llMass <- ptransform(cbind(clon[1],clat[1])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
ulMass <- ptransform(cbind(clon[2],clat[2])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
lrMass <- ptransform(cbind(clon[3],clat[3])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
urMass <- ptransform(cbind(clon[4],clat[4])/180*pi,'+proj=longlat
+datum=WGS84 +no_defs',proj)
wrfLccExtMass<-extent(ulMass[1],lrMass[1],llMass[2],ulMass[2])


According to this the correct extent for mass variables should be:

extent(wrfLccExtMass)
class       : Extent
xmin        : -3575343
xmax        : 3575342
ymin        : -2293058
ymax        : 2306330


hope this is correct now

cheers Chris

Am 24.02.2016 um 05:12 schrieb Agus Camacho:

> With the help of everybody i finally got this to work, so here is a script
> that does the job of reprojecting both, a raster layer obtained from a .nc
> and some locations in order to overplot them using plot.raster or mapview.
> Im using a combination of the advices of Dominik, Michael and Chris.
>
>
> require(ncdf4)
> require(raster)
>
>
> setwd("C:/~")
>
>
> r=raster("geo_em.d01.nc",
>                varname="dctmx")# days of ctmax events
>
> # Set extent and projections of rasters for plotting
> # chris gave me the orig data fom the nc file because i could not install
> gdal
> xmin= -3545999
> xmax= 3546000
> ymin = -2286000
> ymax=2286000
> pr<- "+proj=lcc +lat_1=25 +lat_2=45 +lat_0=38.000008 +lon_0=-100 +x_0=0
> +y_0=0"
>
> wrfLccExt<-extent(xmin,xmax,ymin,ymax)
>
> extent(r) <- extent(wrfLccExt)
> projection(r) <- pr
>
> # get and prepare  urosaurus locations
>
> x=read.csv("C:/Users/Agus/Dropbox/Functional Biogeography -
> Copy/scripts/class 4 maxent model/clean urosaurus records.csv")
> x=x[,1:3]
> colnames(x)
> coordinates(x)=~lon+lat
> proj4string(x)<- CRS("+proj=longlat +datum=WGS84")
> x=spTransform(x, pr)
>
>
> # get prepare and plot wrld_simpl
> require(maptools)
> require(sp)
> data(wrld_simpl)
> namer <- spTransform(subset(wrld_simpl, NAME %in% c("United States",
> "Canada", "Mexico")), prj)
>
> #plot with raster
> plot(r, cex.main=.7,legend=F)
> points(x)
> plot(namer, add = TRUE)
>
> # plot with mapview (cool!)
> m=mapview(r)
> u=mapview(x)
> m+u
>
>
> Thanks to all!
> Agus
>
> 2016-02-23 13:11 GMT-07:00 Agus Camacho <[hidden email]>:
>
>> Thanks for that Dominik,
>>
>> Giving that projection to either the locations, the raster layer generated
>> from the .nc file, or both, still did not work. I keep having locations
>> that should be on land falling far on the sea. Might this be a problem
>> derived from using raster with a file whose original grid distances are not
>> constant?
>>
>> Here is a link with the original file which has the original coordinate
>> data.
>>
>> https://www.dropbox.com/s/qpt5twtunhy3x3x/geo_em.d01.nc?dl=0
>>
>>
>> 2016-02-23 12:07 GMT-07:00 Dominik Schneider <
>> [hidden email]>:
>>
>>> This looks like WRF <http://www.wrf-model.org/index.php> data. I just
>>> dealt with this.
>>> The data is on a sphere as opposed to WGS84 so you need +ellps=sphere
>>> +a=6370000 +b=6370000 +units=m
>>>
>>> +proj=lcc which is usually what wrf is run with.
>>> The tricky part is:
>>> +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> because every WRF run is different (the WRF Preprocessing System
>>> optimizes the projection for the domain).
>>> and then there is probably no shift so you need(?) +x_0=0 +y_0=0
>>>
>>> This gives:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> But, wrf users like to give out lat and  long so you need to assign it:
>>> +proj=longlat +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> and then reproject the lat/long to lcc coordinates using this string:
>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0 +lon_0=-100.0
>>> +ellps=sphere +a=6370000 +b=6370000 +units=m +no_defs
>>>
>>> One word of caution, make sure you received the correct coordinates. Some
>>> variables are run cell center while some are run at cell edge. It looks
>>> like from your .nc file it was made by your collaborator so I assume they
>>> are right.
>>>
>>> That said, another word of caution, I found that the XLAT and XLONG
>>> variables from WRF output aren't very precise. There is a "geogrid" file
>>> from the preprocessing system that has the domain corners, resolution, nrow
>>> and ncol from which you can make a better grid using the native projection
>>> system (in my case it was a 4km grid). I suggest you try to get those.
>>>
>>> I hope this helps... I have to run but wanted to save people too much
>>> head scratching. I can get you running with more help tonight if you need.
>>> Dominik
>>>
>>>
>>> On Tue, Feb 23, 2016 at 11:27 AM, Agus Camacho <[hidden email]>
>>> wrote:
>>>
>>>> Thanks heaps to all for your effort. If I go to another GEOSTAT ill bring
>>>> more giant crab this time.
>>>>
>>>> The creator of the .nc file also looked at this webpage:
>>>> http://www.pkrc.net/wrf-lambert.html
>>>> It seemed like the right proj4 string might be this one:
>>>>
>>>> +proj=lcc +lat_1=25.0 +lat_2=45.0 +lat_0=38.0
>>>>      +lon_0=-100.0 +a=6370 +b=6370 +towgs84=0,0,0 +no_defs
>>>>
>>>> However this projection also does not allow me to adequately plot the
>>>> locations on the raster.
>>>>
>>>> Here is the .nc file. it contains several layers.
>>>>
>>>>
>>>> https://www.dropbox.com/s/yto3linsgom3zi7/results_us_future_output_none_0.nc?dl=0
>>>>
>>>>
>>>>
>>>> 2016-02-23 2:25 GMT-07:00 Michael Sumner <[hidden email]>:
>>>>
>>>>> On Tue, 23 Feb 2016 at 20:09 Roger Bivand <[hidden email]> wrote:
>>>>>
>>>>>> On Tue, 23 Feb 2016, Alex Mandel wrote:
>>>>>>
>>>>>>> I made an attempt at it too. Investigating the original data, I'm
>>>> not
>>>>>>> sure that the projection information supplied is correct for the
>>>> data
>>>>>>> linked. When I load up the data in a unprojected space, the
>>>> coordinates
>>>>>>> don't look at all similar to any Lambert projected data I have,
>>>> they
>>>>>>> actually look like Lat/Lon in some unprojected coordinate system,
>>>>>>> perhaps a different spheroid than expected.
>>>>>> Does anyone have a link to the original data? Is is possible that
>>>> this is
>>>>>> the General Oblique Transformation used by modellers - that is
>>>> something
>>>>>> that feels like longlat but is recentred and oblique? Example at the
>>>> very
>>>>>> end of my GEOSTAT talk last year (slides 81-83):
>>>>>>
>>>>>> http://geostat-course.org/system/files/geostat_talk_150817.pdf
>>>>>>
>>>>>> Roger
>>>>>>
>>>>>>
>>>>> For what it is worth, the General Oblique Transformation is not the
>>>> only
>>>>> example - it's very common for modellers to have a mesh that has the
>>>>> "mostly-properties" of a projection, but is not actually describable
>>>> with
>>>>> standard transform + affine parameters. The main cases that I've seen
>>>> are
>>>>> polar stereographic, equal area or oblique Mercator. Often they really
>>>> are
>>>>> simple transforms and you can reconstruct without loss, but it's not
>>>>> usually possible to tell without exploration. It's an interesting
>>>>> dis-connect to see code that builds a mesh with certain properties,
>>>> then
>>>>> only stores longitudes and latitudes - when it could be done with
>>>> standard
>>>>> tools and be stored and used much more efficiently.
>>>>>
>>>>> (I've seen Lambert Conformal Conic and Lambert Azimuthal Equal Area
>>>>> terminology conflated in this context too. )
>>>>>
>>>>> I'm also interested to explore the original data.
>>>>>
>>>>> Cheers, Mike.
>>>>>
>>>>>
>>>>>
>>>>>>> -Alex
>>>>>>>
>>>>>>> On 02/22/2016 10:17 PM, Frede Aakmann Tøgersen wrote:
>>>>>>>> Hi
>>>>>>>>
>>>>>>>> I tried to make it work but I had to give up. I wanted to
>>>> reproject
>>>>> the
>>>>>> Lamberth conformal conic coordinates to long-lat but it didn't work.
>>>>>>>> Perhaps someone can see what I did wrong. Here is what I did
>>>> (data in
>>>>> R
>>>>>> binary format and figure in png format both attached):
>>>>>>>> library(raster)
>>>>>>>> library(maptools)
>>>>>>>> data(wrld_simpl)
>>>>>>>>
>>>>>>>> r <- raster("raster.grd")
>>>>>>>> projection(r)
>>>>>>>> ## > NA
>>>>>>>>
>>>>>>>> uro <- read.table("clean urosaurus records.csv", h = TRUE, sep =
>>>> ",")
>>>>>>>> coordinates(uro) <- ~lon+lat
>>>>>>>>
>>>>>>>> ## Set projections for the 3 data sets
>>>>>>>>
>>>>>>>> ## Lamberth's confocal conic projection with given parameters
>>>>>>>> crs(r) <- "+proj=lcc +lat_0=38.0 +lon_0=-100 +lat_1=25.0
>>>> +lat_2=45.0
>>>>>> +ellps=WGS84"
>>>>>>>> projection(r)
>>>>>>>>
>>>>>>>> ## Assume that lon, lat are geographical coordinates (degrees
>>>> decimal)
>>>>>>>> proj4string(uro) <- CRS("+proj=longlat +datum=WGS84")
>>>>>>>>
>>>>>>>> ## wrld_simpl is in geographical coordinates
>>>>>>>> proj4string(wrld_simpl)
>>>>>>>>
>>>>>>>> ## Make figure in png format
>>>>>>>> ## Of course plotting data with 2 different projections will give
>>>>>>>> ## some distortions
>>>>>>>> pdf("uro.png")
>>>>>>>>
>>>>>>>> plot(r)
>>>>>>>> points(uro)
>>>>>>>> plot(wrld_simpl, add = TRUE) # World will be clipped to extent of
>>>> 'r'
>>>>>>>> dev.off()
>>>>>>>>
>>>>>>>> extent(r)
>>>>>>>> ## class       : Extent
>>>>>>>> ## xmin        : -131.4368
>>>>>>>> ## xmax        : -68.56323
>>>>>>>> ## ymin        : 12.35567
>>>>>>>> ## ymax        : 50.26619
>>>>>>>>
>>>>>>>> ## Reproject the raster to long-lat
>>>>>>>> ## This doesn't work (collapsed domain)
>>>>>>>> rp <- projectRaster(r, crs = "+proj=longlat +datum=WGS84 +no_defs
>>>>>> +ellps=WGS84 +towgs84=0,0,0")
>>>>>>>> ## Because
>>>>>>>>> extent(rp)
>>>>>>>> ## class       : Extent
>>>>>>>> ## xmin        : -100.0015
>>>>>>>> ## xmax        : -99.68557
>>>>>>>> ## ymin        : 37.70658
>>>>>>>> ## ymax        : 38.00046
>>>>>>>>
>>>>>>>> ## Save data in R binary format
>>>>>>>> save(list = c("r", "uro", "wrld_simpl"), file = "uro.RData")
>>>>>>>>
>>>>>>>>
>>>>>>>> Yours sincerely / Med venlig hilsen
>>>>>>>>
>>>>>>>> Frede Aakmann Tøgersen
>>>>>>>> Specialist, M.Sc., Ph.D.
>>>>>>>> Plant Performance & Modeling
>>>>>>>>
>>>>>>>> Technology & Service Solutions
>>>>>>>> T +45 9730 5135
>>>>>>>> M +45 2547 6050
>>>>>>>> [hidden email]
>>>>>>>> http://www.vestas.com
>>>>>>>>
>>>>>>>> Company reg. name: Vestas Wind Systems A/S
>>>>>>>> This e-mail is subject to our e-mail disclaimer statement.
>>>>>>>> Please refer to www.vestas.com/legal/notice
>>>>>>>> If you have received this e-mail in error please contact the
>>>> sender.
>>>>>>>>
>>>>>>>>
>>>>>>>> -----Original Message-----
>>>>>>>> From: R-sig-Geo [mailto:[hidden email]] On
>>>> Behalf Of
>>>>>> Agus Camacho
>>>>>>>> Sent: 22. februar 2016 19:20
>>>>>>>> To: [hidden email]
>>>>>>>> Cc: r-sig-geo
>>>>>>>> Subject: Re: [R-sig-Geo] adapting spatial points and wrld_smpl to
>>>> a
>>>>>> reference system implicit in a .nc file
>>>>>>>> Thanks Alex, but the locations still fall in the sea when i plot
>>>> them
>>>>>> using
>>>>>>>> your recommended Solution. I looked at the sites you proposed and
>>>> they
>>>>>> have
>>>>>>>> other values for lat_1, lat_0, etc..
>>>>>>>>
>>>>>>>> 2016-02-22 11:04 GMT-07:00 Alex M <[hidden email]>:
>>>>>>>>
>>>>>>>>> On 02/22/2016 09:50 AM, Agus Camacho wrote:
>>>>>>>>>> Dear all,
>>>>>>>>>>
>>>>>>>>>> Im trying to overlap these points:
>>>>>>>>>>
>>>> https://www.dropbox.com/s/awdclg4cvsdngej/clean%20urosaurus%20records.csv?dl=0
>>>>>>>>>> and a wrld_simpl object:
>>>>>>>>>> library(maptools)
>>>>>>>>>> data(wrld_simpl)
>>>>>>>>>>
>>>>>>>>>> Over this raster layer
>>>>>>>>>>
>>>> https://www.dropbox.com/sh/qcw174tgogpnz7s/AAByDc3TeyFe3W4nEqTFix6Oa?dl=0
>>>>>>>>>> This rastr comes from a .nc file without a reference system. The
>>>>>> author
>>>>>>>>> of
>>>>>>>>>> that .nc file gave me the following data about the .nc.
>>>>>>>>>>
>>>>>>>>>> The projection is *Lambert conformal conic* projection
>>>>>>>>>> CEN_LAT = 38.0
>>>>>>>>>> CEN_LON = -100.0
>>>>>>>>>> TRUELAT1 = 25.
>>>>>>>>>> TRUELAT2 = 45.
>>>>>>>>>>
>>>>>>>>>> However, despite i have gone through many sites in the
>>>> internet, i
>>>>>> cant
>>>>>>>>>> figure it out:
>>>>>>>>>>
>>>>>>>>>> a) if that is all the data i need to set a reference system for
>>>> my
>>>>>> points
>>>>>>>>>> and the wrld_simp object.
>>>>>>>>>>
>>>>>>>>>> b) how to change a typical CRS object with such data
>>>>>>>>>>
>>>>>>>>>> Ex.CRS ("+proj=lcc+lat_0=38.0+lon0_2=-100+ellps=WGS84")
>>>>>>>>>>
>>>>>>>>>> Where do i enter the TRUELAT and CENLAT values?
>>>>>>>>>> Are there any site that explains easily what the fields in the
>>>> CRS
>>>>>> mean
>>>>>>>>> and
>>>>>>>>>> how to change them?
>>>>>>>>>>
>>>>>>>>>> Thanks in advance.
>>>>>>>>>>
>>>>>>>>> https://github.com/OSGeo/proj.4/wiki/GenParms
>>>>>>>>> https://trac.osgeo.org/proj/wiki/GenParms
>>>>>>>>>
>>>>>>>>> I believe:
>>>>>>>>> +lat_0  = CEN_LAT   Latitude of origin
>>>>>>>>> +lat_1  = TRUELAT1   Latitude of first standard parallel
>>>>>>>>> +lat_2  = TRUELAT2   Latitude of second standard parallel
>>>>>>>>> +lon_0  = CEN_LON   Central meridian
>>>>>>>>>
>>>>>>>>> proj strings are defined by the proj4 libary. It's website listed
>>>>> above
>>>>>>>>> and the associated mailing lists or gis stackexchange would be
>>>> the
>>>>>>>>> places to get help on it.
>>>>>>>>> https://lists.osgeo.org/mailman/listinfo/metacrs
>>>>>>>>>
>>>>>>>>> It often helps to browse similar projections on
>>>>>>>>> http://spatialreference.org/
>>>>>>>>> http://epsg.io/
>>>>>>>>>
>>>>>>>>> Enjoy,
>>>>>>>>> Alex
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> R-sig-Geo mailing list
>>>>>>> [hidden email]
>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>> --
>>>>>> Roger Bivand
>>>>>> Department of Economics, Norwegian School of Economics,
>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>> voice: +47 55 95 93 55; fax +47 55 95 91 00
>>>>>> e-mail: [hidden email]
>>>>>> http://orcid.org/0000-0003-2392-6140
>>>>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>>>> http://depsy.org/person/434412
>>>>>> _______________________________________________
>>>>>> R-sig-Geo mailing list
>>>>>> [hidden email]
>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>> --
>>>>> Dr. Michael Sumner
>>>>> Software and Database Engineer
>>>>> Australian Antarctic Division
>>>>> 203 Channel Highway
>>>>> Kingston Tasmania 7050 Australia
>>>>>
>>>>>          [[alternative HTML version deleted]]
>>>>>
>>>>> _______________________________________________
>>>>> R-sig-Geo mailing list
>>>>> [hidden email]
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>
>>>>
>>>>
>>>> --
>>>> Agustín Camacho Guerrero.
>>>> Doutor em Zoologia.
>>>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>>>> Biociências, USP.
>>>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>>>> São Paulo - SP, CEP: 05508-090, Brasil.
>>>>
>>>>          [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>> --
>> Agustín Camacho Guerrero.
>> Doutor em Zoologia.
>> Laboratório de Herpetologia, Departamento de Zoologia, Instituto de
>> Biociências, USP.
>> Rua do Matão, trav. 14, nº 321, Cidade Universitária,
>> São Paulo - SP, CEP: 05508-090, Brasil.
>>
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
12