texts: Add information from a raster nc-file to a polygon shapefile

texts:
DATA
Download from this website the following ZIP file: NUTS_2006_60M_SH.zip and save all the NUTS_RG_60M_2006 files. In what follows I refer to this as "NUTS3 regions".
Download from this website the following nc-file: cru_ts3.23.1991.2000.pet.dat.nc.gz This is a raster with data all over the world.
Desired output
The nc-file is a raster of the world with values in different colors.

The goal is now to translate these data to a NUTS3 polygon regions and save it back as a nuts3-shapefile with the data added to the polyons.

Software
I can work with R, but since I couldn't make it work in R, I switched to ArcGIS. If the data are added to the NUTS3 shapefile, I can continue working with it in R.
Any solution in R or ArcGIS is welcome, as long as the output is 1 polygon shapefile with all the values added.
Problem
I can open the nc-file in ArcGIS (toolbox, multi-dimension tools, make netCDF Raster Layer). I can also overlap the nuts3-file with the raster. But I can't merge the data together in one shapefile (by taking the average of the underlying rasters that overlap with one polygon) and drop all the data that I don't need.
I found some examples on how to do this on Google, but only one of them explained how to start from an nc-file. As far as I understood, I first have to right click on the data-layer to export the data as a raster (as it is now only an nc-file and not yet a shapefile). However, this is not working since I get the following error: "Output file format with specified pixel type or number of bands or colormap is not supported."


texts:
I am studying bird movements and have tracked their locations, ie from flower 1 to flower 2 to flower 3 etc. Each point is mapped out. I am trying to work out the relative direction for each movement of the bird. in the picture below the angles are A1 and A2. so A1 is the relative direction of the bird travelling from N2 to N3 based on the original direction of N1 to N2. I have over 1000 movements which are copletely random with up to 50 flowers per map. is there a way to run the movements with the outcome relative angle?

At this point the bird movements are stored in a spreadsheet. I am unsure of how to incorporate them for this purpose.

P.S. im relatively new to gis so knowledge is limited


texts: QGIS: dissolve with statistics (summarize) for several attribute fields in graphical modeler

texts:
I'm trying to dissolve a shapefile by one attribute field "HYBAS_ID" and I need the other fields to be summarized simultaneously. I've worked a lot with ArcGis and the dissolve tool there had the possibility to do so, but I need the model to be in QGIS. The QGIS dissolve tool does not have the option at all and the gdal dissolve tool does'nt let me choose if I want to summarize, it just has the options to do all the operations at once and only for one field (at least as im aware, maybe the additional creation options can do this, but I dont know how).

I'm aware of the "dissolve with stats" plugin, but I have to add the operation as a tool into a model. As said, I have to dissolve the field "HYBAS_ID" and need to summarzie the fields "AgriET", "NabesET", "WR_t", "DWD", "Wliv" and "Wind".

Is there a possibility to include this into a model? Unfortunately I'm not frequent in any coding, so if it's the only possibility, i'd be grateful if you could explain it a little bit.


texts: Adding other data to points in ArcMap document?

texts:
I'm a user (not producer) of ArcGIS documents in my work and I need some help asking experts with the following:

I have some Microsoft Word documents and some pictures that describe characteristics of buildings that I have represented in ArcMap 10.1 document as a points in a shape.
I need to have access to the data and pictures when I select any of these points (buildings) clicking on ArcMap project and read this info when needed.

Can I add this both type of data in ArcMap or should I create a separated database and link ArcMap and this DB in any way?

Should I create a PDF file joining MS Word and images to have only one file and access them?


texts:
I am trying to install Lizmap on a server that runs on a portal:
portal.domain.com

The server itself runs on the portal with a prefix(?) added to the url  that the server is unaware of to my understanding.
portal.domain.com/prefix

So when I install Lizmap on the server I can access the landing page (no css) but any generated links will point to portal.domain.com/lizmap.
I did not find any solution to this problem.
How can I include this prefix in the installation of lizmap?


texts: Find AB_CSRS.gsb for NAD_1983_TO_WGS_1984_8 Transformation

texts:
I am trying to change my WGS84 data to NAD83 for Alberta data. My research tells me that I need to use the NAD_1983_TO_WGS_1984_8 Datum Transformation which is in the AB_CSRS.gsb. I can not find this gsb file for download on the internet. I did manage to find and obtain the NTv2_0.GSB but that seems to only be for NAD27 To NAD83.

Does anyone know where I can get the AB_CSRS.GSB file?

Thanks


texts:
I want compile GDAL library for use in C# program. I use Microsoft Visual Studio 2013. First, I compiled GDAL C++. As a result I have DLLs and libraries (No problems) . Next, I tried to compile with instruction, because I want use GDAL in C# program. But, I don't understand what should I do in this step.

  Creating the SWIG interface code
  For creating the interface execute the following command.

  nmake /f makefile.vc interface

  Previous Step make's error when regenerates the *PINVOKE.cs files ,so
  skip to next step ,it will work

I tried to skip this step, but always had the same exception (shown below). As result I don't have DLLs to use in C#...

    ogr\OsrPINVOKE.cs(192,10): error CS0111: Type 'OSGeo.OGR.OsrPINVOKE' already
        defines a member called 'OsrPINVOKE' with the same parameter types
ogr\OsrPINVOKE.cs(188,10): (Location of symbol related to previous error)
NMAKE : fatal error U1077: '"C:\Program Files (x86)\MSBuild\12.0\bin\csc.EXE"' :
 return code '0x1'
Stop.

What should I do to compile GDAL to C#?


texts: Assigning field value based on location without join

texts:
I'm interested in updating a field in geodatabase feature classes based on their location in relation to another feature class in the geodatabase, but without having to perform a join (and therefore creating another feature class).

I've got a polygon feature class ("boundaries") and then multiple point and line features that fall within those boundaries ("fieldData," "streams," "origins"). All of the features, once loaded to the database, have a common field labeled "siteCode" that I would like to populate programmatically rather than manually doing a Select by Location and calculating the fields (or creating extra files like performing a join). Better yet, I would like this to occur as the point and line files are loaded to the geodatabase, but that's not necessary. None of the files have this information stored in their attributes prior to being loaded to the database, but the siteCode is the first two letters in the name of each source shapefile.

As it stands now, I select a feature in the "boundaries" feature class, Select by Location on each of the other files one at a time, and populate the fields using Calculate Field.

I am open to any suggestions for workflow improvement here! My experience is limited with Python but I am interested in learning. If people respond with Python suggestions I would just ask that they include as much information as possible! I feel like this must have been answered a million times but I'm seeing a lot about using a join and not many ways to avoid it.

Also: There is this posting that appears to be accomplishing something similar but I'm too much of a novice to understand exactly what it's doing and if it will do what I need. If someone can confirm that some variation on that code will work for me I will dive in and figure it out!

    #Make a feature layer for the neighourhoods
arcpy.MakeFeatureLayer_management(neighbourhood, "neighbourhoodLayer")
arcpy.MakeFeatureLayer_management(parcel, "parcelLayerTest1")
#Setup cursor to update new NBHD field
with arcpy.da.UpdateCursor(parcel, (parcelIDField,nbhdField)) as parRows:
    for par in parRows:
        # Create a query string for the current patrol zone
        #parcelIDstring = par[0]
        #queryString = '"' + str(parcelIDField) + '" =' + "'" + str(parcelIDstring) + "'"

        #Select parcels based on location contained within neighbourhood
        arcpy.SelectLayerByLocation_management("parcelLayerTest1", "HAVE_THEIR_CENTER_IN", neighbourhood)

        nbhdTable = arcpy.SearchCursor(neighbourhood)
        for nbhd in nbhdTable:
            nbhdValue = nbhd.getValue("Neighbourh")
            if nbhdValue:
                par[1] = nbhdValue
                parRows.updateRow(par)

Thanks in advance for any assistance.


texts: Creating geodatabase (and st_geometry) in Oracle Express (XE) for home/personal use?

texts:
I'd like to set up a database environment at home that mimics what I have in the office (Oracle enterprise geodatabase/ESRI's st_geometry). I want to do this for spatial database training purposes in an Oracle/ESRI environment.

I plan to set up ArcGIS Desktop Advanced (Personal Use license) as the client and Oracle XE 11g R2 (64-bit) as the database. I'll create a geodatabase in the Oracle XE database, and add a st_geometry user-defined type.

From what I've read, it looks like this is possible, in theory . The Oracle database requirements for ArcGIS 10.3.x document says that Oracle 11g R2 (64-bit) 11.2.0.3 is supported by ArcGIS desktop. However, it doesn't explicitly say that the express edition (XE) of Oracle is supported. Also, it has occurred to me that I've never heard of anyone else doing this.

Since I will be purchasing new hardware and software solely for this purpose, I'd like to know ahead of time -- will this work? Can I create a geodatabase (with st_geometry) using ArcGIS Desktop Advanced (Personal Use license) in an Oracle XE 11g R2 (64-bit) database?


texts: No height information in GoogleEarth paths?


texts: different browsers referencing different css files, arcgis javascript cmv

texts: When I open a page in chrome, and click on a feature the popup style is coming from one css file, and when I do the same thing in firefox, the popup looks completely different and is referencing a different css. I am very new to developing in javascript and don't know very much yet. Why would this occur? I started out with the CMV template https://github.com/cmv/cmv-app which has a main.css file that is downloaded with the source code, and the index page also references esri.css on the web. I found the popup style was coming from the esri.css and wanted to make some changes, so save it locally and changed my html file to point to the local version. Once I did this, I started getting the funky popup in firefox. How can I correct this? I notice there are some settings for the .esriPopup in the main.css, but when I tried to modify these nothing changed.


texts: Making Python Add-In using tools built in ModelBuilder?

texts:
I have made a set of tools in a toolbox using ArcGIS 10.2 using the model builder.

I would like to distribute them as a tool bar, and so created a ArcGIS Python Add-In I was hoping that I could export the tools as python scripts and copy the code into the Add-In. However, it does not seem to be as simple as that.

Is it possible with Python Add-Ins to make a window for entering in parameters and run a script in the way tools/toolboxes work?

Update

I found a useful blog post.

This is my code

import arcpy
import pythonaddins
import os.path

class multiplebenefits(object):
    """Implementation for ArcGISAddin_addin.multiplebenefits (Tool)"""
    def __init__(self):
        self.enabled = True
        self.checked = False

    def onClick(self):
       # name of toolbox without tbx extension
       toolboxName = "mytoolbox"
       # name of tool to be executed
       toolName = "mymodel"
       # create string with path to toolbox
       toolboxPath = os.path.join(os.path.dirname(__file__), toolboxName + ".tbx")
       # call geoprocessing tool
       pythonaddins.GPToolDialog(toolboxPath, toolName)

But I get an error message:

Failed to open tool mymodel (mytoolbox.tbx)

TypeError: GPToolDialog() takes at most 1 argument (2 given)


texts: Way to analyze size of mbtiles file or pg database to find largest objects to optimize final size

texts:
Is there a way to figure out what OSM data is taking up the most space in my final mbtiles file?

