texts: Google Earth Engine: Error for mapping expression function over feature collection

texts:
I'm having an issue when trying to create a function to map over a feature collection.

I'm trying to create SST thresholds for several regions where over a certain SST value is 1, under another SST, value is 0 and in between is a linear function between those cut offs.

I tried this with a single image first and it looked like it had worked. Link here

However, I'd like to apply this to a featurecollection I've created , so I wanted to create a function of this and then map it over the collection. However I'm getting an error

Invalid argument specified for ee.List(): mean.

I think it is because I'm not defining mean correctly in the expression. I've tried lots of different options and can't figure it out. The code I'm using is here

var sst = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI').select('sst')
.filterDate(startDate, endDate)
.filterBounds(reefs_buffer)
.map(function(image){return image.clip(EEZ)}); //reefs_buffer

var sstMean = sst.reduce(ee.Reducer.mean())

print('sstMean', sstMean)

var reefs_sstMean = sstMean.reduceRegions({
  reducer: ee.Reducer.mean(),
  collection: reefs_buffer,
  scale: 30});

 print('reefs_sstMean', reefs_sstMean)

function sstHS_function(feature) {
   var sstReefs = feature.select('mean');
   var sstHS = sstReefs.expression(
"(b('mean') > 29.5) ? 1" +
 ": (b('mean')  < 24) ? 0" +
 ": 2/11*mean - 48/11 ",
 {'mean' : mean});
   return ee.Feature(sstHS);
 }

 var sstHS_collection = reefs_sstMean.map(sstHS_function)
 print(sstHS_collection)
 reefs_sstMean.map(sstHS_function)


texts: Converting rectangle bounding box from N, S, E, W decimal to a topmost NE, and bottom most SW rectangle

texts:
I want to combine two different database applications with a search, and in order to do that I need convert a rectangle box from its N, S, E, W decimal coordinates to a a box that is amde up from a top NE and bottom SW point.

Any advice on what calculation would help me get started?


texts: Alias for postgis function is really slow

texts:
I have created an alias for st_expand function like this

CREATE OR REPLACE FUNCTION expand(
geom geometry,
value double precision)
RETURNS geometry AS
$BODY$
    DECLARE my_geom geometry;
BEGIN
    my_geom =  st_expand(geom,value);
    return my_geom;
END;
 $BODY$
LANGUAGE plpgsql VOLATILE;

This function on running my sql query returns in 1.5s, but if I directly use st_expand function - the return is in ~ 30ms.

Do note, I have to do an alias, since RDS Postgres doesn't allow me to create expand same as st_expand since they restrict creation of functions of language C.

The time to change all our sql scripts(legacy) to call st_expand and not expand is way too long. Hence creating the alias, changing to st_expand is not an option for now.

Is there any change that I can do to ensure the time is still ~30ms and not 1.5s since there's no logic done here.


texts: Unable to connect PostgreSQL to ArcMap 10.1

texts:
I'm trying to connect my local PostgreSQL database to ArcMap 10.1 but for some reason I'm unable to connect to the server. Any ideas of what might be causing the issue. For the instance option in ArcMap I'm using the name of my server which is PostgreSQL 9.3 but i get "the unable to connect to data server" error. I have also tried localhost:5432 but i get the same error. Attached there is a pic showing the error the inputs and PostgreSQL server properties.
Thanks for the help


texts: Polygons from Line layer grouped by value

texts:
Is there a method to create polygons from a line layer, grouped by data values?

I have a line shapefile layer with +/-1000 records, about half of the records have a unique group value that designates them by a geographic area.  I'm currently symbolizing the lines based on the group that they belong to, but I'd also like to be able to create a layer with a polygon for each group.


texts: using where clause for selecting more than one value in the field in arcpy.Select_analysis

texts:
I have a problem in using the Select_analysis or Select by attribute, what I need t do, is to extract all the residential uses from a land use feature class, and make a feature class for residential, but there are more than one value for residential uses like (''Residential and driving school'', ''mixed commercial Residential'', ''mixed residential with other uses'', 'residential')
 As I use the following script:

uses_FieldForLandUses= arcpy.GetParameterAsText(4) # field of land use
Values_of_Residential_uses_in_Land_use_field= arcpy.GetParameterAsText(5) # type: string (multi value)- all the values for residential
Output_Residential_uses_parcels= arcpy.GetParameterAsText(6) # type: feature class
#Variables for residential uses
valueList = [x.strip() for x in Values_of_Residential_uses_in_Land_use_field.split(";")]

FieldInfo = arcpy.ListFields(inputLandUseParcels_FeatureClass,uses_FieldForLandUses)[0]
fieldType = FieldInfo.type

if uses_FieldForLandUses  == '#' or not uses_FieldForLandUses:
    Field = "LAND_USE" # provide a default value if unspecified

# Process: Select

if fieldType == 'String':
    where = "{} IN ('{}')".format(arcpy.AddFieldDelimiters(inputLandUseParcels_FeatureClass,uses_FieldForLandUses), "', '".join(valueList))

arcpy.Select_analysis(inputLandUseParcels_FeatureClass, Output_Residential_uses_parcels, where)

but when selecting more than one value for residential, I get the following error:

ExecuteError: ERROR 000358: Invalid expression land_use IN (''Residential and driving school'', ''mixed commercial Residential'', ''mixed residential with other uses'', 'residential')
Failed to execute (Select).


texts:
I am working with a Mac OSX El Capitan, QGIS 3.4.5 Madeira and Python 3.6.1.

I have a question which I believe is essentially a python one but applied into QGIS.

I intend to program some algorithms from the Saga tools to run various files. So, for instance I want to clip a raster with a polygon and so I run the following:

processing.run("saga:cliprasterwithpolygon", { 'INPUT' : '/Users/.../wheat_HarvestedAreaHectares.tif',
'OUTPUT' : '/Users/.../wheat.sdat',
'POLYGONS' : '/Users/.../africa_countries.shp' })

The raster files are agricultural data from all the world and I am clipping it to a shape file of the African continent to do some computations afterwards with zonal statistics.

Naturally, I would love to find a way to run this code for several crops. I am imagining something as creating a list with the names of the crops and looping through the code with the name of the crop as a macro.

Something like:

crop = [maize, wheat, rice]
for crop in crops:
    processing.run("saga:cliprasterwithpolygon", { 'INPUT' : '/Users/.../`crop'_HarvestedAreaHectares.tif',
        'OUTPUT' : '/Users/.../`crop'.sdat',
        'POLYGONS' : '/Users/.../africa_countries.shp' })

However, I've read elsewhere that Python does not work with macros, and similar questions around stackexchange does not seem to address this issue directly.

Any hint?


texts: How to calculate the average polygon shape in R (adehabitatHR)?

texts:
I have a code that generates polygons from a MCMC model. For each iteration in the model, I get a polygon of the 95% credible interval for the area. The idea here is that if I generate 10 polygons, I'm able to compute the area for each and then calculate the average area for all polygons.
So,

is it possible to compute the average polygon from all the polygons generated?
Would this "average polygon" be of equivalent size to the average area of all polygons?

Here is the code:
# 1. Packages
library(adehabitatHR)
par(mfrow=c(1,1))

kernel.poly.all = NULL
# 2. Assigning values
set.seed(13456789)
for (potatoe in 1:10) {
  # 3. Empty Dataframe
  points <- data.frame(ID = double())
  XY_cor <- data.frame(X = double(),
                       Y = double())

for(i in c(1:100)){

  if(i >= 50){points[i, 1] <- 1}
  else {points[i, 1] <- 2}
  XY_cor[i, 1] <- runif(1, 0, 100)
  XY_cor[i, 2] <- runif(1, 0, 100)}

# 4. Transform to SpatialDataframe
coordinates(points) <- XY_cor[, c("X", "Y")]
class(points)

# 5. Domain
x <- seq(-50, 150, by=1.) # resolution is the pixel size you desire
y <- seq(-50, 150, by=1.)
xy <- expand.grid(x=x,y=y)
coordinates(xy) <- ~x+y
gridded(xy) <- TRUE
class(xy)

# 6. Kernel Density
locs = SpatialPoints(points)
kud_points <- kernelUD(locs, h = "href", grid = xy)
# image(kud_points)

kernel.ref <- kernelUD(xy = locs,
                       # grid = xy  ,
                       grid = 250,
                       extent = 10,
                       # same4all = TRUE,
                       h = "href")
area.save = kernel.area(x = kernel.ref,
                           unin = 'm', unout='m2',
                           percent = 95)

kernel.poly <- getverticeshr(kernel.ref,
                             percent = 95,
                             unin = 'm', unout='m2')
kernel.poly.all = c(kernel.poly.all,kernel.poly)
if(potatoe==1){
  plot(kernel.poly,
       border = "black",
       lwd=1,
       axes = TRUE,
       xlim = range(x),
       ylim = range(y))
}
if(potatoe>1){
  plot(kernel.poly,
       border = "black",
       lwd=1,
       add = TRUE)
}
print(potatoe)
}

The object kernel.poly.all contains the 10 iterations of loop, so a total of 10 polygons.

The code above was derived from this question and adapted to make a reproducible example.


texts:
My task is to calculate the average ground survey level for each tree in my dataset, using the 3 closest survey points (optionally within a 3m radius of the tree). I want to write this in a Spatial VIEW (or a table with triggers) so that any time I add a new tree (or new survey points, though that's unlikely) the values update themselves.

Here is the core of my VirtualKNN() CREATE statement to show what I'm getting at, minus the 3m radius:

CREATE view tree_avglev_k AS
SELECT
   t.geometry AS geometry,
   t.pk_uid AS pk_uid,
   t.treeid AS treeid,
   AVG(l.lvl) AS avg_lvl
FROM
   (SELECT
         k.*,
         t.pk_uid
      FROM
         knn k,
         trees t
      WHERE
         f_table_name = 'levels'
         AND ref_geometry = t.geometry)
   a
   LEFT JOIN
      levels l
      ON (l.ogc_fid = a.fid)
   LEFT JOIN
      trees t
      ON (t.pk_uid = a.pk_uid)
GROUP BY
   a.pk_uid
ORDER BY
   t.treeid

But as of Jan 2018 QGIS 2.99/3 still does not support the VirtualKNN() function in Spatialite 4.4.0, so it can't process the information (see feature request)

Is there some other more efficient way to work around this (without using PostGIS/any other client-server database)?

What I have tried (as reference):

I have a long, ugly alternative using PtsDistWithin() and Virtual Index. It's not the most efficient (still new to SQL) but on ~500 trees and ~5000 survey points, it loads in about 2 seconds in spatialite-gui (bit slower than the VirtualKNN() query)

But when loading it in QGIS it is particularly slow - and buggy (despite selecting unique integer column in DB Manager, it loads the geometry but won't load attributes). When I add other constraints (like using only levels within a certain area, again using VirtualIndex) it's practically unloadable.

Maybe I need a different approach (e.g., create normal view with just data, join to geometry separately)

However the main thrust of the question is not so much to fix my alternative approach, but to ideally find a solution as concise and zippy as VirtualKNN() for QGIS.


texts: Download geo-referenced photos from ESRI collector app

texts:
I have ESRI Collector app on an iPad. I have used it to mark up a number of features on site, complete with photos and videos.

My end goal is a kmz of georeferenced features complete with photos in Google Earth (the videos I'll just have to provide a link or similar).

I am stuck at the first hurdle and cannot work out how to export the features with photos into either arcmap or Google Earth. I can export the collector feature class as a shapefile or various other formats but none have the photos linked as far as I can tell.

How can I accomplish this download of georeferenced photos?


texts: Finding Polygons which do not touch any other polygon in a single oracle database table

texts:
I currently have an oracle table of approximately 100,000 records - all of which are polygons.
I want to find all the polygons which do not touch any other polygon within the same table.

There are a few previous posts which use the ST_TOUCHES query but this relies on you having two geometries to compare against from different tables. Unless my understanding is incorrect?
Where as what I want is a query which looks at the one table and one Geometry column and looks to see if a polygon touches another.


texts: What is the correct way to assign GeoPandas a CRS from Rasterio?

texts:
I'd like to be able to create a new GeoPandas object and inherit the CRS from Rasterio, and save a shapefile. The current projection format is seen as invalid when writing an ESRI driver.

From here, it would suggest GeoPandas requires CRS in dict format:
https://geopandas.org/projections.html and Reprojecting data with GeoPandas using 'to_crs'?, though this hints as dataset.crs.data (How can I superimpose a geopandas dataframe on a raster plot?).

It seems to be known that you must pass the crs_wkt directly to geopandas.to_file method:
GeoPandas to_file() saves GeoDataFrame without coordinate system

simple psuedo code

boxes = geopandas.GeoDataFrame()
with rasterio.open(raster_path) as dataset:
    boxes.crs = dataset.crs.to_dict()

boxes.to_file(<filename>,driver="ESRI Shapefile",crs_wkt=boxes.crs)

yields

TypeError: invalid crs_wkt

The specific projection i'm trying is:

dataset.crs.to_dict()

{'proj': 'lcc',
 'lat_1': 43,
 'lat_2': 45.5,
 'lat_0': 41.75,
 'lon_0': -120.5,
 'x_0': 399999.9999999999,
 'y_0': 0,
 'ellps': 'GRS80',
 'units': 'ft',
 'no_defs': True}

Though I've played around with others, so I think its a workflow issue.

Rasterio version: 1.1.2
GeoPandas version: 0.6.2
Fiona version: 1.8.13

As an added comment, opening the raster in R and extracting the proj 4 string

> proj4string(a)
[1] "+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999999999 +y_0=0 +ellps=GRS80 +units=ft +no_defs"

and assigning the string into boxes.crs above yields no error.


texts: How to fit index report into ArcMap layout?

texts: I want to fit report into map. I'm creating grid and in order to index it I'm creating report. View > Create Report. It gives me a new window that includes texts. I want to fit texts for index points and ways. How can I fit index into my map ?


texts: GeoServer container need to test locally?

texts:
GeoServer container need to test locally

I did some tests with geoserver and openlayers but it did not work with the layers I have on my machine. I did not need the container to use locally calling my WMS layers using openlayers. My javascript pages are in the openlayers folder on my C :/ are not in JSF project.
My question is I need the container to use locally in my tests?

    <!DOCTYPE html>
<html lang='en'>
<head>
 <meta charset='utf-8' />
 <title>My OpenLayers Map</title>
 <script type='text/javascript' src='OpenLayers.js'></script>

 <script src="http://maps.google.com/maps/api/js?sensor=false&v=3.2"></script>

 <script type='text/javascript'>

 var map;
 var mercatorProjection = new OpenLayers.Projection('EPSG:900913');
 var latLongProjection = new OpenLayers.Projection('EPSG:4326');
 var bounds = new OpenLayers.Bounds(
                -51.2809219360352, -30.2445888519287,
                -51.0207977294922, -29.9661273956299
            );

 function init() {

  map = new OpenLayers.Map('map_element',{
    maxExtent: bounds,
    maxResolution:  0.0010877400636672,
    units: 'm',
    //allOverlays: true,
    projection: mercatorProjection,
    displayProjection: latLongProjection
 });

 var google_streets = new OpenLayers.Layer.Google(
    "Ruas",
    {numZoomLevels: 20}
    );

var wms_layer = new OpenLayers.Layer.WMS(
'Linhas Ônibus',
'http://localhost:8080/geoserver/Teste/wms',
{layers: 'Teste:onibus', transparent: true},
{isBaseLayer: false,
    opacity: 0.7}
);

//Adiciona as camadas ao mapa
map.addLayers([google_streets, wms_layer]);

var point = new OpenLayers.LonLat(-51.22,-30.08);
point.transform(new OpenLayers.Projection("EPSG:4326"),
    map.getProjectionObject());
map.setCenter(point, 10);
 //Camada de controle que vai mostrar as camadas no mapa
 map.addControl(new OpenLayers.Control.LayerSwitcher({}));

 //Verifica se o mapa tem um ponto central e o extende a sua extensão máxima
 if(!map.getCenter()){
 map.zoomToMaxExtent();

 }
 }

 </script>
</head>

<body onload='init();'> <!-- Chama a função js init() -->

<!--Elemento HTML onde o mapa é exibido -->
 <div id='map_element' style='width: 800px; height: 800px;'>
 </div>

</body>
</html>


texts: I am using python 2.7 now in want mapnik to be included in python. I have use exe file to include other library file like GDAL, and other. But for Mapnik can any one suggest me any exe like that so it can be easy or library file but how to install it, as i found many help not usefulness. PLEASE  any one help me out in this issue


texts: LonLat Points to Grid of Polygons in PostGIS

texts:
I'm going to preface this with it's my first time using PostGIS.

I'm starting out with a 5 km x 5 km grid cell in which I know the center point lon/lat (SRID=4326). I'm trying to create a polygon for each of these grid cells in PostGIS using this command:

    create table lonlat_poly as
    select ST_MakeEnvelope(lon-2.5 km, lat-2.5 km, lon+2.5 km, lat+2.5 km)
    as newlonlat from lonlat;

where lonlat contains the two columns lon and lat.

The potential pitfall is that the simple subtraction and addition of 2.5 km is not going to be correct due to projection issues (right?). I'm sure there's a more precise way of doing this, but I haven't thought or stumbled upon anything yet. Thanks in advance for any hints you can provide!


texts: How to cut data inside GeoJSON polygon from adjacent ASCII grid files?

texts:
I have a set of Digital Terrain Model ASCII grid files; they're the output of a LIDAR survey over our city, and so cover adjacent areas. Now I need to "cut" a polygon out of these files; this polygon is in GeoJSON format and can overlap across more than one of these DTM files.

Then the cut area should be put into a new ASCII grid file, because it needs to be passed to an existing program our team created, with the area outside the polygon.

Here is a diagram to give a more visual description.

How can I execute this operation?

I have GDAL installed, and I have admin rights, so I can install any software I need. Of course it would be better if it was Open Source, but the most important things is that I will have to automate this operation next, so something that works from the command line or has a Python API would be ideal for me.


texts: Feature selections not being honored when exporting subset of polygon .shp (ArcPy)

texts:
I want to write a for-loop where I select random subsets of a polygon shapefile (10% of the 600 total polygons) and save each subset to a new shapefile, named subset1, subset2, etc.

My code successfully runs and outputs shapefiles with the correct name; however, each output shapefile is identical and contains all of the polygons in the original shapefile, not the subset selected prior. It doesn't seem to honor my feature selections the way that manually right-clicking Data > Export Data > Export (Selected features) does.

#define variables
Polygons600 = r"C:\DataFolder\600Polygons.shp"
out_path = r"C:\OutputFolder"

# defines subset selection tool, ignoring any current selections in ArcMap
def SelectRandomByPercent (layer, percent):
    #layer variable is the layer name in ToC
    #percent is percent as whole number  (0-100)
    if percent > 100:
        print "percent is greater than 100"
        return
    if percent < 0:
        print "percent is less than zero"
        return
    import random
    fc = arcpy.Describe (layer).catalogPath
    featureCount = float (arcpy.GetCount_management (fc).getOutput (0))
    count = int (featureCount * float (percent) / float (100))
    if not count:
        arcpy.SelectLayerByAttribute_management (layer, "CLEAR_SELECTION")
        return
    oids = [oid for oid, in arcpy.da.SearchCursor (fc, "OID@")]
    oidFldName = arcpy.Describe (layer).OIDFieldName
    delimOidFld = arcpy.AddFieldDelimiters (layer, oidFldName)
    randOids = random.sample (oids, count)
    oidsStr = ", ".join (map (str, randOids))
    sql = "{0} IN ({1})".format (delimOidFld, oidsStr)
    arcpy.SelectLayerByAttribute_management (layer, "", sql)

for i in range(1,100):
    SelectRandomByPercent("My Layer Title", 10)
    arcpy.FeatureClassToFeatureClass_conversion(Polygons600, out_path, "Polygons600_subset" + str(i) + ".shp", "", "", "")

Other questions/answers on this site recommend using FeatureClassToFeatureClass_conversion or CopyFeatures_management to subset shapefiles, claiming they honor feature selections when they export. I have tried both in my loop (and I've also tried creating a .lyr version of my input .shp and using that as the selection layer) to no avail.

Does this problem have to do with how I've selected my subset polygons (the SelectRandomByPercent tool)? For instance, if I select a subset of my shapefile using a SQL statement in the FeatureClassToFeatureClass argument it exports just fine. How can I make FeatureClassToFeatureClass or CopyFeatures recognize my selections?


texts: Why python function not work in custom function in QGIS 3.2?

texts:
in QGIS atlas I use a custom function to find the first file containing texts :

import os, fnmatch
from PyQt5.QtCore import QFileInfo
from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def findphoto(nb, src1,src2,feature, parent):
    lnb=0
    val='NA'
    filen = QFileInfo(QgsProject.instance().fileName())
    Dir_report=filen.absolutePath()
    for root, directories, filenames in os.walk(Dir_report):
        for filename in filenames:
            if fnmatch.fnmatch(filename, '*' +str(src1) + '*' + str(src2) + '*.jpg' ):
                lnb=lnb+1
                if lnb==nb:
                    val=os.path.join(root,filename)
                    break
    return val

In QGIS 2.18 it works fine (just change PyQt5 by PyQt4), in QGIS 3.2 it return "NA" value.

In QGIS 3.2 if I try in python mode, it works fine:

import os, fnmatch
from PyQt5.QtCore import QFileInfo
from qgis.core import *
from qgis.gui import *

def findphoto(nb, src1,src2):
    lnb=0
    val='NA'
    filen = QFileInfo(QgsProject.instance().fileName())
    Dir_report=filen.absolutePath()
    for root, directories, filenames in os.walk(Dir_report):
        for filename in filenames:
            if fnmatch.fnmatch(filename, '*' +str(src1) + '*' + str(src2) + '*.jpg' ):
                lnb=lnb+1
                if lnb==nb:
                    val=os.path.join(root,filename)
                    break
    return val

why?


texts: How can I identify the pixel with the maximum value in a grid?

texts:
Possible Duplicate:
  How to find the locations of the highest values in a raster?

How can I identify the pixel with the maximum value in a grid and get its coordinates? Any VBA code please?


texts: How to update CRS r object from PROJ4 to PROJ6

texts:
I have a script which contains several r objects which are CRSs that I use for transforming spatial objects. After updating R and many packages, there is a warning which wasn't there before.
LCC <- CRS("+proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs")

Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition

I understand that this is due to a change in how spatial information is stored... PROJ4 to PROJ6 upgrade and "Discarded datum" warnings
..but beyond that, I don't understand what this error means or if I need to do anything in response.
Any thoughts?


texts: How to calculate percentage area cover using Focal Statistics

texts:
I am trying to produce raster outputs in ArcGIS for later use in MaxEnt to model species distributions, and want to produce variables with statistics on the percentage cover of different habitat polygons (i.e. woodland) and linear features (i.e treelines or hedgerows) at different spatial scales (e.g within a moving analysis window).

I have a fairly rudimentary understanding of ArcGIS. I feel that Focal Statistics tool should be able to achieve what i want if i convert my layer inputs to raster format, however percent is not listed as an option in the Statistics field.


texts: How to change the GeoServer port number in the "Mac OS X installer" version

texts:
There is a version of GeoServer available at http://docs.geoserver.org/latest/en/user/installation/osx_installer.html which is very simple to install on a Mac, only requiring the application to be dragged to the Applications folder and opened.

The only options within the GeoServer application itself are Start/Stop, and GeoServer defaults to using port 8080.

Is it possible to change the port number that is used by GeoServer in this version?


texts: How to burn values from a shapefile into raster with gdal.Rasterize?

texts:
I have this code to rasterise a point shapefile into raster TIFF.
vector_ds  = ogr.Open('stations.shp')
shp_layer = vector_ds.GetLayer()
shape_file = 'stations.shp'
output_raster = os.getcwd()+'/rasterized_points.tif'
pixel_size = 0.01
xmin, xmax, ymin, ymax = shp_layer.GetExtent()
ds = gdal.Rasterize(output_raster, shape_file, xRes=pixel_size, yRes=pixel_size,
                    burnValues='data',outputBounds=[xmin, ymin, xmax, ymax],
                    outputType=gdal.GDT_Byte)
ds = None

It works, but the result doesn't have values from 'data' column. Parameter 'attribute' of RasterizeOptions doesn't help neither. Can somebody explain, how I should make the Rasterize burn the right values (0,1,2,..) from shape.


texts: Error when trying to convert dxf to shapefile?

texts: I m using QGIS 2.18. I am getting the following error when I try to convert dxf to shapefile. The other software like online convertor and arcgis give no errors in converting the same shapefile. The program crashes.


texts: Reading pre-10.0 File Geodatabase (GDB) using OGR?

texts:
I have been tasked with processing USGS's National Hydrography Dataset. However, it appears that the only available format is FileGDB (version 9.2-9.3.1). Looking at OGRs support for FileGDB, it looks like they only support those created by ArcGIS 10 and above. Does anybody have any ideas how I would be able to view and export this data to a more usable format?

If anybody is curious about the data, it can be directly downloaded through here.


texts:
Using python I want to convert coordinates from one spatial reference system to another. I usually use the osgeo library for this which works very well. I have however come to the situation that I need to convert to a UTM zone further away from the equator (say UTM zone 29U for example). For this zone I can not find an EPSG code on spatialreference.org.

Why is it that I can only find codes for zone either north or south of the equator (29N and 29S)?


texts: Using Python Script and Model Builder to create JSON output

texts:
So I have a script that creates a selection based off of a user's entry and stores select attributes of that entry into a python dictionary.

Ideally, the user would give a street segment ID  and the code would return the IDs of the nearest features within a certain distance with specific attributes, and store all of this into a dictionary. This dictionary is then converted into JSON using (arcpy.SetParameterAsText(1, json.dumps()).

The goal is to publish this as a processioning service. The code works and so I tried to recreate this via model builder in ArcToolBox. This also runs too but when uploaded as a service, it fails.

Is there anyone that can give insight to why this is happening? I am assuming that I may need to re work my model.

Attached is the resulting Server messages. The chosen ID is 5705. Which again, works perfectly when just running the script and when running the script as a toolbox on ArcMap. However on the server it fails.


texts: Service areas based on time to travel raster layer in ArcGIS

texts:
I have a series of raster layers which contain values for the time to travel (in minutes by motor) to a series of health facilities. I also have the GPS coordinates of each of the healthcare facilities. I also have a raster layer which contains the values for the travel time to the nearest (minimum time to travel) healthcare facility.
I need to create a travel time based "catchment area" or "service area" around each health facility. Ideally, every pixel within each catchment area will be assigned to the healthcare facility it takes the least time to reach, and the catchment areas will not overlap, see image below for what I would like to produce

I have looked at the suggestion to a similar query here, and tried the Path Allocation Distance tool:
Input source data: My GPS coordinates for healthcare facilites
(the "source field" the defaults to my X coordinates, with no option to use the heathcare facility ID
Input cost raster: the time-to raster to nearest healthcare facility.
However the areas overlap as, see image below.
What are the best tools and methods to use in Arcmap 10.8.1?


texts:
There are a few similar Questions which mention this issue but sadly not helping me out with my problem.

I want to generate an output out of my tool but am not sure which phrase and Datatype to use.

Additionally while importing the Script to ArcMap it just says failed to execute without further information. In Pycharm its stopping at the processing of merge and I am not sure where the mistakes are.

# Import arcpy module
import arcpy
from arcpy import env
arcpy.env.workspace = arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput = 1
arcpy.CheckOutExtension("Network")

# Local variables:
Streets_ND = "Streets_ND"
Route = "Route"
Route__2_ = Route
Stores_to_merge = "Stores_to_merge"
MainStore = "MainStore"
Stores = "Stores"
Stores_select_shp = "Stores_select"
Route__3_ = Route__2_
Solve_Succeeded = "false"
Route__4_ = Route__3_
Route_Directions_html = "Directions.lyr"

# Process: Make Route Layer
arcpy.MakeRouteLayer_na(Streets_ND, "Route", "Minutes", "USE_INPUT_ORDER", "PRESERVE_BOTH", "USE_TIMEWINDOWS", "", "ALLOW_UTURNS", "Oneway;RestrictedTurns", "USE_HIERARCHY", "", "TRUE_LINES_WITH_MEASURES", "")

# Process: Select_MainStore
arcpy.Select_analysis(Stores_to_merge, MainStore, "OBJECTID = 1")

# Process: Merge MainStore+Stores_to_merge
arcpy.Merge_management([Stores_to_merge, MainStore], Stores)

# Process: Select_createPick-List
multival = arcpy.GetParameterAsText(1)  # 'Store_1;Store_2;Store_3'
vals = multival.split(";")

where_stub = """"NAME" IN ({})"""
Stores = ','.join(["'{}'".format(x) for x in vals])
where = where_stub.format(Stores)

arcpy.Select_analysis(Stores, Stores_select_shp, where)

# Process: Add Locations
arcpy.AddLocations_na(Route, "Stops", Stores_select_shp, "Name Name #", "5000 Meters", "", "Streets SHAPE;Streets_ND_Junctions NONE", "MATCH_TO_CLOSEST", "APPEND", "NO_SNAP", "5 Meters", "INCLUDE", "Streets #;Streets_ND_Junctions #")

# Process: Solve
arcpy.Solve_na(Route__2_, "SKIP", "TERMINATE", "", "")

# Process: Directions
arcpy.Directions_na(Route__3_, "HTML", Route_Directions_html, "Miles", "REPORT_TIME", "Minutes", "en", "NA Desktop", "")

arcpy.SaveToLayerFile_management(Route_Directions_html, "Directions.lyr")


texts: Selecting by location AND attribute using one query

texts:
I am using ArcMap 10.8 and ArcGIS Pro 2.7
I have two feature classes: one containing currently uploaded county address points (cad_n) and the other containing updated county address points (d_n). A large majority of these points are the same, but slightly offset. In order to create a single feature class containing a combined, full set of data, I am trying to delete overlapping points from the cad_n file.
I have a Python code that first selects features from cad_n that are within 500 feet of d_n features, and then selects a subset of those features that share an address with d_n:
import arcpy

arcpy.env.workspace = "C:/Users/carson.jones/Desktop/Data/DRIZCAD.gdb"

locq = arcpy.SelectLayerByLocation_management('cad_n', 'WITHIN_A_DISTANCE', 'd_n', '500 Feet')

locattq = arcpy.SelectLayerByAttribute_management(locq, 'SUBSET_SELECTION', 'T_Address IN (SELECT T_Address FROM d_n)')

While this is close to what I want, I would rather make a selection of cad_n features within a distance of a d_n feature that has an identical address.
Is there a way to combine these queries into one selection?


texts: Creating FeatureLayer from CIMFeaturelayer in ArcGIS Pro?

texts:
I've been asked to develop a strategy to work with Layers that have not yet been added to a Map.

From what I can tell, one cannot simply create a new layer as one can with ArcObjects. One must use a method on the LayerFactory, all of which add a layer to the map upon creation.

So instead of using Layers, I thought I would have my analysis engine deal with layer definitions, e.g. CIMFeatureLayer, which can exist outside of a map.

For example, my analysis engine would generate CIMFeatureLayers and then, once I'm ready to present the layer, I would use something similar to LayerFactory.CreateFeaturelayer to create a layer from the CIMFeatureLayer.

I've searched but cannot find any method that does this.

Clearly Pro does this behind the scenes - if I save a .lyrx file and look at it in a text editor, I see that it is a json serialization of a CIMLayerDocument with CIMFeatureLayers within it.  How can I do what Pro does when it reads the .lyrx file?

Update

I also tried this code, but SetDefinition doesn't seem to do anything.

private void Test()
{
    var def = MapView.Active.Map.Layers[0].GetDefinition() as CIMFeatureLayer;
    var fLayer = LayerFactory.CreateLayer(def.FeatureTable.DataConnection,MapView.Active.Map) as FeatureLayer;
    fLayer.SetName("my layer");
    // name remains "my layer" Setdefinition has no discernable effect
    fLayer.SetDefinition(def);
}


texts: Default Value for Field with Python Expression in Function Editor

texts:
I have a line shapefile in QGIS with several columns and the values for many of those columns could be automatically calculated and updated depending on certain other fiels.

Now, I came across the default value option in the layer properties and I tried to build a Python function that automatically generates the value, but for some reason I do not manage to automatically get the value from the field that will serve as the argument for the function.

So, after saving the function, I try the expression

defaultPrice("Type")

I just get an error, although this is how I would usually apply functions when I manually update a field (like lower("Type") or whatever).

So, I have the field "Type" which has the values Gardentrench, Asphalt, Brick Pavement, As Built, HDD, Crossing and Grass. I have another field "PricePerMe" where I assign the cost of the "Type" per Meter. Not even calculated, just the price per one meter.

This is the code that I wrote:

from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def defaultPrice(Type, feature, parent):
    if Type == 'Gardentrench':
        price =  12.5
    elif Type == 'Grass':
        price =  15
    elif Type == 'As Built':
        price = 1
    elif Type == 'Asphalt':
        price =  30
    elif Type == 'Brick Pavement':
        price =  100
    elif Type == 'HDD':
        price =  250
    elif Type == 'Crossing':
        price =  90
    else:
        price = 0
    return price

In the expression field, I would like to trigger this function by

defaultPrice("Type")

and update this value. In the attribute table when I update a field value or values of a column, it does display me down a correct preview value, but after running it, there are no changes in the "Type" column:

In the default value expression editor I get this:

What am I doing wrong?

Looked at some examples that I found, but couldn't really apply any of them to my problem, although I am sure it's very trivial.

Edit:
Quick note: I know how to use virtual fields and that would do also, but they are not saved in the shapefile, only in the project file and I would like to learn something new.


texts: Obtain distance of each raster pixels from anoth feature/ proximity tool failing

texts:
I have a raster with population densities for Vietnam. Each pixel contains 1 population count. For the purpose of my analysis, I'd like to find the distance of each of the pixels to a line feature, representing the coastline as shown in the picture.

I have read a similar question in an other thread: Calculating distance to points in QGIS and managed to rasterize my coastline (even though there is a big gap in the southern part of the coastline).

When I try to use the proximity tool I get the following message on the log: The following layers were not correctly generated. I tried to save in different areas my output map but the error message keeps on coming back.
What I need is just the distance of each pixel from that coastline, similarly to what the "distance to nearest hub" tool does with vector data.


texts: Populating the empty xy fields in point layer

texts:
I have a number of points that form a line, with xyz values for each 100th meter, aka Chainage Number (this is a railway tunnel, chainage number 2900 = 2900 meters from a train station).

But I want xyz values for every meter and have therefore (in excel) added 99 rows below each 100th Chainage Number, and linearly calculated the z values in the empty rows  based on the z values for every 100th meter. In the added table image the Chainage Numbers are in field called "Inbound_Ch", not sure why the xy values are "1" in what I call empty rows.

How can I calculate the xy values in the empty rows?
I use arcgis 10.3.


texts:
I have created a MultiListBox control which is populated with a list of all currently loaded TAB file names in MapInfo. The user selections are then read from the MultiListBox and added to another array (so I can use these in an INSERT command later on).

However the problem I am getting is when I select records from my MultiListBox. As an example I have selected 3 TAB files (3 records) from my MultiListBox but my array size is only 2. If I select 5 records from my list my array size is only 4 and so on. When I select just the 1 record from my list my array size is saying 0.

Here are some code snippets to help.

Dim i as integer
Dim j as integer

i = NumTables()

ReDim TableArray(i)

For j = 1 to i
TableArray(j) = TableInfo(j, TAB_INFO_NAME)
Next

Print "Table Array no. of records: " + UBound(TableArray)

' code to read the values selected in the MultiListBox
Dim iSel, j As Integer

j = 0
iSel = ReadControlValue(1)

Do Until iSel = 0
j = j + 1
Redim SelectionArray(j) ' re-size our array. Now = 1
SelectionArray(j) = iSel
iSel = ReadControlValue(1)

Print "Selection Array items selected: " + SelectionArray(j)

Loop

Print "MultiListBox no. of selected records: " + UBound(SelectionArray)


texts: Validating user input against a range using tool validator in arcpy script?

texts:
I have a script that takes user input. For the "Project ID" input field, the user can only enter a 4 character ID. The ID may include letters and numbers, but is typically a 4-digit number. I'm trying to use the script's ToolValidator to make sure the user enters only 4 characters before the script will run. I am using the correct parameter index #.

I've tried several different code variations, but this is an example of my last attempt. With the following code, no matter what value I enter, the "Project ID" field returns blank with no error message.

Also, the "Project ID" parameter is a STRING data type in the event we encounter a project ID that contains a letter and is a required field. I would normally just check the length of the input, but couldn't find a len() method to use. I figured I could test if the entry was less/greater than 10000 to obtain the 4 digit restriction.

How should I proceed?

import arcpy
class ToolValidator(object):
  """Class for validating a tool's parameter values and controlling
  the behavior of the tool's dialog."""

  def __init__(self):
    """Setup arcpy and the list of tool parameters."""
    self.params = arcpy.GetParameterInfo()

  def initializeParameters(self):
    """Refine the properties of a tool's parameters.  This method is
    called when the tool is opened."""
    return

  def updateParameters(self):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""
    if self.params[6].altered and not self.params[6].hasBeenValidated:
      if self.params[6].value >= '10000':
        self.params[6].value = ""
    return

  def updateMessages(self):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    if self.params[6].value >= '10000':
      self.params[6].setErrorMessage ("Project number can only have 4 positive digits. Please enter a 4-digit project number.")
    else:
      self.params[6].clearMessage()
    return


texts: QGIS import osm data via OSM-plugin and sqlite does not only import selected tags

texts:
I have the following problem: I want to use the OSM plugin to import data from a sqlite database that contains OSM data. In the 'Export OSM topolgy to spatialite' dialog I choose the database file and the polygon geometry. Then I click the 'Exported tags' button. This shows me the list of tags together with the count of features. In my case, I select admin_level, which has 267 features. Now I press ok. After the export, I would expect that only 267 features are inserted in the added layer. But when I check the feature count of the new layer it has about 52000 tags, which is pretty much every polygon. When I had a look at the attribute table the admin_level field has the value NULL for every feature.

When I check the 'Not Null' checkbox, the new layer is empty and has no features at all, so it seems that all values are NULL in some sense. I'm not sure to what attribute the 'not null' conditions applies. If I would know, I could cross-check it in the database.


texts: How to be multiple publishers on ArcGIS Server

texts:
I have and account named "sde" and "publisher" on my ArcGIS 10.5.1 geodatabase.

I have two account on ArcGIS Server named "siteadmin" and "publisher".

So I have published the REST services on Server using "sde" db user fom ArcMap. My MXD sources are using "sde" user informations. And ArcGIS server services data source uses "sde" user accounts.

When "publisher" user wants to publish his own MXD file from ArcMap, an error occured following:

data source must be registered with the server.

Because his db connection user informations are different.

In this stuation, will all publisher use "sde" user?


texts:
I have n objects and each object has m positions.

Every position is defined by x, y (in millimeters) and a (unix) timestamp.

{
  objId: 'a-b-c',
  x: 123,
  y: 567,
  date: 123456789,
}

Now I want to calculate how long an object didn't move.

My first algorithm always updated a list of "last" positions for every object if it got a different position. If the new position wasn't different, I increased the standing duration for that object.

Problem I encountered:

The system which tracs the objects (almost) never delivers the same coordinates even if the object didn't physically move.

x:123 -> x:120 -> x:124 -> x:123 -> ...

Should I just round to centimeters or decimenters or are there "standard" GIS-algorithms to solve this?


texts:
I am a bit new in p\Pthon programming in QGIS. Now, I am trying to develop a QGIS Plugin about Suitability Analysis using Ordered Weighted Averaging process. So in brief, I want to make a plugin to produce a single map overlayed from several maps/layers and also assigned with user-based preference weight.

The idea, the plugin will provide the user to input several layers, and the user can assign their preference weight numbers.

Now, I am only able to create the UI and have succeeded to manage how to input the data. But I am stuck what is the next step? How to design the weight for each layers, and how to perform the overlay analysis itself?


texts: Overlaying Aviation Airspace Boundaries in QGIS 3D Map View

texts:
I'm a new user of QGIS 3.12.3 Bucuresti (and GIS in general).
I have loaded OS Terrain 50, overlaid with a mix of Google Satellite and OpenStreetMap which produces a 3D Map view that I am very happy with.
However, on trying to load obstructions such as masts, church spires etc. represented by columns up to the altitude (above MSL/Geoid) and protection surfaces represented by polygons in 3D space the obstructions and surfaces just appear to be draped over the surface of the terrain mapping with no elevation.
I started with obstruction data in CSV format, latitude, longitude and elevation which produced circles on the map surface.  Am I missing something simple here please?


texts: When ESRI tool "Calculate Grid Convergence Angle (Cartography)" is used to calculate line bearing or azimuth, is the bearing referenced to True North or Grid North?


texts: How to add multiple locations to google map?

texts:
Well I am not a js developer, and I want to build a google map with multiple locations.
I have a styled google map which looks like this, and I like to add some locations.
After some research I found this code works

var locations = [
                ['Bondi Beach', -33.890542, 151.274856, 4],
                ['Coogee Beach', -33.923036, 151.259052, 5],
                ['Cronulla Beach', -34.028249, 151.157507, 3],
                ['Manly Beach', -33.80010128657071, 151.28747820854187, 2],
                ['Maroubra Beach', -33.950198, 151.259302, 1]
            ];

Then I added it to my javascript(below), but it just doesn't work.

//<![CDATA[
$(function(){
jQuery(document).ready(function () {
    var locations = [
                ['Bondi Beach', -33.890542, 151.274856, 4],
                ['Coogee Beach', -33.923036, 151.259052, 5],
                ['Cronulla Beach', -34.028249, 151.157507, 3],
                ['Manly Beach', -33.80010128657071, 151.28747820854187, 2],
                ['Maroubra Beach', -33.950198, 151.259302, 1]
            ];
    var map;
    var centerPosition = new google.maps.LatLng(35.335293,161.323242, 16.050432,33.815918);
    var style =
 [{"elementType": "labels","stylers": [{ "visibility": "off" }]},
  {"featureType": "landscape","stylers": [{ "color": "#222222" }]},
  {"featureType": "road","stylers": [{ "color": "#000000" }]},
  {"featureType": "poi","stylers": [{ "visibility": "off" }]},
  {"featureType": "water","stylers": [{ "color": "#000000" }]},
  {"featureType": "administrative.country","elementType": "geometry","stylers": [
  { "visibility": "on" }, { "color": "#2c2c2c" }]}]
    var options = {
      zoom: 2,
      center: centerPosition,
      mapTypeControl: false,
      scaleControl: false,
      zoomControl: false,
      navigationControl: false,
      streetViewControl: false,
    };
    map = new google.maps.Map($('#map')[0], options);
    map.setOptions({
        styles: style
    });
    var image = {
        url: 'http://icons.iconarchive.com/icons/glyphish/glyphish/16/07-map-marker-icon.png',
        origin: new google.maps.Point(0, 0),
        anchor: new google.maps.Point(12, 59)
    };
    var shadow = {
        url: 'https://dl.dropboxusercontent.com/u/814783/fiddle/shadow.png',
        origin: new google.maps.Point(0, 0),
        anchor: new google.maps.Point(-2, 36)
    };
    var marker = new google.maps.Marker({
        position: centerPosition,
        map: map,
        icon: image,
        shadow: shadow
    });
});
});//]]>

How can I add these locations to my map?


texts: IElement, extract Element from GroupElement nested in GroupElement

texts:
As per title.

I am trying to extract Element from nested GroupElement. I can not find any code sample for that. Please refer to one or point me in to right direction, how to do that.

It sounds simple but turned out to be nightmare. I cannot figure out how whole mechonism of Group in the Group works. Any advice?


texts: How to Query in a web map for distance between 2 points

texts:
I am new to web mapping and am looking out for a solution to my issue.
I have been working on creating a Web Map for Paying guest location from my workplace and am looking to find the distance to each PG locations. However I am not sure how to do this? Can anybody help?

I am using GeoServer, Postgresql and Leaflet.
To do the query is there any plugins available or do I have to write a Javascript code?


texts: How do I turn on auto increment function for a field that already has incrementing numbers in it?

texts:
In ArcMap, I've got an attribute table with two object ID fields. One is active and populates every time I make a new feature, and the other, which is older, is not. They do not match exactly because there were entries removed earlier in time, so the active field goes up to 925 and the older one goes up to 1042, but there are only 925 features. I want to reactivate this older field so that it populates like the other. I.e. next time I create a feature, Object_ID will be 926 and the older will be 1043.

I know the Esri autoincrement() Python script to start a sequential numeric list, but I'm not looking to start a whole new one. I just want to turn on one that already exists so that it adds 1 to the value with each new entry.


texts: Deleting multiple polygons at once in ArcGIS 10.1?

texts: I have a cloud mask shapefile that accidentally identified buildings as clouds.  I want to delete this mass of polygons in Arcmap 10.1 by highlighting all of them and pressing delete.  Then re-saving the shapefile.  Right now I can only highlight one at a time through clicking and pressing delete.  It is taking a long time.  Is there some settings I need to change?


texts: Creating landslide assessment with dynamic rainfall data using Google Earth Engine?

texts:
I am new to GEE and I am trying to create a landslide assessment which dynamically changes depending on the rainfall data from our source and shows the map on our website.

Currently, my processes are uploading all susceptibility maps which already processed using ArcGIS and calculate it, I can get the exact same result in GEE. However, since changes depend on the rainfall data, I would NOT want to calculate the rainfall data using IDW in ArcGIS and upload again to GEE just to get a new landslide map.

Is there a way to upload the rainfall data from our website and calculate the IDW in GEE together with other susceptibility maps so I could get a new landslide map?


texts: How can I load pre-computed Google Earth Engine datasets into the stand-alone Google Earth client?

texts: The Google Earth Engine website lists a number of featured datasets.  When I click on one of the time lapse datasets, it opens in my browser (Firefox on Linux), but for the "Precomputed datasets", I get "The Google Earth Plugin is currently only available on Windows and Mac OS X 10.6+.".  That's okay, I have the standalone client.  But do I load the examples listed on the website — in particular, the two Global Roadless Area datasets — into the stand-alone Google Earth client?


texts: Can a spatial View in ArcSDE be viewed in ArcCatalog?

texts:
Can a spatial view in ArcSDE created by "sdetable create view" command be viewed as feature class in ArcCatalog, or is it just another plain table?

I have a feature class table and an attribute table (non spatial) stored in ArcSDE database. They can be related by a command column. I plan to create a spatial view in ArcSDE joining those tables and also publish the view as a feature service in ArcGIS Server. Is it possible? And, if the view can be published as a feature service, can feature operations such as Add, Edit and Delete Features be applied to all the view fields?


texts: ArcGIS Online Web Tile Layer for DroneDeploy in QGIS

texts:
I cannot seem to make this work and cannot find an answer on the site.
DroneDeploy has an app to "seamlessly generate template links for your maps to use as web tile layers in ArcGIS Online." see https://support.dronedeploy.com/docs/arcgis-online-web-tile-layer
Is there a way to load these tiles into QGIS?
I've tried loading via Layer > Add Layer > ADD ARCGIS Mapserver Layer
but I get the following Error:
Failed to retrieve service capabilities
network error: error downloading


texts: Can ArcMap open shapefiles made by QGIS?

texts:
I have created many EPSG shapefiles in QGIS and would like to import them into ArcMap for analyzing.

Can ArcMap open shapefiles made using QGIS?


texts: Can one convert text annotations to a layer in QGIS?

texts: I know text annotations can be turned on & off with Ctrl+T... but is there any way to export them as a layer?


texts:
I'm trying to calculate the average Julian day of year for greenness decrease from the MODIS MCD12Q2 dataset. Following is the code and error, which comes from having to normalize the julian day within each year using the apply_offset() function. There are several similar questions on here with this error, but I can't figure out how to fix it here.

// The values of MCD12Q2 are "Days since Jan. 1, 2000"
// This is the offset needed to convert to the relative Julian day for the year
var all_doys = ee.List.sequence(1,366);
var apply_offset = function(image){
  var day_offset = image.date().difference('2000-01-01','day');
  return image.subtract(day_offset).multiply(-1);
}

var all_years = ee.ImageCollection('MODIS/MCD12Q2')
            .select(['Onset_Greenness_Decrease1'])
            .map(apply_offset);

print(all_years)
var mean_decrease = all_years.reduce(ee.Reducer.mean());
print(mean_decrease)

var viz_options = {
  min: 0,
  max: 365,
  palette: ['0f17ff', 'b11406', 'f1ff23'],
};
Map.addLayer(
    mean_decrease.select('Onset_Greenness_Decrease1_mean'), viz_options,
    'Mean Decrease');

The error:

Mean Decrease: Tile error: Expected a homogeneous image collection,
but an image with an incompatible band was encountered.
Mismatched type for band 'Onset_Greenness_Decrease1':
Expected: Type<Float<-65169.0, 366.0>>.
  Actual: Type<Float<-64804.0, 731.0>>.
          Image ID: null
This band might require an explicit cast.


texts: How to fill gaps in a multipart polyline

texts:
I want a function that can fill small gaps between parts of a multipart polyline.

The polylines represent events on roads and are linear, not branched. The parts may be separated by small or large gaps. I want to connect any parts that are separated by a small gap less than some tolerance, while leaving the large gaps. No part is shorter than the gap tolerance.

The parts are (probably) well ordered and directed (meaning that geographically adjacent parts are adjacent in the IGeometryCollection, and parts all point in the same direction), but bonus points for an implementation that can handle misordered or flipped parts.

I'm working in ArcObjects, C#, .NET 4.0, ArcGIS 10.1SP1.

Update2:
Some example test data that, when filled with a gap tolerance of 0.25, should yield 3 parts and length 7.0910109853324:

//first two parts have no gap (thus they could really be a single part,
//but this might happen. They'll get merged in the output)
gcoll.AddGeometry(MakePath(MakePoint(0.0, 0.0), MakePoint(1.0, 0.0)));
gcoll.AddGeometry(MakePath(MakePoint(1.0, 0.0), MakePoint(2.0, 0.0)));
//third part has small gap and is flipped.
//it is also set at an angle to check that we don't remove a vertex.
gcoll.AddGeometry(MakePath(MakePoint(3.0, 0.0), MakePoint(2.1, 0.1)));
//gcoll.AddGeometry(MakePath(MakePoint(2.1, 0.1), MakePoint(3.0, 0.0)));
//fourth part has large gap
gcoll.AddGeometry(MakePath(MakePoint(5.0, 0.0), MakePoint(6.0, 0.0)));
//fifth part has small gap
//it is also set at an angle to check that we don't remove a vertex.
gcoll.AddGeometry(MakePath(MakePoint(6.1, 0.1), MakePoint(8.0, 0.0)));
//sixth part has large gap
gcoll.AddGeometry(MakePath(MakePoint(9.0, 0.0), MakePoint(10.0, 0.0)));


texts:
I'm doing research on chlorophyll using GCOM-C / SGLI L3 Chlorophyll-a concentration (V2) satellite data. And I want to make two UI designs.
First design I want to use Select to select the month that you want to see the information.
But I had a problem because I could not make the select functions according to the desired concept.
Here I attach the source code that I have created.
Could you help me please to solve my problem?
Your answer will help me in research.
var initIndex = 'Month';
var indexUrl = ui.url.get('index', initIndex);
ui.url.set('index', indexUrl);

var Month = ['Januari','Februari','Maret','April','Mei','Juni','Juli','Agustus','September','Oktobor','November','Desember']
var select = ui.Select({
  items: Month, value: ui.url.get('index'), style : {stretch: 'horizontal'}
});
select.onChange (optionChange)

// Set a place holder.
select.setPlaceholder('Pilih Bulan')

function setParams() {
  ui.url.set('Month', MonthSelect.getValue());
}


texts: Filter through geometries of a feature collection in Google Earth Engine to export as a shapefile

texts:
I'm trying to export a shapefile of administrative units out of GEE, but I get the error Shapefiles cannot contain multiple geometry types; found 'LineString', 'Polygon'
Here is my code and the feature collection I'm trying to export:
// Create variable with FAO Admin 2 boundaries, adm2
var adm2 = ee.FeatureCollection('FAO/GAUL/2015/level2');

// filter adm2 to only keep the municipality of interest, Ceiba
var ceiba = adm2.filter(ee.Filter.eq('ADM2_NAME', 'Ceiba Municipio'));

// add this layer to the map and center there
// Map.addLayer(ceiba, {}, 'Ceiba');

// // Export the ceiba shapefile to Drive
Export.table.toDrive({
collection: ceiba,
description:'faoadmin2_ceiba_only',
fileFormat: 'SHP'
 });

// "error: Shapefiles cannot contain multiple geometry types; found 'LineString', 'Polygon'."

The problem is that the Geometry of the feature is a GeometryCollection. Within that geometry collection there are 118 geometries, and all but one are Polygons. Geometry 0 is a LineString. In fact, it looks like a mistake to me (see visualization below), a random line off of one polygon, so I would like to delete the LineString geometry type or to filter it out.
// explore Ceiba in the console
print(ceiba);

// There are 118 geometries, all but one are polygons. Geometry 0, is a LineString
// Visualize this line in the map
var line_in_ceiba = ee.Geometry.LineString([[-65.5991688684573,18.21315366178916], [-65.5991688684573,18.2131492026881]]);
var point_in_ceiba = ee.Geometry.Point([-65.5991688684573,18.21315366178916]); // just for visibility
Map.addLayer(line_in_ceiba , {color: 'FF0000'}, 'line in Ceiba');
Map.addLayer(point_in_ceiba, {color: 'FF0000'}, 'start point in Ceiba');
Map.centerObject(line_in_ceiba, 20);

So far, I've tried setting the geometries as a property based on this answer, adding it as a new column using map, and trying to filter by that property
ceiba = ceiba.map(function(f) {
  return f.set('geo_type', f.geometry().geometries())
});

print(ceiba);
print(ceiba.filter(ee.Filter.eq('geo_type', 'Polygon')));

However, this doesn't quite work because the new column "geotype" now has a list of geometries, not individual elements.
How may I delete or filter out the LineString geometry type?
My main issue is that I need some of the FAO GAUL 2015 data as a shapefile, If there's an alternative way to export (or even to access it outside of GEE), I'm open to that too!


texts: Moisture Stress Index - how important is seasonal influence in a deciduous forest?

texts:
I want to calculate moisture stress index for a deciduous forest using the formula from the article 'Detecting vegetation leaf water content using reflectance
in the optical domain' (Ceccato et al. 2001) (SWIR/NIR). However, my forest browns/loses its leaves in autumn. This leads to the NIR decreasing, which leads to lower MSI values (meaning less water stress). However, this does not seem an accurate reflection of reality, as my area actually records far more rainfall in autumn than in other seasons. Also, the true colour images show that in autumn the grass fields are more yellow/bare in summer, and thus summer should be more dry and have higher MSI values.

How to deal with this issue? Is the MSI from the linked article simply not suited for deciduous forest, and are there any suggestions for moisture stress index that are accurate in all seasons?

I've added a picture showing the seasonal differences (values for nir and swir are the original sentinel values, so not divided by 10000), and a graph showing the monthly rainfall in mm (clearly summer is drier)


texts: Hot Spot Analysis using ArcGIS for Desktop?

texts:
I have a point shape containing locations of cars and a line shape containing roads.
Now I want to identify hotspots/clusters.

I want to find hot spots where for example the speed of the cars is very high/ very low. How can I do this in ArcGIS for Desktop? Should I use the Hot Spot Analysis tool here?


texts: Slope regions from elevation points using QGIS 2.0

texts:
I am using QGIS to create a layer of slope region polygons (e.g. less than 30% slope, 30-50%, more than 50% categories). My data comes from point elevations and I interpolated using SAGA b-spline approximation. I then used the Raster Terrain plugin to calculate slope. However, after using the raster calculator and polygonizing the categories of interest, I still have a problem in terms of how neat the borders are. I read about region growing, but am not sure if I have enough information to do so/how to do it.

This above image resembles what I'm aiming for. The bottom image is what I currently have.


texts: Clip line vector by shape fails - error 1005

texts:
I'm experiencing problems with the Clip tool in QGIS 2.18.1.

I'm trying to clip the intersection between some polygons and lines from another shapefile.

Every time, I get a "1005 error". I tried to read the logs, but I'm not familiar enough to solve this.

Both of my shapefiles are in the same projection, and have been cleaned (GRASS v.clean).

Here is the log of the process:

  Uncaught error while executing algorithm
              Traceback (most recent call last):
                File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\core\GeoAlgorithm.py", line 203, in execute
                  self.processAlgorithm(progress)
                File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\algs\qgis\Clip.py", line 71, in processAlgorithm
                  source_layer.crs())
                File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\core\outputs.py", line 326, in getVectorWriter
                  crs, options)
                File "C:/PROGRA~1/QGIS2~1.18/apps/qgis/./python/plugins\processing\tools\vector.py", line 610, in init
                  uri = GEOM_TYPE_MAP[geometryType] + "?uuid=" + unicode(uuid.uuid4())
              KeyError: 1005

EDIT:
It would seem that the problem comes from my line shapefile, as the process works with another file. How can I identify what's causing the problem with these infos?

FINAL EDIT:
As Steven Kay pointed out, the error code expressed that my geometry was MULTILINESTRINGZ, which was not supported by the Clip for some reason. By using the "Convert geometry type" tool from QGIS's Geoprocessings, I could convert the geometry to LINESTRING and finally do the Clip.


texts: Getting first and second row with ArcPy searchcursor?

texts:
I have a shapefile with only 4 rows in it. I already could filter it so I have only 2 left. I need these sorted by the 'length' field, then take the first and second row to save in to individual files. Here is the code for sorting:

import arcpy
fc = 'C:\\Users\\test.shp'
fields = ['length', 'FID']
expression = arcpy.AddFieldDelimiters(fc, "length") + " > 1000"
for row in sorted(arcpy.da.SearchCursor(fc, fields,)):
    print('{0}, {1}'.format(row[0], row[1]))

For saving in to the new file:

arcpy.FeatureClassToFeatureClass_conversion(?, outLocation,
                                            outFeatureClass,)

How can I make the search cursor to show first only the lower value row, then the higher one? (order doesn't matter)

How can I select the number 1 row, then the number 2 row?


texts: Calling qgis:postgisexecuteandloadsql from standalone script. DataItemProvider is empty

texts:
I'm trying to execute sql and load layer from script. whcih works fine with qgis python console but standalone script fails.
There may be a proper way to initialize postgresql db settings but so far I have no data how to solve.
import os, sys
# Append QGIS Python library to python search path
sys.path.append(r'C:\Program Files\QGIS 3.18\apps\qgis\python')
# PyQt5
sys.path.append(r'C:\Program Files\QGIS 3.18\apps\Python37\lib\site-packages')
#processing
sys.path.append(r'C:\Program Files\QGIS 3.18\apps\qgis\python\plugins')
#set qt5 path
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH']=r'c:\Program Files\QGIS 3.18\apps\Qt5\plugins'

import processing
from qgis.core import QgsApplication,QgsDataItemProvider
from processing.core.Processing import Processing

QgsApplication.setPrefixPath(r"C:\Program Files\QGIS 3.18\apps\qgis", True)

showGui=True
qgs=QgsApplication([], showGui)
QgsApplication.initQgis()
Processing().initialize()

from qgis.core import QgsProviderRegistry

#here standalone script returns empty array.
print(QgsProviderRegistry.instance().dataItemProviders('postgres') )

#exception thrown.
dbResult= processing.run("qgis:postgisexecuteandloadsql",{"DATABASE":"mydb","SQL":"select * from table1","ID_FIELD":"id1"})
QgsApplication.exitQgis()


texts:
I'm having a hard time installing grid_renderer on Ubuntu Maverick (10.10). I've followed these instructions but when I compile the mapnik branch I always get this error (http://pastebin.com/Txc5Wt9W):

[...]
scons: *** [bindings/python/mapnik_python.os] Error 1
scons: building terminated because of errors.


texts: Adding KML layers to ArcGIS JS API v4 (beta)?

texts:
I have been building a web application using ArcGIS JS API v4.0.Beta1. I am having difficulties finding out how to add KML layers to the scene. How do I add the KMLs? I have tried adding a new feature layer with the URL from the host server and I am having no luck.

Example of what I've tried that doesn't work:

var featureLayer = new FeatureLayer({
url: "[insert KML url here]",
mode: FeatureLayer.MODE_SNAPSHOT

Any ideas?

Additionally: Here is my JSFiddle if you wanted to take a look for yourself

Here is my current .js as of 8/5/2015:

var map, view;

require([
  "esri/Map",
  "esri/views/SceneView",
  "esri/layers/KMLLayer",
  "dojo/dom",
  "dojo/on",
  "dojo/domReady!"
  "dojo/parser",
  "dojo/dom-style",
],
function (Map, SceneView, KMLLayer, domStyle, parser) {

  //Create map
  map = new Map({
    basemap: "topo"
  });

  //Create SceneView
  view = new SceneView({
    map: map,
    container: "viewDiv",
    camera: {
      position: [-105.581, 41.305, 3000],
      tilt: 45
    }
  });

    parser.parse();

    var kmlURL = 'http://uwyo.maps.arcgis.com/sharing/rest/content/items/457a85b101ba4fb6a159d4ea16c75153/data';
    var kml = new esri.layers.KMLLayer(kmlURL);
    map.addLayer(kml);
      kml.on("load", function() {
      domStyle.set("loading", "display", "none");
      });

});


texts:
I have :

FeatureClass (Structure)
column (Type)
column (OID)

I am trying to get the (Type) for specific OID.
I tried to use this code but it dose not work.

IFeatureWorkspace featureWorkspace = (IFeatureWorkspace)workspace;
IFeatureClass testFeatureClass = featureWorkspace.OpenFeatureClass("Structure");
IQueryFilter queryFilter = new QueryFilterClass();
queryFilter.SubFields = "Type";
queryFilter.WhereClause = "OBJECTID = " + feature.OID;
IFeatureCursor fCursor = testFeatureClass.Search(queryFilter, true);

I don't know what is the wrong, any help would be appreciated.


texts:
I have a KML file that I'm trying to open on QGIS. However, when I add them in QGIS I obtain a raster and a polygon but the raster do not fit in the polygon.
I guess the polygon contain geovalue while the raster do not. How can I manage that?


texts: Searching free data source: Shapes of counties worldwide

texts:
I want to create a map with openlayers, that contains a vector layer, with all "counties" of all countries of the world. "County" doesn't mean any strict political term here. I just want to split the earth into areas of ~500 - ~2500 km2 (or even much larger in sparsely populated areas), while these areas should have names, that are common to the people living there.

Now I'm searching data sources.

The data should be free
Each county information should contain the name in one of the languages spoken in this county.
It is not important, that the borders are exact - it is just used for a visual selection of some counties, the user wants to do something with
It must be complete for each country, that is available (no areas without a "county").
For now, it would fit to have those data available for Germany, Austria and Switzerland, but it would be fine to find a solution, that would allow us, to add other countries without having another data import and maintainance workflow for each country.


texts:
I've been trying to follow the documentation for cartodb for a while to figure this out, but haven't found it very clear. All I'm trying to do is hide one layer and then show another layer upon button click.

However, I have no idea of the syntax and the documentation wasn't much help.

$(document).ready(function() {

var map2 = new L.Map('attribute_map', {
center: [39.8282,-98.5795],
zoom: 4,

})
var map2;
$('#load_2').click(function() {

$('#load_2').css("display","none");
$ ('#attribute_map').css("height","520px");

 var mapURL = 'https://cartodb-basemaps-{s}.global.ssl.fastly.net/light_nolabels/{z}/{x}/{y}.png'
L.tileLayer(mapURL, {

}).addTo(map2);
 cartodb.createLayer(map2, 'https://jonmrich.cartodb.com/api/v2/viz/50366df6-afed-11e4-ab33-0e853d047bba/viz.json', { https: true,legends: true,cartodb_logo:false,layerIndex:1 })
.addTo(map2)

});
$('#load_3').click(function() {

 cartodb.createLayer(map2, 'https://jonmrich.cartodb.com/api/v2/viz/0a63ff26-ad47-11e4-bebb-0e4fddd5de28/viz.json', { https: true,legends: false,cartodb_logo:false, })
.addTo(map2)
 });
  });

Can someone point me in remotely the right direction? In this example, I'm trying to add visualizations based on different tables. I also tried doing it for multiple layers in the same visualization, but couldn't figure that out either.


texts: A problem with event firing upon mouseout

texts:
I'm trying to create a point map that displays information about a point in an info box upon user mouseover and changes the styling of the point. On mouseout, I want the info box and marker styling to revert back to empty/default styling. On click, I want the map to zoom-in to the point.

As written, the code changes the style of a marker upon mouseover. Upon click, the map zooms-in on the point. However, the feature does not revert back to its original styling upon mouseout. Additionally, the info box does not display any information about the feature. The zoom when clicked does work.

I'm very new to leaflet and javascript.

Below is the code that create the info box in which I will display information about the feature moused over:

    // Map control to present informaiton
    var infobox = L.control ({
        position: 'bottomright'
    });
    infobox.onAdd = function (e) {
        this._div = L.DomUtil.create('div', 'info');
        this.refresh();
        return this._div;
    };
    infobox.refresh = function (properties) {
        this._div.innerHTML = '
Observation Information
';
        if (typeof (properties) != 'undefined') {
            this._div.innerHTML += 'Refernce #: ' + properties.RefNum + '
'
            + 'Group: ' + properties.Group + '<br/'
            + 'Plaintiff Type: ' + properties.PlaType + '
'
            + 'Click to zoom.';
                } else {
                    this._div.innerHTML += 'Hover an observation for information.';
                }
        };
        infobox.addTo(map);

Here is the code that reads in the geoJson file with the point features:

    //create layer for all observations, load geoJson object with all observations, define the marker type, bind pop-up using onEachFeature, add to map
    var allObs = L.geoJson(allobs, {
        pointToLayer : function (feature, latlng) {
            switch (feature.properties.Group) {
                case "Control": return L.circleMarker(latlng, controlMarkerOptions);
                case "Limited": return L.circleMarker(latlng, limitedMarkerOptions);
                case "Maximal": return L.circleMarker(latlng, maxMarkerOptions);
            }
        },
        onEachFeature : function (feature, layer) {
            layer.on ({
                'mousemove': function(e) {
                    e.target.setStyle({
                        weight:3
                    });
                    infobox.refresh(feature.proprtyies);
                },
                'mouseout': function (e) {
                    allObs.resetStyle(e.target);
                    infobox.refresh();
                },
                'click': function (e) {
                    map.fitBounds(e.target.getBounds());
                }
            });
        }
    }).addTo(map);


texts: Add indefinite number of answers to Survey123

texts:
I am looking to create a Survey123 project for a sampling procedure. The sampling will entail the "surveyor" standing at an apparatus capturing information as is needed.

The issue I am having is the number of records per point of capture will be different each time the "surveyor" visits a location.

The procedure will measure a distance and then a value at that given distance. The distance value will change based on site and the total number of distance intervals will be different.

Does Survey123 allow for this type of dynamic data entry?
I haven't found a an answer through online searches, but this may be due the incorrect use of the terminology.

I haven't been able to find anything that allows for this using the Web Creator,


texts:
I am trying to understand the output of the heatmap tool of QGIS (kernel density). For this reason I created a kernel density map for a single point that has a projected coordinate system (epsg: 21781). I used the following parameters:

radius: 200m
Cell size X and Y: 5m
kernel shape: "Quadratic (biweight)"
output value: "Raw values" and "Scaled by kernel size"
=> I tried both options

On first sight the outputs of the heatmap generated with the output value "Raw values" and "Scaled by kernal size" seem identical as there is no visual difference between them. However, at second sight I realized that the legends differ.

My questions are:

Am I right to assume that "Raw values" refers to the total quantity of sand in each cell and that "Scaled by kernal size" refers to the amount of sand per unit area? (See this link)
What is the unit of the output values "Raw values" and "Scaled by kernal size"  assuming that my points are burglary incidents in a city? Am I right to assume the following:

"Raw values" = number of burglary per cell ?

"Scaled by kernal size" = number of burglary per m2 (my layer unit is meter)

Is the term "kernal size" a synonym for the chosen radius (in my case 200m) or does it refer to the cell area (in my case 25m2)?
As described here I tried to multiply the density values by the cell area and added up these values over the grid in order to get the total sum of the original data. My cell area is 25m2 and I firstly multiplied this value with all cell values of the heatmap produced with the "Scaled by kernal size" output value. Secondly I summed up all multiplication results: value of cell 1 *25m2 + value of cell 2 * 25m2, + ... + value of cell n * 25m2. I expected to get a total sum of one as I created a kernel density map for a single burglary incident (In case my points represented 50 burglary incidents, I expected the total sum to be 50). However, I got a much higher value (41887). Maybe I misunderstood the meaning of "total sum of the original data"?
I expected the sum of all cell values of the heatmap produced with the "Raw values" also to equal one (= total amount of sand = in my case one burglary?). However, this is  not the result I got (1675). Why?


texts: GEE NDVI time series extraction does not populate CSV file columns

texts:
I am trying to extract a time series of sentinel NDVI data. This script (below) takes a feature collection of points, buffers around them, and collects all available imagery 30 days before each feature date.
The output should be a CSV file with two additional columns: "NDVIs" and "dates." Each row in these columns will be an array of NDVI and unix dates that correspond to each feature date.
However, the task does not fill the columns. Instead the resulting CSV has '[]' in both columns. See here:

Any idea what could be going wrong here?
########################################################################################
#################            NDVI at locations                ##########################
########################################################################################

import ee
#ee.Authenticate()
ee.Initialize()
l8 = ee.ImageCollection("COPERNICUS/S2")

table = ee.FeatureCollection("users/ddlawton/LocustDat_500")

def locational(feature):

  end = ee.Number(feature.get('unixdates'))
  start = end.subtract(2592000000) #subtract 30 days (in milliseconds)
##initialize lat and long
  lat = ee.Number(feature.get('latitude'))
  longi = ee.Number(feature.get('longitude'))

##initialize an EE geometry point with lat and long and buffer
  geom = ee.Geometry.Point([longi,lat])
  buffed = geom.buffer(500)

  #start = ee.Date(1496793600000)
  #end = ee.Date(1528329600000)

  #end = ee.Date('2019-06-01')
  #start = ee.Date('2019-05-01')

  ########define date range of interest

  fCollectp = l8.filterBounds(buffed)
  fCollect = fCollectp.filter(ee.Filter.date(start, end))
  #filterdate range

  def addDates(image):
    imagep = ee.Image(image)
    setted = None
    nir = imagep.select('B5');
    red = imagep.select('B4');
    ndviImage = nir.subtract(red).divide(nir.add(red)).rename('NDVI');
    NDVIavgp = ndviImage.reduceRegion(ee.Reducer.mean(),buffed,scale=30,maxPixels=1e9)
    NDVIavg = NDVIavgp.get('NDVI')
    NDVIdateCap = ee.String(image.get('GENERATION_TIME'))

    tempA = ee.List([NDVIavg,-9999])
    tempB = tempA.reduce(ee.Reducer.firstNonNull())#this will filter out a null if the image was not actually over the gauge

    NDVISetted = ee.Image(image.set('date',NDVIdateCap).set('NDVI',tempB))

    return NDVISetted

  mapC = fCollect.map(addDates)
  agg = mapC.aggregate_array('date')
  aggb = mapC.aggregate_array('NDVI')
  agFeat = feature.set('dates',agg).set('NDVIs', aggb)
  return agFeat

mapF = table.map(locational)

task = ee.batch.Export.table.toDrive(
  collection=mapF,
  description='Locust_outbreaks_jan_23',
  fileFormat='CSV'
)

task.start()


texts: GDAL GetGeoTransform Documentation -- Is there an oversight, or what am I misunderstanding?

texts:
The documentation for GDALDataset::GetGeoTransform says some things that are rather confusing.

Fetches the coefficients for transforming between pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space.
Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];

In a north up image, padfTransform[1] is the pixel width, and padfTransform[5] is the pixel height.
The upper left corner of the upper left pixel is at position (padfTransform[0],padfTransform[3]).
The default transform is (0,1,0,0,0,1) and should be returned even when a CE_Failure error is returned, such as for formats that don't support transformation to projection coordinates.

If this is a matrix, it makes sense that Pixel_X is at 1,1 and Pixel_Y is at 2,2.
Here is my question:
So should the equation for Yp be:
Yp = padfTransform[3] + L*padfTransform[4] + P*padfTransform[5];

in which the positions of L and P are swapped?
Or if I am confused about something, can someone help me get unconfused?


texts: Customize header for coordinates in CSV export

texts:
Is it possible to customize the header for csv export from QGIS?
I use "Export -> Save Features as" to create a csv with the names and coordinates for a list of points.
The default when choosing Geometry AS_XYZ is "X,Y,Z,name" - would it be possible to get "easting,northing,elevation,name" instead?


texts: Choosing between QGIS and GRASS for simple thematic maps with different layers

texts:
I was working with GIS and remote sensing in a university and professional environment some ten years ago (ArcInfo etc), so I have some experience even though somewhat outdated.

Now I am looking for an GIS application for a small farming project and I found GRASS and QGIS projects which both sound very interesting.
However, I do not really get the difference between the two programs.

Which one do I use best for simple thematic maps with different layers?


texts: Openlayers 3.18.2 vector tiles - horizontal and vertical line artifacts

texts:
When rendering a vector tile layer in Openlayers 3.18.2 i get horizontal and vertical line artifacts on some zoom levels. When zooming in or out they vanish or appear on other places. I tried to use different render modes, with no effects. It is also browser independent (tested on Chrome and Firefox)

Layer Code:

var vector_tile_layer = new ol.layer.VectorTile({
        source: new ol.source.VectorTile({
          attributions: attribution,
          format: new ol.format.MVT(),
          tileGrid: ol.tilegrid.createXYZ({maxZoom: 20}),
          tilePixelRatio: 16,
          overlaps: false,
          url: 'http://localhost:8080/data/{z}/{x}/{y}.pbf'
        }),
        style: function(feature) {
          var kgnr = feature.get('KG_NR');
          if (kgnr % 3 == 0) {
            color = '#000000';
          }else if(kgnr % 3 == 1){
            color = '#91AA00';
          }else if (kgnr % 3 == 2){
            color = '#F1A4C3';
          }

            style = new ol.style.Style({
              stroke: new ol.style.Stroke({
                color: color,
                width: 1
              })

            });

          return style;
        },
        options: {
            'renderMode': 'vector'
        },
      });

Sample image (all straight vertical and straight horizontal lines are artifacts):

Any ideas anyone?


texts:
When trying to 'Geocode' a CSV file I receive this answer:

UnboundLocalError: local variable 'attributes' referenced before assignment
Traceback (most recent call last):
File "/home/user/.local/share/QGIS/QGIS3/profiles/default/python/plugins/mmqgis/mmqgis_dialogs.py", line 1384, in run output_file_name, not_found_file_name, self.mmqgis_status_callback)
File "/home/user/.local/share/QGIS/QGIS3/profiles/default/python/plugins/mmqgis/mmqgis_library.py", line 2087, in mmqgis_geocode_web_service notwriter.writerow(attributes) UnboundLocalError: local variable 'attributes' referenced before assignment
Versión Python: 3.6.8 (default, Jan 14 2019, 11:02:34) [GCC 8.0.1
20180414 (experimental) [trunk revision 259383]]
Versión QGIS: 3.8.1-Zanzibar Zanzibar, dcd95cc648

I've already read a lot of posts and this, but I don't understand it.


texts: I have been working on a QGIS plugin. I have run into some problems and to solve it, one of the solutions would be to add a column to the WFS layer. However most of the times, I have not the permission to edit the WFS layer. I was wondering whether it would be possible to create a kind of temp attribute to each feature. It does not have to be persistent, so I only needed at runtime. Would something like this be possible?


texts: Is there any efficient way to convert GeoJSON to WKT?

texts:
I want to use ST_GeomFromText() to load vector data to PostGIS. However, I don't know how to convert from Geojson type to WKT.
N.B. I know I can use OGR2OGR to upload the GeoJSON directly, but I just want to confirm if there is any way to convert it or not.


texts:
I have some code that will take a ArcGIS 10.2.2 feature class and post a copy of it in a different geodatabase.

My code is set up so the end user can pick the feature class they want the model to run (getparameterastext).  My problem is the code doesn't name the new feature class the same name as the original feature class.  My code "hard codes" a name to the new feature class being posted.  What do I need to change the my code so the new file will name the new feature class the same name of the original feature class the end user pics?  This is what my code looks like.

# Import arcpy module
import arcpy

# Local variables:
maint_FC = arcpy.GetParameterAsText(0)
sdr_staging_gdb = "Y:\\data\\dev\\staging\\sdr_staging.gdb"
staging = "Y:\\data\\dev\\staging"

# Process: Feature Class to Feature Class
arcpy.FeatureClassToFeatureClass_conversion(maint_FC, sdr_staging_gdb, 'maint_FC', "", "POL_DIST \"POL_DIST\" true true false 3 Text 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,POL_DIST,-1,-1;INSPECTOR \"INSPECTOR\" true true false 50 Text 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,INSPECTOR,-1,-1;INSP_NBR \"INSP_NBR\" true true false 3 Text 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,INSP_NBR,-1,-1;COMMENT \"COMMENT\" true true false 50 Text 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,COMMENT,-1,-1;DATE_MOD \"DATE_MOD\" true true false 36 Date 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,DATE_MOD,-1,-1;DATE_ADDED \"DATE_ADDED\" true true false 36 Date 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,DATE_ADDED,-1,-1;MOD_BY \"MOD_BY\" true true false 50 Text 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,MOD_BY,-1,-1;ADDED_BY \"ADDED_BY\" true true false 50 Text 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,ADDED_BY,-1,-1;Shape_STArea__ \"Shape_STArea__\" false false true 0 Double 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,Shape.STArea(),-1,-1;Shape_STLength__ \"Shape_STLength__\" false false true 0 Double 0 0 ,First,#,Database Connections\\UGGSSSQL_GISPUB_basemap.sde\\gispub.BASEMAP.gisdata\\gispub.BASEMAP.TL_TEST_2,Shape.STLength(),-1,-1", "")

# Process: Feature Class To Shapefile (multiple)
arcpy.FeatureClassToShapefile_conversion(maint_FC, staging)


texts:
I'd like to import Coordinates from an Excel list, for this reason I saved the excel sheet as an .csv ant tried to import it via "Layer --> Add delimited text Layer". The points are imported, but with the wrong Coordinates. They are very far away from the right location. They should be in Germany, but they are in the Sea near Africa.

In the table of attributes and in the import dialogue the coordinates are displayed correctly.

I tried:

to exchange Longitude and Latitude.
Save the csv in utf8 (and choose utf8 in the QGIS dialogue)
format the Coordinates in Excel as a Number

Does anybody know what can be wrong?

(in the other question it is only asked how it works not, how to solve my problem)


texts: Convert KML files to vector tiles or images and show them on open street map

texts:
I am using Leaflet KML plugin as below:
https://github.com/windycom/leaflet-kml
Here is full code:

<meta name="viewport" content="width=device-width, initial-scale=1.0,  maximum-scale=1.0, user- scalable=no" />
<link rel="stylesheet" href="http://unpkg.com/leaflet@1.4.0/dist/leaflet.css" />


texts: Subtracting data between features using ArcMap

texts:
I have a layer that contains feature classes that store optimal elevation & existing elevation values for a place. Right now the existing and optimal are stored in different features under the same layer, yellow is existing and red is optimal.
Is there a good tool or method to group them into pairs because I need to subtract the existing from the optimal to create a feature or a field with the result?
I was going to modify the Python script of the Intersect tool and use a higher cluster tolerance to group the features and subtract them.
Attribute Table
ID  Elev_Value  Elev_type
1   69.5         Optimal
2   69.3         Optimal
3   70.1         Existing
4   70.7         Existing
5   70.2         Optimal
6   74.3         Existing
7   73           Optimal
8   77           Existing
9   75           Optimal
10  69.5         Optimal


texts:
Is there a way to change the icon of a geodatabase?

I have a geodatabase from ArcMap 10.2 (have now upgraded to 10.7) and want to distinguish between the old and the newer geodatabases.


texts: I have calculated NDVI and SAVI for landsat-8 image. But I am getting higher value of NDVI than SAVI Can you tell me reason for this NDVI=0.46 SAVI=0.34 Specific band pixel values are toa b4- 0.130 toa b5- 0.365 L= 0.5 taken


texts:
I am developing a web application that brings in a number of tileLayers and geoJSON layers. For the tile layers I have been able to do similar to this JSFiddle in that the layers appear depening on their zIndex value regardless of which layer was called in last.

I thought I would be able to do the same for geoJSONS, however when i add/ remove the layers, the last layer that is added in appears on the top, despite which order they have been written in.

After checking through the documentation I am to understand that the zIndex parameter only applies to tileLayers and at least in version 0.7.5 cannot be used with geoJSON.

I have read that a potential workaround is to use the bringToFront/Back() method, however this does not appear to work and I can't imagine will provide the best solution for reading in a large amount of geoJSONS.

If there is no official method with 0.7.5 is there at least a workaround? I have been bringing in the geoJSONS as below.

var testLayer = new L.geoJson(json_testdata, {
style: testStyle
});
testLayer.addData(json_testdata);

Update: This is a JSFiddle that shows the tileLayers working with a ZIndex, however when the GeoJSON layers are added and removed the last one to be added goes on top. Is it possible to maintain the layers position? i.e. the order it is written or assigning an equivalent to a ZIndex?


texts: How to calculate quadruplet reference (east, west, north, south) of a parcels

texts:
I faced a challenge to calculate which parcels are located in east, west, north, south of a definite parcel?

I mean that when I select a feature, how to calculate that which feature is located in east or west ... of this feature?

In this picture, I want to calculate which feature is located to east and which feature is located in west?


texts: Can I create a reusable text annotation style in QGIS?

texts:
I want to create a fairly simple map with a few locations emphasized using annotations.

I don't like the default annotation formatting:

So, I've reformatted the annotation background color, border color, fonts, font sizes, marker size, etc. to something that I prefer:

How do I re-use this same formatting on several other annotations without having to re-set each of the text annotations' formatting by hand?  I cannot copy and paste one annotation (that would be nice).  I believe the html annotation might be a solution, but I cannot find any good documentation describing how to use it.

An even better solution would be to automatically label each feature using this format, and pulling the content of the label from fields within the feature layer.  Hmm, perhaps I should investigate labels instead of annotation?


texts: Auto wrap legend text based on character count

texts:
In QGIS 2.18 there is an option to wrap text based on a certain character. What I need is to automatically wrap my text on the ' ' (space) character IF the character length exceeds the size of my legend box. Question is: how do I get an expression builder in the supplied warp box?

In this example I warped on the caracter '_'


texts:
I am trying to calculate the value of a number of attribute tables as shown below (Using the graphic modeler). In the first image, I outline the process whereby I have three polygons and a point feature, each with a row of values derived from calculations such as number of points inside polygon and does the point intercept polygon XYZ.

Essentially, I need to add predefined values seen in the table in image 2 to the number in each cell then add these together to get a total for each shapefile and point feature. I then need to calculate the SUM value of all the shapefiles and point features as a cumulative value.

The table of values which are used for the calculation may need to be edited before running the script but will otherwise remain as default.

I have attempted to simplify this in the images below.

Image 1

Image 2


texts: Remove/Delete a feature from OpenLayers vector layer

texts:
I have an OpenLayers layer with several vectors on it. I know the OpenLayers Feature ID (e.g. OpenLayers.Feature.Vector_241) and I would like to know how to remove this feature and only this feature?

I have considered deleting them all and then looping through an array to redraw them, but this seems a little inelegant, and long winded.

I have tried a variety of versions of the folllowing code:

layer.removeFeatures( featureID );

And

layer.removeFeatures( [featureID] );

And

layer.removeFeatures( {feature : featureID });

Etc...

All help is welcome!

Thanks

C


texts: Geodjango distance between two points appears incorrect

texts:
I want to show you something. I am using django+postgis and I am experiencing wrong distance values between two points.

This is my model:

class Location(models.Model):
    name = models.CharField(max_length=255)
    latitude = models.FloatField(blank=True, null=True, verbose_name='Latitude')
    longitude = models.FloatField(blank=True, null=True, verbose_name='Longitude')
    location = models.PointField(blank = True, null=True)

    objects = models.GeoManager()

    def save(self, *args, **kwargs):
        self.location = Point(self.latitude, self.longitude)
        super(Location, self).save(*args, **kwargs)

This is the data:

bogota = Location(name='bogota', latitude='4.61779093242', longitude='-74.0893133545')
bogota.save()
barranquilla = Location(name='barranquilla', latitude='10.9731710819', longitude='-74.8020513916')
barranquilla.save()

This is the query:

ls = Location.objects.distance(bog.location).order_by('distance')
for l in ls: l.name; l.distance.km

This is the output:

Out[29]: u'BOGOTA'
Out[29]: 0.0
Out[29]: u'BARRANQUILLA'
Out[29]: 205.27667079004999

But Barranquilla city is to ~700km from Bogota city, see here.

So, what is wrong? the SRID are the default values.


texts: Mosaic to new raster on rasters with the same name in different locations

texts:
I have national woody vegetation raster files for Australia for years between 1988-2017. These files are separated into spatial tiles with the same name. So, what I would like to do is combine raster files with the same name together. Here is my code:

import os, arcpy

workspace  = r"C:\Users\lawtond\Rasters\Woody"
list_raster= [] # the list must exist before you can append
walk       = arcpy.da.Walk(workspace, type="GRID")
output_location  = r"C:\Users\lawtond\Rasters\Woody"

for dirpath, dirnames, filenames in walk:
    for file in filenames:
        if "woody13" in file.lower():
            list_raster.append(os.path.join(dirpath,file)) # FULL path to each raster

print (list_raster)

arcpy.MosaicToNewRaster_management (list_raster, output_location, "merged_woody_2013.tif","1", "FIRST")

However I get the following error:

Message File Name   Line    Position
Traceback
    <module>    <module3>   18
    MosaicToNewRaster   C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py 13191
ExecuteError: ERROR 000622: Failed to execute (Mosaic To New Raster). Parameters are not valid.
ERROR 000628: Cannot set input into parameter coordinate_system_for_the_raster.

What am I doing wrong here?


texts:
I want to process a stack of Sentinel 1 SLC images (4) to produce tomographic layers for above ground biomass estimation. To do that, I need to generate the interferometric vertical wave number images for further processing.

Is there a software which can help me in this case>?