For example, Connecticut is 54MB after I have made some modifications to the scripts in openmaptiles. This is larger than I want, and I want to know what is taking up so much space. Is there any such method of doing this?


texts: Vegetation Index for dry and sparse vegetated areas?

texts:
I found many different index for the estimation of vegetation cover, but many of them struggle in sparse areas with low vegetation cover, as the soil signals are way stronger than the vegetation signals.

Is there a vegetation index suited for dry and sparse vegetated areas?


texts: runGdal returns "gdalinfo failed - unable to open" [...] ".hdf" file

texts:
I'm trying to use the function runGdalin R on a Mac and it seems that gdalinfo is not able to open the file.

In the terminal, when I entered
gdalinfo /MODIS/MCD12Q1.006/2009.01.01/MCD12Q1.A2009001.h10v05.006.2018054191941.hdf

I got the following:

Driver: HDF4/Hierarchical Data Format Release 4
Files: /MODIS/MCD12Q1.006/2009.01.01/MCD12Q1.A2009001.h10v05.006.2018054191941.hdf
Size is 512, 512
Coordinate System is `'
Metadata:
[...]

So GDAL seems to be working (version of gdalinfo: GDAL 2.1.3, released 2017/20/01). But in R, it doesn't seem to be recognized.

In R, I set up the path to GDAL like this:

MODISoptions(MODISserverOrder = c("LAADS", "LPDAAC"),
             localArcPath = "my.path",
             outDirPath = "my.path",
             gdalPath = "/Library/Frameworks/GDAL.framework/Versions/Current/Programs/",
             quiet = FALSE)

Running this in terminal showed that HDF4 is supported:

gdalinfo --formats
  [...]
  HDF4Image -raster- (rw+): HDF4 Dataset
  [...]

In R, I typed system("gdalinfo --version") and got back GDAL 2.4.1, released 2019/03/15

system("gdalinfo                                                              /path.to.my.hdf")

ERROR 4: `... .hdf' not recognized as a supported file format.
gdalinfo failed - unable to open '... .hdf'.

system("/Library/Frameworks/GDAL.framework/Versions/Current/Programs/gdalinfo /path.to.my.hdf")

This last line one worked (with the full path to gdal)! But how to make it work in runGdal?

Is there a way to make the runGdal function work in R on a Mac?


texts: Render map data faster in Leaflet using MapTiler

texts:
I have a Vuejs app containing a Leaflet map that displays GeoJSON files. Some of these files, such as the file containing polygons for all zip codes in a US state are very large and load too slowly. Our client code pulls this file from an external github repository. How can we achieve this more quickly? I am open to using a paid third party software, such as MapTiler. I see it is possible to create custom raster map tiles with MapTiler, but how can I create a custom map containing zip code polygons? We would like to color each zip code to match a statistic such as the average age in that zip code. We can do this already by fetching the age statistics from our database and mapping them to the zip code in the GeoJSON. How can we still achieve this chloropleth end result, but faster?

Update: I learned about the Leaflete.vectorGrid npm package at  @TomazicM's suggestion. I am able to integrate the example into our app like this...

  var vectorGrid = L.vectorGrid
        .slicer(euCountries.euCountries, {
          rendererFactory: L.svg.tile,
          vectorTileLayerStyles: {
            sliced: function(properties, zoom) {
              var p = properties.mapcolor7 % 5;
              return {
                fillColor:
                  p === 0
                    ? "#800026"
                    : p === 1
                    ? "#E31A1C"
                    : p === 2
                    ? "#FEB24C"
                    : p === 3
                    ? "#B2FE4C"
                    : "#FFEDA0",
                // fillOpacity: 0.5,
                fillOpacity: 1,
                stroke: true,
                fill: true,
                color: "black",
                //opacity: 0.2,
                weight: 0
              };
            }
          },
          interactive: true,
          getFeatureId: function(f) {
            return f.properties.wb_a3;
          }
        })
        .on("mouseover", function(e) {
          var properties = e.layer.properties;
          L.popup()
            .setContent(properties.name || properties.type)
            .setLatLng(e.latlng)
            .openOn(myMap);

          clearHighlight();
          highlight = properties.wb_a3;

          var p = properties.mapcolor7 % 5;
          var style = {
            fillColor:
              p === 0
                ? "#800026"
                : p === 1
                ? "#E31A1C"
                : p === 2
                ? "#FEB24C"
                : p === 3
                ? "#B2FE4C"
                : "#FFEDA0",
            fillOpacity: 0.5,
            stroke: true,
            fill: true,
            color: "red",
            opacity: 1,
            weight: 2
          };

          vectorGrid.setFeatureStyle(properties.wb_a3, style);
        })
        .addTo(myMap);

However, when I try to add a different GeoJSON file (containing polygons of Illinois zip codes), I get the error:

`Leaflet.VectorGrid.js?45d4:1931 Uncaught (in promise) DOMException: Failed to execute 'postMessage' on 'Worker': function (f) {
      return f;
    } could not be cloned.
    at NewClass.initialize at function L.vectorGrid.slicer...

      var IllinoisZips = L.vectorGrid
    .slicer(ILZipCodes, {
      rendererFactory: L.svg.tile,
      vectorTileLayerStyles: {
        sliced: function(properties, zoom) {
          // var p = properties;
          return {

          };
        }
      }
    })

    .addTo(myMap);

And this is how I'm getting the ILZipCodes GeoJSON layer...

   var ILZipCodes = new L.geoJson.ajax(
        "https://raw.githubusercontent.com/OpenDataDE/State-zip-code-GeoJSON/master/il_illinois_zip_codes_geo.min.json"
      );

Why might this be?


texts: QGIS 3.16 - Print layout grid coordinates in DMS?

texts:
I am trying to create a grid with coordinates printed in degrees, minutes, and seconds in print layout of QGIS 3.16. The project is defined to EPSG:32652 WGS 84/UTM Zone 52N.
When I go to Grid > Draw Coordinates > Format, I do not get the usual DMS option. There are only three options:
The following is what I get when I change to CRS to the default of EPSG:4326 - WGS 84, but when I choose the option, the coordinates still do not appear on the map:

What should I be doing to bring the DMS option using my project's CRS? Why am I only getting decimal, decimal with suffix, custom?


texts: Using rdgal on R- error- No GDAL installation found

texts:
I trying to map burnt area using hdf data. For this i have all the neccessary packages, including the rdgal pacakge.
However when i type in my burnt area code, i get the response
I have downloaded rgdal as it is the R version of gdal, and apparently should work fine, however i cant seem to get rid of this error.
I have tried downloaded the gdal_setInstallation but again get errors saying it cannot find the function.
Has anyone experienced this?


texts: Methods of cloud shadow detection for Landsat data?

texts:
Detecting clouds doesn't seem to be much of a problem, rather the challenge is in selecting cloud-shadow classes that don't also include coastal areas, agriculture, etc.  For a project I'm working on, I've been using two methods:

Fmask (https://code.google.com/p/fmask/) is powerful and fully automated, plus it works for TM, ETM+, and OLI/TIRS data.  It is great at catching harder to detect haze and the majority of other clouds/shadows.  The problem I have with it however is that it tends to miss certain cloud types, namely smaller low altitude clouds and their shadows.  I played around with the input parameters and found the best success using default values.

My other method is a semi-automated classification.  I create NDSI, NDVI, and NDWI for a particular scene, stack those as RGB, and perform unsupervised classification in ERDAS Imagine 2013 (k-means, 10 iterations, 0.95 convergence threshold, and usually 256 classes).  After working on several scenes, I preferred to use the stacked indices over classifying individual bands because it took less time and did a slightly better job at separating clouds/shadows from other classes.

In the end, I had fairly good success with detecting and removing clouds by combining these two methods, but ended up losing good data in my classification of cloud-shadows in particular.  Has anyone had success in removing shadow areas that they'd like to share, or perhaps some insight into improving my classification process?


texts: Styling SLD Files for Geoserver XML Editors and Notepad++

texts:
I am currently using Notebook++ to edit my GeoServer Styles (Yep, I like to do them the hard way). This is a Windows Environment.

I also use XML Copy and XML Notepad 2007 Editor but so far I have been unable to to get a DTD/Schema to work with it for SLD and more specifically the GeoServer Extensions to the SLD Specification.

Questions:

Does anyone have a Style for SLD for Notepad++ or know where I can get one?
Are there any better XML Editors out there that will use the SLD DTD?
Where can I get an XLST file for Geoserver SLD and extensions or vanilla SLDs at the very least


texts: Create a new GRASS database in R with CRS projection

texts:
How to create a new GRASS database in R with CRS projection? I am able to create a new one using the following code but it is not projected

require(raster)
require(spgrass6)
initGRASS("/Applications/GRASS-6.4.app/Contents/MacOS", gisDbase="/Users/junf/Documents",
          location='temp',mapset='PERMANENT', override=TRUE)
# gisdbase    /Users/junf/Documents
# location    temp
# mapset      PERMANENT
# rows        1
# columns     1
# north       1
# south       0
# west        0
# east        1
# nsres       1
# ewres       1
# projection  NA

Now I tried using a projected GeoTiff using the following code which is better but the projection is still NA:

r <- as(raster('/Users/junf/Documents/WatershedModel/UpperPulangi/upw_srtm30.tif'), 'SpatialGrid')
initGRASS("/Applications/GRASS-6.4.app/Contents/MacOS", gisDbase="/Users/junf/Documents", location='temp',mapset='PERMANENT', override=TRUE,SG=r)
# gisdbase    /Users/junf/Documents
# location    temp
# mapset      PERMANENT
# rows        3317
# columns     2774
# north       950826
# south       849428.7
# west        690839.6
# east        775638
# nsres       30.56898
# ewres       30.569
# projection  NA

Can someone please help me create a new GRASS database with WGS 84 (EPSG: 4326) projection
or a projection that matches an existing georeferenced image for import in GRASS?


texts:
I've been trying to get v.krige to run in the Grass7 command line, however, when I do so I get the following error:

Error : .onAttach failed in attachNamespace() for 'rgrass7', details:
call: if (class(t0) != "try-error" && is.character(t0) && nchar(t0) >
error: missing value where TRUE/FALSE needed
Checking rgeos availability: TRUE
ERROR: R package(s) rgrass7 missing. Install it/them and re-run v.krige.

However, when I check in R what libraries are installed I clearly show that rgrass7 is installed.

As a note I get a very similar error if I open R inside Grass7 using the following commands:

> library(rgrass7)
Loading required package: sp
Loading required package: XML
Error : .onAttach failed in attachNamespace() for 'rgrass7', details:
call: if (class(t0) != "try-error" && is.character(t0) && nchar(t0) >
error: missing value where TRUE/FALSE needed
Error: package or namespace load failed for 'rgrass7'

So it appears R inside of Grass7 isn't able to load this module, but if I run the same command inside the RGui, I get this:

> library(rgrass7)
Loading required package: sp
Loading required package: XML
GRASS GIS interface loaded with GRASS version: (GRASS not running)

Anyone know what might be going on and how I can get v.krige to work?


texts: raster sytling multiband rendering 0=nodata

texts:
Is there any way to set the pixel value 0 to nodata of a raster stack in the multiband rendering type in QGIS?

I remember that it is possible to check a box in arcgis for that, so maybe i just oversee it in QGIS?

It is pretty annoying to convert all 0 values to nodata just for the style.

The min value for the rgb bands should start at the real min and ignore the 0 values


texts: How to know the GSD (grid size, resolution) in meters of an existing raster when working in Geographic coordinates?

texts:
I would need to know the GSD or grid size (resolution), in meters, of a raster created when working in a Geographic coordinate system in ArcGIS. Is is possible?

Even better, can I set the Cell size of the raster in meters (not in degree units) when I perform the raster creation if I am working in a Geographic coordinate system?


texts: Display GPX file on a GPSmap 62s GPS device

texts:
I have to do a regular survey and for that I have produced a grid on QGIS 3.0 (Girona)
After producing the (vector) grid I save it as a GPX file (as proposed)
I Saw that the GPX_USE_EXTENSIONS should be activated
I also activate FORCE_GPX_ROUTE and FORCE_GPX-TRACK

Then the produced .gpx file is copied into /Garmin/GPX folder of the GPAmap 62 s device

From there I can't see my grid displayed on the GPS screen

Have I missed something ?


texts: Apparently MOBAC is not generating MBTiles correctly

texts:
I am trying to generate MBTiles with MOBAC, to be used in an offline map Android application (trying to use OSMDroid).

Now I generate the atlas, but something goes wrong because I get the error

  android.database.sqlite.SQLiteException: no such table: images (code 1): , while compiling: SELECT tile_data FROM images LIMIT 0,1

This is my MOBAC setup:


texts:
I am trying to create a flood modelling system. The system consists of two parts mainly a river flooding simulation and rainfall simulation.  To do this I am running the latest version of GRASS 6.4. Also I am using the SRTM obtained from CGIAR as my DEM.

I would like to ask what are the most appropriate commands to run a simulation. Ex. Terraflow

Also since we are modelling flooding based on a river, is it advisable to use HEC-RAS? It was created by the same folks who made GRASS. :)

I am more of a programmer rather than an engineer. It is my first venture into simulation systems as part of my degree.

As of now Here is our plan for the rainfall (pseudocode):

int i = 0, j = 0;
int passes, numberOfmesh;
double rainAmount;
double store1[][];
double store2[][];

numberOfmesh = getTotalMeshCells();
passes = computePasses(numberOfmesh());
rainAmount = getRainAmount();

/* rainAmount corresponds to the input of the user in the textfield.
store1 contains the elevation and the number of times an amount of water
has been passed to it (per pixel)
store2 contains the elevation and the total accumulated water (per pixel)
Both store1 and store2 have 'n' number of rows and 2 columns;
ASSUME THAT THE ARRAY HAS ALREADY BEEN POPULATED */

while(i != passes)
{
    while(j != numberOfmesh)
    {
        direct = getDirection(store[j][0]);
        /* direct corresponds to the direction of a mesh cell;
        I don't know what id the datatype for direct; */

        store2[j][6] += rainAmount;
        gotoDirection(direct);

    }
}
computeAccumulated();

//Here is what gotoDirect is supposed to do:

gotoDirect(direct)
{
   /* How can we determine the elevation of the mesh cell to where the
     'direct' points to?
     In the following code, j refers to the elevation of the mesh cell where
     the 'direct' point to.
   */

   store1[j][7] += 1;
   /* increases the number of times that an amount of water
   was passed to it */
}

//Here is what computeAccumulated is supposed to do:

computeAccumulated()
{
    while(j != numberOfmesh)
    {
        store2[j][8] += store1[j][9] * rainAmount;

    }
}


texts: Replacing value of specific pixel using raster mask in ArcGIS Desktop?

texts:
I have two rasters: the first is a basin with all pixels values equal to 0. The second raster (mask) has a single pixel with value = 1. In fact, I need only one raster where all the pixels into the basin must be = 0 and this specific pixel (in the mask) must have the value 1.

How can I do this process using ArcGIS Desktop?

I tried to use raster calculator with the sum "Raster 1 + Raster 2", but it did not work because the result is all the pixels into the basin equal to 1.


texts: ArcGIS Attribute Assistant method not working

texts:
I don't know how many of you are using Attribute Assistant for ArcGIS 10.1, hope many:) I have some difficulties with a couple of methods.

Well, I got the INTERSECTING_EDGE, the GENERATE_ID, and a couple of other simple ones right, but I still can't make it work for other ones.
For example, several methods give me the approximately same error:

  *Field is not specified, empty, or set for all.
                    Trying: CREATE_LINKED_RECORD
                    ERROR: ID or Field to Copy was not found

  INTERSECTING_FEATURE_DISTANCE : ERROR/WARNING: was not found.
  NEAREST FEATURE: ERROR/WARNING: was not found.
  INTERSECTING_FEATURE:  Source Layer not found.
  CASCADE_ATTRIBUTE: ERROR Target layer was not found*

Not sure why "they cannot be found". I hope I wrote the syntax correct.
I can imagine it would be difficult to find the problem for each one of them, but is it somewhere a section on ESRI help resources dedicated to the errors in Attribute Assistant log file, and especially how  can I fix them? Say on a "Field to Copy not found" error what should I do.

So, if anyone is willing to further help me with debugging this I would be grateful.

I tried this on two computers:

1) arcgis 10.0 + attribute assistant addIn for 10.0.

2) arcgis 10.1 +  attribute assistant addIn for 10.1


texts: Random points within a grid cell on a line

texts: I'm trying to create random points on a line within a grid cell. I have a river layer within a state boundary. I have placed a grid over the state boundary. I would like to assign 5 random points along the river line with each grid cell in the state layer.


texts: Which Table is Holding Layers Information and Relationship with Actual Tables in Enterprise Geodatabase

texts:
Can you please let me know which table(s) of Enterprise Geodatabase on top of SQL Server contains contains information about layer ID showing on service URL

"layers": [
  {
   "id": 0,
   "name": "Electric.DBO.Open_Point"
  },
  {
   "id": 1,
   "name": "Electric.DBO.Fuse"
  },
  {
   "id": 2,
   "name": "Electric.DBO.Circuit_Breaker"
  },
  {
   "id": 3,
   "name": "Electric.DBO.Recloser"
  },
  {
   "id": 4,
   "name": "Electric.DBO.Switch"
  },
   ...

and how they are connected (has reationship) with physical tables in Geodatabase?


texts: Python Scripting Displaying Incorrect Names after Join Field

texts:
My problem is that the code is adding numbers to the field names, and this is incorrect.

Problem:

Field should be names VALUE_122, but after running JoinField code, it becomes VALUES_122_123. The values inside are correct.

I have manually gone over the code for the join field, and it does not contain the numbers 123 anywhere.Interestingly, VALUE_123 is not a field in anyone of the 2 tables involved in the joining.

I am using ArcMap 10.1


texts:
Is there a good way to change leaflet basemaps from cartodb.js or leaflet.js? I am able to grab the map object and basemap layer from the layers property in the done() promise.

cartodb.createVis(map, layerUrl, mapOptions)
      .done(function(vis, layers) {
          var map = vis.getNativeMap();
          var basemap = layers[0];
      })
      .error(function(err) {
        console.log("An error occurred: " + err);
      });

Right now I'm thinking I need to be able to grab tileLayer from the basemap object or the map object, then use leaflet's setUrl method as per http://leafletjs.com/reference.html#tilelayer-seturl, but so far no luck there.


texts:
When I install QGIS on windows I get a bunch of custom widgets in  Qt Designer

On OSX I can't figure out how to get those same widgets to appear in Qt Designer. I have set my environment variables according to installation instructions by KyngChaos as so:

export PATH=/Applications/QGIS.app/Contents/MacOS/bin/:$PATH
export PYTHONPATH=$PYTHONPATH:/Library/Frameworks/GDAL.framework/Versions/2.1/Python/2.7/site-packages:/Applications/QGIS.app/Contents/Resources/python

But no luck. Has anyone else gotten this to work?


texts:
I have a table that containing placename regions and another table of roads. I have another table that containing centroid of both placename region and roads.The table structure of third table is placename,Roadname,x,y. I want to update the placename field of the roads which are passing through the corresponding placename region. I used the following query:

Update haryana set placename = placename_region.name from placename_region where haryana.the_geom in (select st_centroid(ch_roads.the_geom) from ch_roads,placename_region where st_contains(placename_region.the_geom,ch_roads.the_geom))

But it went wrong. wrong placenames updated  in placename field. Please help me to solve this


texts: How to import editable layers from other QGIS projects?

texts:
I have a QGIS project that I have spend a lot of time styling over a 80 layers and it seems to have somehow become unstable. I believe that I renamed a folder containing some layer source files which caused QGIS to crash when I reopened the project. After renaming the folder back, the project opens but now will not let me add new layers (says QGIS can't run the Projection Coordinate Selection tool, and nearly every time I use the "Browse" function (in any context), the app crashes. Sometimes I get the Handle Bad layers dialog for two layers, and I can't Browse to find them because it crashes, but then when I reopen it, they are good anyway. Also, I cannot save any Symbol styling as templates without crashing.

Anyway, none of this occurs in a blank project, so I would like to import layers one at a time from my old project into a new project keeping the styling, etc. I see the Embed Layers feature, but it does not allow me to make modifications to the layers afterward. I read about a plugin called ImportLayersfromProject, but I cannot seem to find it on the web anywhere.

I am using QGIS 1.8 on Mac 10.7.5. I would prefer to not upgrade to the latest QGIS until after I complete this big project I am working on. Is there any way of getting the layers imported and being able to edit them too?


texts: Difference in QField Sync's offline editing vs. copy & replace

texts: What is the advantage of configuring a layer in the plugin "QField Sync" for "offline editing" vs just leave it at "copy", do the field work and then copy the edited data from the tablet to the original project folder on my PC? As is understand it, there is no error handling for layers marked for "offline editing" ..so why bother?


texts: ArcSDE PostgreSQL configuration parameters for Spatial Index?

texts:
I wonder why there is not any configuration parameter (DEFAULT Keyword) for spatial indexes in Arcsde for postgres. In ArcSDE for Oracle we have S_STORAGE,S_INDEX_SP_FID,S_INDEX_ALL parameter for configuring spatial index table and its indexes.

When I creating a feature class in postgresql, spatial indexes go into the default tablespace. I don't want this behavior. Am I missing anything?

To elaborate on my question, when I create a feature class in database, the sql that gets executed is (gdo is the user and schema owner):

CREATE TABLE gdo.FC1 (
  objectid  integer NOT NULL,
  name2       varchar(20),
  "name"      varchar(20),
  shape       sde.st_geometry
) WITH (
    OIDS = FALSE
  )
  TABLESPACE gdo_bdata;

CREATE INDEX a19_ix1
  ON gdo.FC1
  USING gist
  (shape);

ArcSDE for postgres, creates spatial tables (FC1), using B_STORAGE (business storage) parameter. But it doens't use any B_INDEX_* parameter for spatial index creation. It seem the user (me!) has no control over spatial index (GIST) creation.

ArcGIS (+ArcSDE) : 10.2.1 and 10.3
Postgresql: 9.2
Geometry Type: ST_Geometry


texts: Error in Converting XY Coordinates to Lat-Long Using r proj4 Library

texts:
I am trying to write a script that will allow my colleagues to easily convert XY coordinates into lat and long. The idea here is to take a spreadsheet with only X and Y available and create a new sheet with latitude and longitude values appended, and to do this without working in a GIS interface at all.

Based on the code provided here, I wrote a script that I thought would do the trick. Most of my colleagues use R, so setting the working directory, reading in the right file, and running the script should be a snap for them. I am testing this with the NYPD's year-to-date complaint data, which already has lat and long but also comes with X and Y coordinates requiring the same projection system (FIPS 3104) as the actual data my colleagues will use:

####libraries####
library(proj4)
library(tidyverse)
library(leaflet)

####working directory####
wd <- ("my working directory")
setwd(wd)

####file location####
file <- ("/NYPD Open Data -- Has State Plane Coords/NYPD_Complaint_Map__no lat lon.csv")

####load data####
xy <- read.csv(paste0(wd,file))%>%
  filter(!is.na(X_COORD_CD),
         !is.na(Y_COORD_CD))

####convert points####
  #assign projection
proj4string <- "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs"
  #project point data
pj <- project(xy, proj4string, inverse=TRUE)
  #create latlon
latlon <- data.frame(xy, lat=pj$y, lon=pj$x)

####check points visually####
leaflet() %>%
  addProviderTiles("CartoDB.Positron", options = providerTileOptions(noWrap = TRUE)) %>%
  setView(lng = -73.9760663, lat = 40.7916269, zoom = 12) %>%
  addCircleMarkers(data = latlon,
                   lng = ~lon, lat = ~lat,
                   radius = 5.0, color="red",stroke=FALSE,fillOpacity = .8

####save output####
write.csv(latlon, file = "test.csv")

Unfortunately, I am not a savvy enough GIS-er to debug the problem I'm having. When I run the script using NYPD's year-to-date complaint data, all of the lat-longs end up being points in Antarctica. Per NYPD's documentation for this data, the X and Y coords should be in the NAD83 New York State Plane-Long Island projection (FIPS 3104), which is the projection string I am using. I have tried swapping X and Y around and setting inverse to false, on the theory that this was a similar error to reversing lat and long when loading data into ArcGIS,but it didn't solve the issue (although it did swap the values from the lat column to the lon column and vice-versa).

I have seen questions about similar errors with pyproj, but don't really understand how to apply the solutions in R (the proj4 project command doesn't have an equivalent argument to pyproj's preserve_units, which seemed to be the problem for pyproj users). I also found this discussion about a bug in the proj4 library that sounds related, but I lack the skill/experience to get any actionable solution from it.


texts: Data View works but Layout View doesn't work

texts:
I created this map on arcGIS 10.2 and was going to finish it in layout view.
The first click on layout toggle gave me a little colored square, second time gave a very tiny square, then nothing.
i tried opening a template then bringing the map into the frame but that did not work either.
Any suggestions?


texts: QGIS 3.4 error: utf-8' codec can't decode byte 0xe8 in position 50: invalid continuation byte

texts:
I have the following error message every time I try to run a GRASS algorithm in QGIS 3.4. I unistalled the software and installed it again but I did not solve anything.

Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 415, in processAlgorithm
Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 372, in executeGrass
for line in iter(proc.stdout.readline, ''):
File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codecs.py", line 322, in decode
(result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe8 in position 50: invalid continuation byte


texts: PostGIS insert failed without raising error

texts:
I'm learning PostGIS through a book "PostGIS in action".

The SQL code written in the book has a lot of issues and I've figured most of them except this one:

Copied from the book chapter:

CREATE TABLE ch03.paris_points(gid SERIAL PRIMARY KEY, osm_id integer, ar_num integer, feature_name varchar(200), feature_type varchar(50));

SELECT AddGeometryColumn('ch03', 'paris_points', 'geom', 32631,'POINT', 2);

INSERT INTO ch03.paris_points(osm_id, ar_num, geom, feature_name,feature_type)

SELECT osm_id, ar_num, geom, tags->'name' As feature_name,

COALESCE(tags->'tourism', tags->'railway','other')::varchar(50) As feature_type

FROM ch03.paris_hetero

WHERE ST_GeometryType(geom) = 'ST_Point'; "

The ch03.paris_hetero contains every column needed. The query returned no error but it just did not insert, nothing has been added into the ch03.paris_points table.

After briefly examining the code, I really can't tell anything that went wrong.


texts: Calculation of surface derivatives at the edge of the raster

texts: I was wondering how does ArcMap10.2 calculate surface derivatives (slope, aspect) at the edge of the raster, since the Horn algorithm is used and uses the 8 neighboring pixels around each pixel. When beyond the raster there are pixels with NoData, do they contribute in the calculation of slope, aspect etc.


texts: Are XY event layers optimized/indexed on-the-fly?

texts:
Oracle 18c; 10.7.1 EGDB; ArcGIS Desktop 10.7.1
I have a non-spatial WORKORDER table that has X & Y columns. The table has 10,000 records.
I’ve created an XY event layer from the table in ArcMap.
The performance of the XY event layer has been surprisingly good, considering the shape field is generated on-the-fly. As far as I can tell, the XY event layer refreshes in the map almost as fast as a feature class does (the difference is negligible). Whereas a spatial view on the same table performs much slower (using an SDE.ST_GEOMETRY or an SDO_GEMETRY column; generated on-the-fly).
Question:
I would like to understand more about how different layer types perform in ArcGIS.
Do XY event layers have any built-in/dynamic optimization mechanisms (like a spatial index) that would account for the surprisingly good performance?


texts: Program to create EMF that conserves transparency

texts: What open source tools allow for the creation of EMFs that conserve transparency in ArcMap?  I tried VeryPDF and Any DWG to Image Converter, neither on actually gave me a vector file that provided transparency.  Adobe Illustrator is a bit pricey when all I want to do is convert CAD blocks to vector symbols to use as ArcMap markers.  Anyone know of a good software for this purpose?


texts: ArcGIS Server 10.6 the field "created user" shows "Esri_Anonymous" when data is sent through collector app

texts:
I have published a service with ArcGIS Server 10.6 as the hosting server. In organisation tab users have been created. Collector app is used to collect site data. Users are signed-in through collector and are able to access the published data and collect site information and send the site data. At database I am able to see the data sent from site but the field "created user" shows "Esri_Anonymous" for all records.

How do I know which record is sent by which user?


texts:
I am using a query where the following line written

CREATE INDEX trgm_idx ON public.roads_r USING gin (name gin_trgm_ops);

Its returning the following error

  operator class "gin_trgm_ops" does not exist for access method "gin"

How can I overcome this?


texts: How to combine vector layer and wms layer based on some constraints?

texts:
I would like to group  wms and vector layer together but with some constraints in openlayers. The requirement is, if lets say vector layer 1 is selected then its corresponding wms layer should also be visualized as the merged layer instead of showing two different overlays. I tried to load the wms layer by checking if vector layer is visible or not using changelayer event with evt.property === "visibility". Is there any way I can do the grouping without using events check?
Thanks


texts:
My organization is looking to manipulate the time slider widget to do the following:

Not show the "play button"-- is there a way to turn this feature off?

Show time in a single year (YYYY) format. We are only able to show time in an interval format of YYYY - YYYY,   despite selecting the YYYY format option.

Any thoughts or advice on these items would be great-- thanks!


texts: Cost raster/friction map and data extraction

texts:
I am in an early stage of least cost paths analysis of a wide area. One part of the operation includes simulation of movement across the grid of previously agreed size of cells and extent of the grid (for reference one cell is 5X5 km).

I have easily enough created the grid (as vector - polygon file), however i would now need to provide friction  values as a attribute input for each cell (each cell is a separate polygon entry), based on the above the sea level values from an underlying DEM model (SRTM in this case). In short this would, ideally, look like this:
- if a 5X5km cell contains more than 50% of surface that is of certain above-the-sea level, i would attribute a value to it from a devised chart of values that describe the friction value for this cell.

Can someone suggest how to proceed with this?

I am quite uncertain how to satisfy the "above 50%" condition.
While i expect to be working with raster data, majority of work will be done while working in table format, with values extracted from various cost raster.

UPDATE:
Following the comments bellow, i worked through the DEM dataset in the following manner: using answers found here:How to reclass a raster with reclassify grid values in QGIS? i reclassified the dem raster into having a number of classes i needed (0 for elevation <=0, 1 for elevation range from 1 to 499 meters and so on - 4 in total). I had then used gdalwarp to obtain a desired resolution of the raster (equal to the resolution of the vector grid - i kept on to the vector, since it seemed to be good way to get the desired data from .shp attribute table into excel, as i need it to be in that form during one of the process steps). In such way 1 grid cell would get the attribute value of the dominant class i.e. elevation range within its boundaries. I went further with calculations of percentages of each elevation class contained within each grid cell, but i need to work on it further - comments on answers bellow point out to the direction im going.
So - in order to get the dataset i need, it is necessary for me to get for an each cell grid (5X5 kilometers size) percentage of each elevation range im having within its boundary (<=0, >= 1 and <=499, >=500 and <=1000m and so on) as a attribute table data, ideally as separate columns. Can this be done?


texts: Extract subset of postgis data to sqlite using SQL query

texts:
I have a large vector mapping dataset held in postgis.  This consists of about 20 individual tables in a single schema.  I want to extract a subset of each of these tables into a standalone sqlite database using ogr2ogr.  This code from the postgis website does this and this works ok for transferring the entirety of each table:

    ogr2ogr --config PG_LIST_ALL_TABLES YES --config PG_SKIP_VIEWS YES -f "SQLite" mydb.sqlite -progress PG:"dbname='mydb' \
  active_schema=myschema schemas=myschema host='localhost' port='5432' user='myuser' password='mypass' " -lco LAUNDER=yes \
  -dsco SPATIALITE=yes -lco SPATIAL_INDEX=yes -gt 65536

However, I would like to only export a subset of each table based on an field called 'Tile_Name'.

I assume I can use the keyword '-clipdstwhere' to extract data based on an attribute, for example where the 'Tile_Name' field is NT, NY, or NZ?

In a standard definition query I would use Tile_Name IN (NT, NY, NZ), but this obviously isn't the correct syntax as I get the error "ERROR 1: No field definitions found for 'IN', is it a table?".

There are also options for clipping the data by a boundary or bounding box - might this be an easier way?


texts:
I am embedding a map I made in Tilemill/Mapbox onto a Wix website. Items of the map such as legend, title, and description are not automatically included in the embed. I learned how to add the legend by simply adding ,legend, in the code below. However, I do not know how to add the title and description of the map as you see it here: https://a.tiles.mapbox.com/v3/yolivens.hgck7b4h/page.html?secure=1#4/40.85/-94.48

Here is my current  code:

What would I need to do to add the title and description? I am above all concerned about the title.


texts: How do I change the extent of a raster file in ArcGIS in order to remove nodata values?

texts:
I am currently working with a netcdf file which I export to a GeoTiff and then layer over a 3D globe, the problem is a get a white line showing up at the dateline, straight through the pacific. When I look at the data I see that the final line of pixels on the east (right) side of the netcdf and GeoTiff file is set to nodata, hence the white line. My question is how do I remove this line of nodata pixels and rescale (I'm not sure this is the right terminology) the raster so that its extent is exactly -180, 180, -90, 90?

Here are the details of what I am doing. I am opening the netcdf file with ArcGIS using the "Make netcdf raster layer". I extract the relevant data and export to a GeoTiff. I get the following extent: -179.999946594238, 179.846244812012, -90.0677337646484, 90.0219573974609 with 4320, 2160 columns and rows and a cell size of (x,y) 0.083297729, 0.083343506. I however want the extent to be exactly -180, 180, -90, 90 with 4320, 2160 columns and rows and a cell size of (x,y) 0.0833333, 0.0833333. I have tried changing the extent in the layer properties, as well as, changing the raster size or cell size when exporting the GeoTiff. Nothing seems to work.


texts:
I am software developer and have been working on database
and Java as major platform for development. I am using Mapinfo for GIS analysis.
Now, the system needs to show trouble tickets in Mapinfo that are achieved from
Maximo software. How is it possible to connect these two application?


texts: Georeferencing photographs using ArcPy with ArcMap

texts:
I have 563 photographs that have been taken as part of a site survey. The only location information that I have is a single GPS location (lat,long) at which the photograph was taken. I am looking for a Python Script (ArcPy) that will use the GPS location of each Photograph (Point Feature Class) and then based on a given width and height georeference each photograph. It's important to read the raster properties to determine if the image is portrait or landscape.

Can you provide sample code in ArcPy?


texts: CARTO/Leaflet: How to improve map performance?

texts:
I'm having problems getting my map to load data consistently quickly and smoothly. The buttons on my map execute sql queries to display data per category and by year. When I click between them, especially the radio buttons and year options, the data is sometimes displayed quickly, but most of the time (for me at least) there will be a pause for 2-3 seconds before the map refreshes, or an 8+ second pause at its worst.

Because I want to eventually publish my map for public use, I really want to fix this issue first - how can I speed up or optimise my map's performance to prevent these long waits between reloads?

Link to the map: https://jsfiddle.net/Sidders/w9rzhLbr/


texts:
I have set pgrouting and inserted my data with osm2pgrouting .

I would like to know how can I find out which way is one way.
Also I read that if I assign the reverse_cost of an one way to -1 only the this road is considered as one way, is it true?


texts: Embed / Iframe on Google Streetview, disable click but keep the drag

texts:
This is a question due an older post from Cesare.

In my case I want to disable click but keep the drag (so you can explore the 360-degree camera) in an embedded (with Iframe) map on Streetview.

Here is the map I am trying to embed and is from the International Space Station.

I tried with Cesare solution (JSFiddle here) but (maybe because is the International space station and there are some legal rights) it shows a black screen with the note: Image may be subject to copyright. Therefore now I am trying with iframe. Here is my code on JSFiddle.

<!DOCTYPE html>
<html>
  <head>
    <meta charset="utf-8">
    <title>Street View controls</title>
   <style>
      /* Always set the map height explicitly to define the size of the div
       * element that contains the map. */
      #map {
        height: 100%;
      }
      /* Optional: Makes the sample page fill the window. */
      html, body {
        height: 100%;
        margin: 0;
        padding: 0;
      }
    </style>
  </head>
  <body>
      <div class="map-container">
        <iframe src="https://www.google.com/maps/embed?pb=!4v1518736193836!6m8!1m7!1szChzPIAn4RIAAAQvxgbyEg!2m2!1d29.56028532748453!2d-95.08539144618427!3f213.12!4f0.18000000000000682!5f0.7820865974627469" style="margin-top:-100px" width="100%" height="1000" frameborder="0"allowfullscreen >
   </iframe>
        <style>
        .map-container {
        </style>
  </body>
</html>


texts: Create grid where points is centroid? QGIS

texts:
I want a grid around my points(100 x 100 m), where points are center. Creating a grid from

Vector --> Research Tools --> Vector Grid

creates grid, but not centered around my points.

How do I do this?


texts: Earth Engine: Can I export a feature collection as a tif?

texts:
I am using Random Forest in Earth Engine to classify land use from Sentinel 2 images. I want to eventually have a raster file of the classification with each pixel classified using the RF model. At the moment my classification is a feature collection, can I convert/export this as a tif file?

Here is the code and the link to the code:

ee_code

In the code below "classified" is the clssification feature collection.

Map.addLayer(guan_ca);

var vis = sen2.filter(ee.Filter.calendarRange(2017,2019,'year'))
    .filter(ee.Filter.calendarRange(01,01,'month'))
    .filterBounds(guan_ca)
    .map(function(image){
    return image.clip(guan_ca);
});

var vis = vis.mean();

var trueColour = {
         bands: ["B4", "B3", "B2"],
         min: 0,
         max: 3000
     };

     // Add image to map
Map.addLayer(vis, trueColour, "true-colour image");
// Define false-colour visualization parameters
var falseColour = {
    bands: ["B8", "B4", "B3"],
    min: 0,
    max: 3000
};

// // Define false-colour visualization parameters
var falseColour = {
    bands: ["B8", "B4", "B3"],
    min: 0,
    max: 3000
};

// Add image to map
Map.addLayer(vis, falseColour, "false-colour image");

// // Create the training data

var newfc = ca.merge(water).merge(ca).merge(baresoil).merge(crops).merge(urban).merge(clouds).merge(shadows).merge(forest).merge(lightcover);

print(newfc);

var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8'];
var training = vis.select(bands).sampleRegions({
    collection: newfc,
    properties: ['landcover'],
    scale: 30
});

print(training);

// Train the classifier, in this case, I use the Random Forest algorithm
var classifier = ee.Classifier.randomForest(100)
   .train({
           features: training,
           classProperty: 'landcover',
           inputProperties: bands
});

// Run the classification
var classified = vis.select(bands).classify(classifier);
print(classified);

// Display the classification
Map.centerObject(newfc, 11);
Map.addLayer(vis,
{bands: ['B2', 'B3', 'B4'], max: 0.3},
'Sentinel 2 image');
Map.addLayer(classified, {min: 0, max: 8, palette: ['126CC6', 'EBEB07', '613D05', '2CBE17', 'B9BFB8', 'FFFFFF', '000000', '076105', 'D47606']},
'classification');
Map.addLayer(newfc);

Export.table.toDrive({
  collection: newfc,
  description: ('landuse_trainingData'),
  fileFormat: 'KML',
});

Export.table.toDrive({
  collection: classified,
  description: ('land_use_cassification'),
  fileFormat: 'KML',
});

// Error estimation
var trainAccuracy = classifier.confusionMatrix();
print('Resubstitution error matrix: ', trainAccuracy);
print('Training overall accuracy: ', trainAccuracy.accuracy());

var withRandom = training.randomColumn('random');

// We want to reserve some of the data for testing, to avoid overfitting the model.
var split = 0.7;  // Roughly 70% training, 30% testing.
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));

// Trained with 70% of our data.
var trainedClassifier = ee.Classifier.randomForest().train({
  features: trainingPartition,
  classProperty: 'landcover',
  inputProperties: bands
});

// Classify the test FeatureCollection.
var test = testingPartition.classify(trainedClassifier);

// Print the confusion matrix.
var confusionMatrix = test.errorMatrix('landcover', 'classification');
print('RF testing error matrix', confusionMatrix);
print('RF testing accuraccy', confusionMatrix.accuracy());
print('RF testing accuracy', confusionMatrix.array());

var matrix = ee.Feature(null, {matrix: confusionMatrix.array ()});


texts: How to quickly show areas on a map (with OpenStreetMap?)

texts: Yesterday I wanted to quickly indicate some areas (lat/lon rectangles) on a map, like OpenStreetMap. And I couldn't figure out how to do it. Any suggestions? It doesn't have to be interactive, just a bitmap image of the map, with the areas showing as rectangles, would be fine.


texts: Can anyone tell me why my field calculator is returning this error?

texts:
All I am trying to do is pull into a text field, a section of another text field via:

fieldB = !fieldA![some_index:other_index]

However I keep receiving the error below and I can't figure out why. My python syntax is correct and I have tried it in IDLE on an example string with the exact contents of a row in my feature class and it works fine.

I am using ArcGIS 10.0. The error is:

Executing: CalculateField AA_TCPs_Complete_060613 Cross_rd_1 !PopupInfo![5:10] PYTHON_9.3 #
Start Time: Fri Jun 07 13:31:03 2013
ERROR 000539: <type 'exceptions.SyntaxError'>:invalid syntax(<expression>, line 1)
Failed to execute (CalculateField).
Failed at Fri Jun 07 13:31:03 2013 (Elapsed Time: 0.00 seconds)


texts: Combining Hover and Click events in OpenLayers

texts:
I can specify event callbacks in an OpenLayers layer:

var vectorA=new OpenLayers.Layer.Vector("A");
vectorA.events.on({
  "featureselected": callback
});

and I can specify event callbacks in a control's constructor:

highlightControl = new OpenLayers.Control.SelectFeature(vectorA, {
  renderIntent: "temporary",
  hover: true,
  highlightOnly: true,
  eventListeners: {
    featurehighlighted: hoverFeatureCallback,
    featureunhighlighted: unhoverFeatureCallback
  }
});

Can I mix them properly?

Basically, one layer needs both select and hover events, so I'd have two controls, one that hovers, and give each their own callbacks (like the second example above.)

AND I have multiple layers that need select events, so in that case I'd define the events in the layers themselves, like the first example. However, I can't use the first example with hover and select because, in OL's infinite wisdom, hover events are really just select events.

How can I have two layers respond to hover events and one of them also respond to select events?


texts:
Can anyone direct me to GIS data for malaria mortality rate?

I am looking for countries in South Asia, but not looking for country level numbers but smaller level numbers like district level.


texts:
I'm using ArcGIS 10.3, and I have 20 shapefiles, each representing 1 year worth of data over the same study area, that contain postal code point data, with each point having some NDVI values (growing season mean, growing season peak etc.).

What I would like to do is a time series analysis to see how vegetation changed over the course of 20 years, and if possible to identify areas where the change was significantly different than what would normally be expected. The reason for this is that in one year there was a large storm which wiped out a good chunk of the vegetation in my study area. So I want to identify those areas that were hardest hit by the storm.

I've looked online and through the GIS stackexchange and while I've found some examples of what I want to do (I think) like here:

http://r-video-tutorial.blogspot.com/2016/07/time-series-analysis-in-arcgis.html

or maybe like this:

https://www.umesc.usgs.gov/management/dss/curve_fit.html

They are always using raster data, which I don't have access to.


texts: Creating KMZ files from QGis using GarminCustomMap plugin

texts:
I am new to GIS and QGIS, and really lost with some basic QGIS concepts.

I have installed the GarminCustomMaps plugin, and by usign the default plug-in settings I am able to generate a single tile of the area that is visible on the screen in QGIS.

What I'd like to do instead is generate a KMZ file that contains multiple tiles covering an area I select in QGIS.

My source in QGIS is a WMTS layer from the French geoportail site:
   https://wxs.ign.fr/<access key>/geoportail/wmts?SERVICE%3DWMTS%26REQUEST%3DGetCapabilities
   The source is multi-resolution - as I zoom in I get different sets of images.

Here are some questions I couldn't figure out:
- what is the zoom factor of the plugin used for exactly? If I use anything by the default value of 1.0, I end up with a blank KMZ file (checking in Basecamp)
- The plug-in seems to export whatever is visible on the screen. If I make a larger area visible, the layer switches to lower scale maps (e.g. 1/50,000 instead of 1/25,000): is there any way to display a larger area without the scale changing?

Is there any way to define an arbitrary extent to be exported by the plug-in?

I am probably not using the right terminology, but I have gone through all the documentation I could find and I am still pretty confused.


texts: How do you set the extent of a CartoDB 'embed' map?

texts:
Using the basic CartoDB 'embed' option, how is the extent stored / saved / set for a map?

For example:

A CartoDB map (formerly 'visualization') has three potential pages on my blog / website (WordPress) for each State I want to show from my map.

I can't seem to figure out how to use the same 'map' for 3 different extents (views) of the same map.

There is no option to set the extent, and it seems to be set to the last place you panned/zoomed to on the map itself - with no option to save.

Am I missing something or is there a better workflow for this?

And since I'm using WordPress, I can't use CartoDB.js, and the same problem seems to apply to using the 'simple link' for the WordPress plugin.


texts: Copying the shape of a polygon but paste it on a different location in QGIS

texts: I am labeling solar panels and often they have exactly the same shape. Is it possible to copy the shape of a polygon and paste it on another location in my shapefile?


texts: Map not showing while following tutorial for QGIS

texts:
1I downloaded and installed QGIS 3.4 and tried to load an example shapefile layer as described in this tutorial using the provided shapefiles.

But the images are not showing up. I tried to uninstall/reinstall but it didn't help. What am I doing wrong?
Here is what I'm getting.
[]
Here is the attribute table


texts: Importing OSM/Google Maps to QGIS but removing/excluding extraneous areas to limit geographically use of OSM?

texts:
I am looking to place an OpenStreetMap/Google Maps layer underneath the rest of my shapefiles for referencing purposes.

Let's say for instance, that I am looking at the Swiss Cantons and cities, and I want the Street Maps underneath these layers. What I have been doing is importing OpenStreetMap from the OpenLayers plugins. However, this imports a world map and rendering is slow.

Is there a way to import the OSM/Google Maps layer and then exclude extraneous or superfluous areas?

If I am only interested in Switzerland I don't necessary need Uganda, Botswana, Kazakhstan and Indonesia as well.


texts: Extract the county (or village) level data of forest cover change from Google earth engine

texts:
I am currently working with the Hansen data (Global Forest Change) in Earth Engine. And I would like to extract the county (or village) level data of forest cover change in China from Google earth engine.

The problems is that I'm not sure if I can extract the county (or village) level data by using google earth engine, also I don't know what I need for specifying the boundaries if I can extract the county level data. I've tried to read most of the GUIDES and use the Featurecollection to specify but couldn't make it.


texts: Use GeoExt to edit Features attributes on a Map

texts:
I'm new to geoext and the documentation doesn't really help me.

I've got a map with points on it. These points have attributes like 'name'.
I coded a map with geoExt that has a vector layer which is poplulated by features(points) from a server. I can add features to the vector layer (not to the server) and select them.

Now I'd like to select the features and edit their attributes using a window showing up. Saving them and in the next step saving them to the server.

I've been trying to do this for a few hours now, but I can't get used to ext js/geoExt "philosophy" and I can't find examples that help me do what I wan't.

Could somebody give me advice how to proceed?

Edit:
i worked on the code and i'm now able to read the attributes of a feature an show it in a popup.
i still like to know how to edit the values and save it to the server

Here's my (new) existing code:

var app, items = []

var osmLayer = new OpenLayers.Layer.OSM();
var vectorLayer = new OpenLayers.Layer.Vector("Vector Layer");
var map = new OpenLayers.Map({
    projection: new OpenLayers.Projection("EPSG:900913"),
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    // todo: benutzerdefinierbar
    units: 'm',
    controls: [
        new OpenLayers.Control.ZoomPanel()
    ]
});

Ext.onReady(function() {
    app = new Ext.Viewport({
        layout: "border",
        items: items
    });
});

// registering a event that calls popup creation
vectorLayer.events.on({
    featureselected: function(e) {
        createPopup(e.feature);
    }
});

// create a popup showing the properties of a feature
var popup;
function createPopup(feature) {
    // so kann man properties schreiben
    // feature.attributes.placename = "testtestPlacename";
    popup = new GeoExt.Popup({
        title: "Features Properties",
        location: feature,
        width:200,
        // steht unter dem grid
        //html: "hier steht html text",
        maximizable: false,
        collapsible: false,
        unpinnable: false,
        items: //new Ext.grid.PropertyGrid({
            new Ext.grid.EditorGridPanel({
            xtype: "editorgrid",
            //title: "testTitle",
            source: feature.attributes,
            autoHeight: true
        })
    });
    popup.show();
}

// populate vectorLayer with features from remote source
var store = new GeoExt.data.FeatureStore({
    layer: vectorLayer,
    proxy: new GeoExt.data.ProtocolProxy({
        protocol: new OpenLayers.Protocol.HTTP({
            url: '...',
            format: new OpenLayers.Format.GeoJSON({})
        })
    }),
    fields: [
        {name: 'title', type: 'string'},
        {name: 'description', type: 'string'}
    ],
    autoLoad: true
});

items.push(new GeoExt.MapPanel({
    region: "center",
    map: map,
    title: 'xyz',
    layers: [
        osmLayer,
        vectorLayer
    ],
    tbar: buildMapToolbar(map, vectorLayer)
}
));


texts:
When I booted up my windows machine this morning my menus had switch from English to what I think may be Russian and I can't read the menu options. I have tried unistalling, rebooting and then reinstalling but this has not helped. Can someone please advise me how to change them back.

Thank you


texts: how to save drawable entities with leaflet.draw/draw polygons into my Postgis database with getting also the coordinates

texts:
I am very new to creating webmaps with Leaflet and I am trying to save my shapes I draw into my postgis db. I would prefer to be saved as .shp or any default format (i.e. json) with its geometry field. I seek a step-by-step code like a small demo.

the script till now:

var drawnItems = new L.FeatureGroup();
    map.addLayer(drawnItems);

    map.on('draw:created', function (e) {
        var type = e.layerType,
        layer = e.layer;
        drawnItems.addLayer(layer);
    });

    polygon_options = {
        allowIntersection: false,
                drawError: {
                    color: '#b00b00',
                    timeout: 1000
        },
                showArea: true,
                metric: false,
                repeatMode: false,
        shapeOptions: {
            stroke: true,
            color: '#6e83f0',
            weight: 4,
            opacity: 0.5,
            fill: true,
            fillColor: null,
            fillOpacity: 0.2,
            clickable: true
        }
    }

    function drawPolygon(){
            var polygonDrawer = new L.Draw.Polygon(map, polygon_options);
        polygonDrawer.enable();
    }


texts: change the direction of the marker of a line in QGIS

texts:
I have a map with three lines.  Each has a marker on it.

I would like to reverse the direction of the topmost marker.  So I open Style dialog and enter the expression "direction" * 180.  The value of direction of the topmost line is 1.  Others are 0.

The result is shown below.

Two issues here:

The angles of all markers are the same.
The direction of the topmost marker is not reversed.  (Is there a way to access a default value of 'Angle' property?)


texts: Ignoring NA values when computing variogram in R?

texts:
I am trying to compute a variogram with the R gstat package:

 > variogram(yield_cleaned~1, data= newdf)

but it turns out:

 Error in .local(obj, ...) : NA values in coordinates

My dataset has NA values. How can I tell the program to ignore them?


texts: average weighted value by region with use of Voronoi in Qgis

texts:
I have a vector layer of Voronoi polygons. This layer has three columns -longitude, latitude, Z-value.

I need to cut polygons by mask of administrative area and then I need to calculate average-weighted Z-value (AWZ) for every administrative area with use of formula:

AWZ = SUM(Z-value*Square_of_polygon_inside_area)/SUM(Square_of_polygon_inside_area)

How should I solve this problem?


texts: Is it possible to create a QGIS Windows Installer from my own source on Ubuntu?

texts:
I would like to create a standalone Windows installer for QGIS from the modifications I have made locally to my source in Ubuntu.  I have done various things like change the splash image and added some functionality here and there.  I can build and successfully run this version on Ubuntu.

I am able to use the creatensis.pl script to create a Windows Installer executable...

cd <source>/ms-windows/osgeo4w
./creatensis.pl qgis

but it seems to use the qgis master version as its source.  When I install the resulting executable on Windows, it shows none of the changes that I made.


texts: Testing spatial autocorrelation for presence-absence records?

texts:
I want to test the spatial autocorrelation for presence-absence data records, I know that the Moran's I test cannot be used here since I have categorical data.
I was thinking about the join-count statistic, but is there any work that applied this statistic to the presence-absence data?
I have data points and I don't know how to use this statistic here.


texts: How do I Convert from LayerDescription.SelectionFeatures to GraphicsLayer?

texts:
How do I convert (in a clear and concise way) from

ESRI.ArcGIS.ADF.ArcGISServer.LayerDescription.SelectionFeatures

to

ESRI.ArcGIS.ADF.Web.Display.Graphics.GraphicsLayer

I want the selected features and the attributes.

side note:
I got the the selection features via FIDSet returned from a SOAP based MapServerProxy.QueryFeatureIDs query using a buffer geometry to pull intersecting features.


texts: Coordinate transformation results are different between Proj4JS and pyproj

texts:
I try to make a coordinate converter using Proj4JS for mobile app, I use this source code. Before I tested conversion for WGS84 geographic and WGS84 UTM, there was not problem with it. I tested conversion from TUREF TM zone:27 to ED50 TM zone:27 that input coordinates are Easting:507768.06 Northing:4354677.64. For proj4js output is; Easting:507819.315 Northing:4354847.596 But when I tested it for epsg.io or pyproj output was Easting:507815.93 Northing:4354857.26.

Edit: I use proj4 text parameters for pyproj and proj4js.

TUREF / TM27: +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
ED50 / TM27: +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=500000 +y_0=0 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs

What am I missing?


texts:
Given a point and few lines, how should I go about finding the nearest line to the point?

I am aware of Finding nearest line to point using ArcGIS Desktop (ArcObjects/ArcPy)?, but that is in ArcGIS and uses its functions.

The "duplicate" question doesn't answer this question properly, it's just a pointer to MMQGIS.


texts:
I want to display the labels for each feature, but I get the "Uncaught ReferenceError: Style is not defined" error. How do I fix it correctly?

     var map = new ol.Map({
        target: 'map',
        layers: [
            new ol.layer.Tile({
                source: new ol.source.OSM()
            })
        ],
        view: new ol.View({
            center: ol.proj.fromLonLat([102.291163, 12.809588]),
            zoom: 9
        })
    });
    var labelStyle = new Style({
        text: new Text({
            font: '12px Calibri,sans-serif',
            overflow: true,
            fill: new Fill({
            color: '#000'
            }),
            stroke: new Stroke({
            color: '#fff',
            width: 3
            })
        })
    });

    var tambonStyle = new Style({
        fill: new Fill({
            color: 'rgba(255, 255, 255, 0.6)'
        }),
        stroke: new Stroke({
            color: '#319FD3',
            width: 1
        })
    });

    var style = [tambonStyle, labelStyle];

    var layer_source = new ol.source.Vector({
        url: 'geojson/tambon_wgs84.geojson',
        format: new ol.format.GeoJSON()
    });

    var layer_tambon = new ol.layer.Vector({
        source: layer_source,
        style: function(features) {
            labelStyle.getText().setText(features.get('TAM_NAM_T'));
            return style;
        },
        declutter: true
    });

    map.addLayer(layer_tambon);


texts: Extract one Feature from Feature Collection

texts:
I have a feature collection consisting of three individual features, each of which represents an AOI in an area. Now I want to get each feature separately to load the appropriate Sentinel2 data for it. First of all I load the FIRMS (MODIS) Data and convert day to date. Then I clip the mosaic to my AOI.
This is my code:

  var FIRMS_filtered = ee.ImageCollection('FIRMS')
 .filterDate('2018-01-01', '2018-12-31')
 .filterBounds(AOI)
 .select('T21')
 .map(function(image) {
 var date = ee.Date(image.get('system:time_start')).format("YYYY-DDD");
 date = ee.Date(date);
 return image.set('date', date);
 })
.map(function(image) {
var ID = ee.String(image.get('system:index'));
 return image.set('ID', ID);
})
.map(function(image) {
 var date = ee.Date(image.get('system:time_start'));
 return image.set('date_time', date);
 });

Then I create a Vector from the result and a Buffer.

  //Creation AOI mask for vectorization
  var T21 = FIRMS_filtered.mosaic().clip(AOI).select('T21');
  var zones = T21.gt(0);
  zones = zones.updateMask(zones.neq(0));

  //Creation of vector
  var vectors = zones.addBands(T21).reduceToVectors({
  geometry: AOI,
  crs: T21.projection(),
  scale: 500,
  geometryType: 'polygon',
 eightConnected: false,
 labelProperty: 'zone',
 reducer: ee.Reducer.mean()});

//Function for creating buffer
var buffer = function(feature) {
return feature.buffer(2000);};
//Creation buffer
var bounds_AOI = vectors.map(buffer);

The result is a Feature Collection with 3 Features. Now I want to get each of this Features single.
How can I do this?

Here is the code

https://code.earthengine.google.com/def19f2aed544401711895df71be76d7


texts: Inputting JSON file as variable in ModelBuilder using ArcGIS?

texts:
I am using ArcGIS 10.3 and trying to create a Geoprocessing service using ModelBuilder, the input of my GP service is in JSON format.

Is it possible to set up the ModelBuilder to accept the JSON input as a variable?


texts: Get latitude and longitude differences by clicking in two locations on the map

texts:
Is it possible in QGIS 3.x to click in two places on the map and to get back the X and Y distances between these two locations (latitude and longitude differences).

I was hoping there is a plugin that does that but couldn't find it.

This would be very helpful.


texts: gdal.Warp ignores some pixels when reprojecting it

texts:
I have a tif in EPSG3035, pixel size 100x100m (very important for the study).
I use gdal.Warp to reproject the tif to EPSG4326 for getting lat/long of pixels along with their values. i am forcing the function to preserve the array area since some calculations follow.
The tif statistics are the following:
in_raster sum = 149580
out_raster sum = 136220
When visualizing the rasters (reds are out_raster) I notice there are some pixels not taken in the out_raster
The code I am using is the following:
def cs2wgs84(input_raster, output_raster):
raster = gdal.Open(input_raster)
band = raster.GetRasterBand(1)
band.SetNoDataValue(0)
band.FlushCache()

rwidth = raster.RasterXSize
rheight = raster.RasterYSize

woptions = gdal.WarpOptions(width = rwidth ,height = rheight,
    dstSRS='EPSG:4326',overviewLevel = 'NONE', dstNodata="0", resampleAlg = "near")
gdal.Warp(output_raster, input_raster, options=woptions)
band = None
raster= None

The results of the raster statistics are the following:
in_raster:
in_raster statistics:
Results:
{'MAX': 2267.493408203125,
'MEAN': 10.610063969353307,
'MIN': 0.0,
'OUTPUT_HTML_FILE': 'C:/Users/diana/AppData/Local/Temp/processing_tliQAH/a43f435d8d294ad3be80a681046b20bc/OUTPUT_HTML_FILE.html',
'RANGE': 2267.493408203125,
'STD_DEV': 81.50870519877743,
'SUM': 149580.68183994293,
'SUM_OF_SQUARES': 93655802.2197855}

out_raster statistics:
Results:
{'MAX': 2267.493408203125,
'MEAN': 9.321215707322825,
'MIN': 0.0,
'OUTPUT_HTML_FILE': 'C:/Users/diana/AppData/Local/Temp/processing_tliQAH/b68b22ce98064a5b9d20e616ef06b6fa/OUTPUT_HTML_FILE.html',
'RANGE': 2267.493408203125,
'STD_DEV': 76.70612274249875,
'SUM': 132062.9841413498,
'SUM_OF_SQUARES': 83356209.21407527}


texts:
I am curious if there is a way to link a PostGIS table with a Google Sheets?

Both tables (Post and Google) have a unique ID that can be used as a linkage.  I am hoping to be able to build a view that will display information from both tables into a single view.


texts: How to extract phenology metrics from every singe pixel in an image with TIMESAT

texts:
I want to extract phenology metrics of every pixel in an image. I have 92 MODIS images (4 years and 23 images per year) of a zone, and I want to extract the metrics data (as a text or ASCII file) to use it in a model to map a crop in that zone.

I used .tif images but had some problems to load the image correctly so I transformed them to ENVI (with QGIS), the image load correctly in TSM_imageview but I can not process the list of the time series in TSM_imageview neither in TSM_GUI, I tried using other format (IDRISI, .img) but it does not work.

When I load the image list in TSM_imageview I only receive a message in the command window:

   Error using dataread
   Trouble reading integer from file (row 1, field 1) ==> ÿ_9

   Error in strread (line 48)

   Error in openfilelist

   Error using waitforallfiguresclosed (line 20)
   Error while evaluating UIControl Callback

Then when I tried to load the list in TSM_GUI, I write the informations (number of rows and columns, years) and the number of images per year (that supposed to load automatically) disappears, when I write it manually the number of years disappears.

I do not know what happens, anyone can give me an advice?


texts: Field calculation based on the field name by using arcpy

texts:
The shapefile table is like this.

I want to do the sum of cells when their field names meet the condition. In this case, I want to sum the cells to the TotalCT field if the field name contains the string "CT". The manual way is to do it like this.

But this is only a very small example. I wonder how can I realize it with Python or Arcpy code?


texts: Precipitation or snow cover raster data for Scandinavia

texts:
I am looking for free precipitation or snow cover data in Scandinavia mountain during the winter 2011-2012. The WorldClim dataset is not enough acute. I would like 500 m scale at least. There is NASA's MODIS, but MOD10 (MOD10A2 for instance) are not available. Does anybody know where I can find these kind of data?
Thanks.


texts: Is it due to irregularity of geoid, we cannot convert the height between spherical datum and local datum by cs2cs?

texts:
Is it due to irregularity of geoid, we cannot convert the height between spherical datum and local datum (which uses locally defined geoid and height) by cs2cs?

What I mean is that, by using cs2cs, it only accepts x and y / lat and lon as input. The z / height data is not accepted for the conversion program. Does it mean that z / height cannot be directed converted?

What if I convert between spherical data, can it be converted directly by using projection formula?


texts: Calculating quantile values over a stack of images with grass gis. Help needed on r.series

texts:
I need to calculate the values for the 25%, 50%, 75% and 90% quantile over a series of rasters (time series of ndvi values from landsat).
I am exploring the r.series command in Grass GIS, but the documentation is very poor and it doesn't seem to fit my needs.

I would need to calculate out of the overall distribution of values in the 12 monthly ndvi images the indicated quantile. I just need the values, not a map.

Any idea on the matter? Anyone can actually explain me what method=quantile actually does in r.series on grass gis?


texts: Converting geometry stored in SQL Server to Latitude/Longitude using ArcObjects?

texts:
I do not want to use ArcMap or any other ESRI tool except ArcObjects to accomplish this task. Another open source solution would also be fine.

I need to convert geometry data from ArcGIS using SQL server to WGS84. I have some code to do a specific projection but I don't know how to apply this to an unknown projection without manually changing the code. As I said the data I need to convert is in an ESRI SQL server. I've looked through the SDE tables and see basically how ESRI stores the data. My question is what table holds the EPSG for the stored projection for a given feature class table so that I can properly provide the esriSRGeoCSType? What esriSRProjCSType do I use for standard lat lon (WGS84) North America?

My goal is to write a simple routine that you can point to a ESRI SQL server, the user selects the feature class table they want to convert and the routine spits out the data along with the geometry projected to WGS84. The original projection is unknown and will have to be retrieved from the ESRI SQL server. Once the user selects the table, the selection is stored the process will continue to run as a service and the process will happen on a timed interval so it needs to be done without further user interaction. Maybe I'm going about it all wrong, I'm a little lost at this point.

UPDATE - RESOLVED
Here is the working code to reproject one coordinate system to another. In this case I'm using NJ state plane to WGS84 Lat and Long.

            int inEPSG = 3424;
            int outEPSG = 4326;

            ISpatialReferenceFactory2 srFactory = new SpatialReferenceEnvironmentClass();
            ISpatialReference ProjectedSR = srFactory.CreateSpatialReference(outEPSG);
            ISpatialReference GeographicSR = srFactory.CreateSpatialReference(inEPSG);

            double xCoord = 547326.03872482851;
            double yCoord = 541395.50404749066;

            IPoint geographicPoint = new PointClass();
            geographicPoint.PutCoords(xCoord, yCoord);
            geographicPoint.SpatialReference = GeographicSR;

            geographicPoint.Project(ProjectedSR);

            var _gcsLat = geographicPoint.Y;
            var _gcsLong = geographicPoint.X;

To get the current EPSG for a particular geometry column in an ESRI geo-enabled database you can retrieve it from the table SDE_geometry_columns. It has a FK relationship to SDE_spatial_references using the srid column. The column auth_srid in SDE_spatial_references will have the EPSG used for that column's coordinate system.

As I thought it was very simple.


texts: How to find and remove invalid geometry in a multilinestring python geodataframe?

texts:
I have original data as a 700000's rows excel, with geometry data in hexwkb format.
I want to plot it and make some spatial operations on it like intersect, distance and buffer, on a jupyter notebook.
Half of the geometries are null. Part of the remaining ones have apparently invalid geometries.
I had succeed in plotting the first 100 lines, apparently with fine geometries, but I'm struggling to plot all the data.
This code worked for lines with no geometry problems:
Data example
id (int64) | geom (object)
1 | '0105000020E610000001000...
import pandas as pd
import geopandas as gpd
from shapely import wkb
from shapely.geometry import MultiLineString

df=pd.read_excel(r'trips.xlsx')
trips= df.dropna(subset=['geom'])

# create column with converted geometry

def convert_to_hexwkb(geomString):
    geometry = wkb.loads(geomString, hex=True)
    return geometry
trips['geometry']=trips['geom'].map(convert_to_hexwkb)
trips.drop('geom', axis=1, inplace=True)

# convert to GeoDataFrame

gtrips=gpd.GeoDataFrame(trips)

print(gtrips['geometry'][0])
>> MULTILINESTRING ((-9.1474504470825 38.744022369385, -9.147529602050801 38.744121551514, -9.14...

# show trips map

gtrips.plot()
>> got a nice map here

But when I run this to the entire data I get the following error given from the shapely -> geos library, during conversion of string to hexwkb:

WKBReadingError: Could not create geometry because of errors while reading input.
IllegalArgumentException: point array must contain 0 or >1 elements

I couldn't find how to apply Geopandas.GeoSeries.is_empty function to find empty geometries, because it was not a geodataframe at that point.
Same for is_valid function from shapely library, since the field was not geometry type.
I've tried to delete exceptional rows with this:
for i, item in enumerate(trips['geom']):
   try:
    trips.at[i,'geometry'] = wkb.loads(item, hex=True)
      #wkb.loads(item,hex=True)
   except:
    trips.drop([i], inplace=True)

And got the errors:

IllegalArgumentException: point array must contain 0 or >1 elements
WKBReadingError: Could not create geometry because of errors while reading input.
another exception occurred: KeyError: '[791] not found in axis' - This related with how the for cycle  was constructed probably.

If I run this for debug:
def convert_to_hexwkb(geomString):
    try:
        geometry = wkb.loads(geomString, hex=True)
        return geometry
    except:
        return np.nan

trips['geometry']=trips['geom'].map(convert_to_hexwkb)

I get these errors, multiplied for several lines:

ParseException: Unexpected EOF parsing WKB
ParseException: Premature end of HEX string

Any clue of what should I do to clean geometries of my dataset?


texts: How does the Phantom 3 4k UAV assign the coordinates/geotag to each photo?

texts: I read that the Phantom 3 used to use a barometric reading but now it uses GPS. Does it take a new reading from the GPS/Glonass satellites every time it takes a photo or does it base every photo from the initial coordinate it established when it is turned on before flght? Furthermore, does it take the pitch, role, yaw into account when calculating the coordinate for the photos?


texts: Projection for polygon intersection confined into one hemisphere?

texts:
I have polygons where each vertex is a (lon,lat) geographic coordinate. I want to find out if they intersect or contain each other. The polygons are guaranteed to never exceed a hemisphere (so they'll never 'wrap around' beyond a plane that cuts the Earth in half), but they can vary dramatically in scale.

I think one way this might be done is by transforming the polygons using some specific projection, doing the math for intersection in cartesian coordinates and then 'unprojecting' back to geographic coordinates.

Is there a projection I can generate that is fast to apply, even if its at the cost of somewhat not-so-accurate (~50m-100m) results?

I'm not using any database tools or large libs here, just C++ as I need a small, portable set of code. It would be handy if I could come up with the projection on the fly based on the geographic coordinate polygons rather than having to refer to a list of 'proper' SRSs.


texts:
Our team is trying to do a Line of Sight analysis and we couldn't figure out what when we run it, the resulting lines drape along the surface until we read the How Ling of Sight Works

This question is similar to 3D Analyst Line of Sight: what am I doing wrong? and was never answered.

We have a TIFF DEM, we have points that have Z values, and the construct sight lines are similar to what we need. However, how can we have our TIFF DEM act as our Multipatch? There are no obstructions in our study other than the terrain itself.


texts:
I have a problem with ArcSde and FME. I try to import some line features to ArcSDE Oracle. When I start the script in FME, I have some errors: impossible to insert linestring in ArcSde because the scale is very small. The length of linestring is equal at 100 meters.

But I have the others geometry with the same length and there is no problem to insert theses features in ArcSDE.
So, i have no idea on this problem.
Is it possible to clean geometry with FME ? or a correction of geometry ?
(I try Geometry Filter, OGC validator, but nothing).
Thanks

Error:
*Line Feature not allowed in table 'CHEMN'. Migration attempted
Error while generating simple line shape for Table 'CHEMIN'. Scale Factor may be too small.:SDE_ERROR_CODE(-155) Linestring or poly boundary is self-intersecting
Error while setting the shape to 'sde30_nil' for feature in table 'CHEMIN'. Shape Type Not Allowed - Change Table Def.:SDE_ERROR_CODE(-35) Shapes of this entity type not allowed in this layer*


texts:
I am installing grass71 on centos6 , i have the following message after running grass71

[bachir@nexus ~]$ grass71
Starting GRASS GIS...
g.gisenv: /lib64/libc.so.6: version `GLIBC_2.14' not found (required by /usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/lib/libgrass_gis.7.1.svn.so)
g.gisenv: /lib64/libc.so.6: version `GLIBC_2.14' not found (required by /usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/lib/libgrass_datetime.7.1.svn.so)
Traceback (most recent call last):
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/gui/wxpython/gis_set.py", line 31, in <module>
    from core import globalvar
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/gui/wxpython/core/globalvar.py", line 29, in <module>
    from core.debug import Debug
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/gui/wxpython/core/debug.py", line 72, in <module>
    Debug = DebugMsg()
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/gui/wxpython/core/debug.py", line 37, in __init__
    self.SetLevel()
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/gui/wxpython/core/debug.py", line 43, in SetLevel
    self.debuglevel = int(grass.gisenv().get('WX_DEBUG', 0))
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 952, in gisenv
    s = read_command("g.gisenv", flags='n')
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 458, in read_command
    process = pipe_command(*args, **kwargs)
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 433, in pipe_command
    return start_command(*args, **kwargs)
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 372, in start_command
    if debug_level() > 0:
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 1535, in debug_level
    _debug_level = int(gisenv().get('DEBUG', 0))
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 952, in gisenv
    s = read_command("g.gisenv", flags='n')
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 461, in read_command
    return handle_errors(returncode, stdout, args, kwargs)
  File "/usr/local/grass7.1.svn-x86_64-pc-linux-gnu-08_02_2016/etc/python/grass/script/core.py", line 329, in handle_errors
    returncode=returncode)
grass.exceptions.CalledModuleError: Module run None [b'g.gisenv', b'-n'] ended with error
Process ended with non-zero return code 1. See errors in the (error) output.
ERROR: Error in GUI startup. See messages above (if any) and if necessary, please report this error to the GRASS developers.
On systems with package manager, make sure you have the right GUI package, probably named grass-gui, installed.
To run GRASS GIS in text mode use the -text flag.
Exiting...
[bachir@nexus ~]$ uname -a
Linux nexus 2.6.32-573.12.1.el6.x86_64 #1 SMP Tue Dec 15 21:19:08 UTC 2015 x86_64 x86_64 x86_64 GNU/Linux


texts:
In one project, I am trying to adopt marvellous styles for flow maps created by @underdark.
At this moment I have the following output

Shapefile was created with the usage of Flow maps (Oursins) Plugin.
Lines have different directions, i.e. the starting point of first is the ending point of the second and vice versa. So basically the table includes lines with both directions but without any description of those directions.

However, I faced several issues
1. Varying buffers' sizes
Buffers around points shall inherit parameters from a different layer, i.e. they are not the same size.
How to make them respond to the size of the circles' radiuses? Is there any kind of scaling factor, so to make my buffers "dynamic"?
difference(
  difference(
    $geometry,
    buffer(start_point($geometry), 500)
  ),
  buffer(end_point($geometry), 500)
)

Maybe I to apply layer_property() or dbvaluebyid() or represent_value() if I want to extract the dynamic value from a different layer?
I know that the buffer values in the expression above only produce appropriate results for line datasets which use a CRS in meters and will have to be adjusted for other units.
2. Proper allocation of arrows
I am using an expression for Arrow type > Plain
CASE
    WHEN x(start_point($geometry)) - x(end_point($geometry)) < 0 THEN 1
    ELSE 2
END

but still, I can not properly allocate orange halves of Arrows on the opposite side of blues.
3. The width of the arrows
I want to adjust the line width accordingly to the values of a field "FLUX". For that reason, I am using an expression for Arrow width.
scale_linear("FLUX",0,10,0.1,3) * CASE WHEN  "FLUX" > 1000 THEN 1 ELSE 0 END

How shall I handle those issues?

My desired output should look like this, curved lines can be ignored.

References:

QGIS-resource-collections
Seeking QGIS geometry generator manual?
Details of good flow maps
New style: flow map arrows
Representing network flows?
Selecting features with same geometry in one shapefile


texts: What's wrong with the resolution of rasters?

texts:
I have two world area rasters. The metadata showed as follows:

file1:

Driver
GDAL provider
ENVI
ENVI .hdr Labelled
Dataset Description
D:/temp/2000npp/MOD17A3_Science_NPP_2000.tif
Band 1
Dimensions
X: 4500 Y: 1750 Bands: 1
X : 2250,Y 875
X : 1125,Y 438
X : 563,Y 219
X : 282,Y 110
Origin
-180,80
Pixel Size
0.08,-0.08
No Data Value
*NoDataValue not set*
Data Type
UInt16 - Sixteen bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
-180.0000000000000000,-60.0000472867500037 : 180.0001215945000013,80.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet

file2:

Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
D:/temp/2000npp/2000.tif
AREA_OR_POINT=Area
Band 1
STATISTICS_COVARIANCES=559.9966665827818
STATISTICS_MAXIMUM=164.06666564941
STATISTICS_MEAN=10.482447155964
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=23.664248701
Dimensions
X: 4500 Y: 2250 Bands: 1
X : 2250,Y 1125
X : 1125,Y 563
X : 563,Y 282
X : 282,Y 141
Origin
-180,90
Pixel Size
0.08,-0.08
No Data Value
-3.40282e+38
Data Type
Float32 - Thirty two bit floating point
Pyramid overviews
Layer Spatial Reference System
+proj=longlat +datum=WGS84 +no_defs
Layer Extent (layer original source projection)
-180.0000000000000000,-89.9999959766864777 :  179.9999919533729553,90.0000000000000000
Band
Band 1
Band No
1
No Stats
No stats collected yet

The rows number of file1 is 1750 while file2 is 2250 while both cols number is 4500. And the extent area looks like 80°N to 60°S on Geography showed in the map. But the file2  looks like stretched. I'm puzzle that:
If with pixel size 0.08° and 2250 rows, the map should contain artic pole area and antiartic( 90°N to 90°S ,for 2250*0.08=180). I think the file1 is right (1750*0.08=140 ,from 80°N to 60°S). What's wrong with file2? How to adjust it?

file1 showed in QGIS:

file2 showed in QGIS:

raster lay overlap will show clearly: