texts:
I found an upload example to load WFS data from geoserver into a GeoJson layer in leaflet. I I adapted it to my case. I have read similar topics from this site  and this other site, but I can not understand why this code does not work for me.

<!DOCTYPE html>

   <meta charset="UTF-8"/>
   <meta name="viewport" content="width=device-width, initial-scale=1"/>
   <title>Geoserver</title>

   <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.6.4/leaflet.css"/>

Datasource of the WFS-Layer: http://opengeo.org

I am including JSONP in GeoServer

<context-param>
    <param-name>ENABLE_JSONP</param-name>
    <param-value>true</param-value>
</context-param>


texts: Getting distance between point and line with specific values?

texts:
I have a table named "housenumbers" that contains points with addresses fields, and one of them is "street_name". And I also have another table named "streets", that contains lines that are the geometries of roads, and a field that contains its "name".

I need to get the distance from each point to the nearest street that has the same street name.

Up to now, I have this:

select st_distance(housenumbers.geom,streets.geom),*
from housenumbers, streets
where housenumbers.street_name = streets.name

It's working, but it's showing all the distances from each point to all streets with the same name in the area. What I need is only the nearest one.


texts: Using attributes of polygon in one layer to identify polygons in another?

texts:
I am in the process of developing a web application, I am using open source tools so qgis, postgis, play framework for the web app. And I have come across an interesting problem that I'm trying to work on if someone might be able to point me in the right direction.

What I need to do is this I have one polygon layer A with attributes about the polygons, and polygon layer B with attributes also.

So if a polygon in layer B falls within a polygon in layer A look up the attributes for the polygon in A, now compare those with the attributes of the polygon in layer B to determine if the attribute for the polygon in B is a multiple ( 2x, 3x, 4x, etc ) of the attribute of the polygon in A.

I am going to be using the results of this in the web application for other things.

So I have multiple options, but looking for some direction as to where to start.

Do I need to use something like sextante?

There are several of these operations that I need to do. It Should be repeatable, but only need to be repeated when updates to the layers come out 1 or twice a year! But it will be acting on several states data layers.


texts:
I have a grid of cells with certain values (e.g. let's say from 1 to 6). I tried to create polygons associated to values above 4 with the help of arcpy.NumPyArrayToRaster then arcpy.RasterToPolygon_conversion.

To do this, I used python to store the values and transformed every value > 4 to 255 (0 otherwise) so that the raster create becomes white and black, and the polygons created are based only on the 255 values.

My problem is that the polygons created are "shifted" from the original grid cells:

I have no idea why I get this. The description of the grid, that I used to reconstruct the NumpyArray gave me correct values (e.g. lowerLeft point).

Any suggestion?

EDIT:
Here is my code:

import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")
env.extent = "MAXOF"
env.workspace = "C:/Project/ArcGIS/(grids)"
Coordsystem = "GEOGCS['GCS_ETRS_1989',DATUM['D_ETRS_1989',SPHEROID['GRS_1980',8183294.64,630550.669301,674150.669301,8145094.64]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
arcpy.env.outCoordinateSystem = Coordsystem
arcpy.env.overwriteOutput = True

# Extraction of data from the grid
gridCells = arcpy.da.SearchCursor("C:/Project/ArcGIS/(grids)/gridN1.dbf",["Data"])
gridCells_list = []
for row in gridCells:
    gridCells_list.append(row[0])

# Modification of data in a new list
gridCells_newList = []
for i in range(len(gridCells_list)):
    if gridCells_list[i] > 4:
        gridCells_newList.append(255)
    else:
        gridCells_newList.append(0)

# Grid description
pDesc = arcpy.Describe(gridN1)
pWidth = pDesc.extent.width
pHeight = pDesc.extent.height
lowerLeft = arcpy.Point(pDesc.extent.XMin,pDesc.extent.YMin)
gridCol = int(pWidth / 500)
gridRow = int(pHeight / 500)

# Array creation
gridArray = np.asarray(gridCells_newList).reshape(gridRow,gridCol)
gridArray = gridArray[::-1]

# Raster creation
gridRaster = arcpy.NumPyArrayToRaster(gridArray, lower_left_corner=lowerLeft, x_cell_size=500,y_cell_size=500,value_to_nodata=0)

# Execute RasterToPolygon
arcpy.RasterToPolygon_conversion(gridRaster, gridPolygons, "NO_SIMPLIFY", "VALUE")


texts: Unable to create geometric network. Feature classes unavailable

texts:
When I try to create a geometric network 4 of my 9 feature classes are 'unavailable'. When I click on the 'unavailable' button it states that each feature class has an 'Invalid 'Enabled' field type' (see image below).

This problem has me stumped as I can't find much information on this particular error. Does anyone know what this error is referring to or how to rectify it?


texts: Searching Open Source tool to load raster data into Oracle spatial

texts:
I'm looking for a tool to load raster data into oracle spatial.
I found several commercial products but no open source tool.

Is there anything available in the open source world?


texts: Label polygon table with points data that are inside polygon

texts:
I have polygon data and point data in PostGIS database. I want to add a new column to polygon table that contains some value from point data (one column) if point is inside polygon. If there are multiple points inside polygon I want to add field data from all points in one colum.
Is there any way to do this using postgis?


texts: Does ArcMap have non-breaking space for labels?

texts:
Is it possible to have a non-breaking space in ArcMap? I particularly need it for labels.

&amp; gives me a & but
&nbsp; doesn't work.

Is there any way to make a non-breaking space?

I'm using Arc10, but this question applies to previous versions as well.


texts:
I am trying to clip about 36 rasters by a polygon. However, I am running into an error. It isn't a coding problem with the loop because I receive the same error when I run the function on a single raster. Here is my code:

# Import system modules
import arcpy
import os

# Set environment settings
arcpy.env.workspace = r"H:\GIS Project\reprojected_rasters" #The file where your files are
arcpy.env.overwriteOutput = True
outWorkspace = r"H:\GIS Project\clipped_rasters" #where you want the new rasters to go
arcpy.env.parallelProcessingFactor = "50"
clipping_shp = r"H:\GIS Project\large AF dataset\point_feature_class\clipping.shp"
rasters = arcpy.ListRasters("*","TIF")
print rasters

arcpy.Clip_management ("SND_100_200_EV_N_P_AU_NAT_C_20140801.tif", "131.183300 -37.714100 152.250000 -16.056200", r"H:\GIS Project\clipped_rasters" , clipping_shp, "ClippingGeometry" , "NO_MAINTAIN_EXTENT")

print "Clipped Script ran ok"

However, when I try to run this code, I receive the following error:

Traceback (most recent call last):
  File "<module2>", line 21, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 13594, in Clip
    raise e
ExecuteError: ERROR 000622: Failed to execute (Clip). Parameters are not valid.
ERROR 000800: The value is not a member of ClippingGeometry | NONE.

Here is some information about my data:

all are .tif files
all are projected into the same datum (GDA 1994 Australian Albers)
The clipping polygon was made from using the minimum extent function of a point feature class


texts: Two Separate QGIS Plugins Under one Provider

texts:
I have developed a QGIS processing plugin that has a processing provider of its own. That plugin is working fine. Now I am developing a second plugin and I want to provide it the same provider so that I can have something like this where many algorithms are under Vector Selection Provider:

These two plugins are different and hence I don't want to package them under same plugin but do wants them to have same provider.
So far I have tried giving the second Plugin Provider the same name as first one but I am getting this error:
Error while unloading plugin curve_number_generator
RuntimeError: wrapped C/C++ object of type CurveNumberGeneratorProvider has been deleted
Traceback (most recent call last):
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 434, in unloadPlugin
    plugins[packageName].unload()
  File "C:/Users/abdul/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\curve_number_generator\curve_number_generator.py", line 74, in unload
    QgsApplication.processingRegistry().removeProvider(self.provider)
RuntimeError: wrapped C/C++ object of type CurveNumberGeneratorProvider has been deleted


texts: Show layer name in attribute table in print composer?

texts:
When adding attribute tables in the print composers it is not possible to distinguish between various tables from various layers; each table is generic and does not indicate which layer the table is from.
Is there a simple way to add an additional header or even a caption to indicate the source layer?


texts:
I have an Access database and have some built in queries in the database. The same data is stored in geodatabase with the spatial information. Now I want to follow this process:

Access-> push out query-> open ArcMap mxd -> display results.

Is there a way to do this?

I am using MS access 2010 and ArcGIS 10.1. Any help on this is most welcome.


texts: QGIS Azimuth and Distance plugin length issue

texts:
Relative newbie to QGIS and trying to figure out exactly what I'm doing wrong. A couple days ago I had azimuth and distance working perfectly. I was using it in a different project with different shapefiles and I'm not sure what projection I was using then.

I'm using QGIS Version 2.8.1-Wien and Azimuth and Distance Plugin Version 0.9.3

Now for some reason no matter what distance I put in it draws the line way off of the map. Before I was using default, but I've tried using both default and feet and it seems to make no difference.

I started over from scratch and everything is projected in EPSG:4326 WGS 84 and the project default is set to EPSG:4326 WGS 84. Whenever I draw to a temp layer it auto sets the projection to EPSG:4326 WGS 84 since one is not defined.

Here is a screenshot of what I'm inputting. Doesn't matter if I choose default or feet under units, it will still draw it off into infinity.


texts: Remove a raster from the composer legend, but display it on the map

texts:
I'm working with PyQGIS in QGIS 3.10.
I'm trying to remove a layer from the composer legend if it is a raster (but I want to see the raster on the map).
This is my code:
for i in self.newLayout.items():
    if isinstance(i,QgsLayoutItemLegend):
        legend = i
        legend.setAutoUpdateModel(True)
        for j in range(legend.model().rowCount()):
            modelindex=legend.model().index(j, 0)
            layertreenode=legend.model().index2node(modelindex)
            if layertreenode != QgsLayerTreeGroup :
                if layertreenode.layer().type()==QgsMapLayer.RasterLayer:
                    layertreenode.setItemVisibilityChecked (False)
        legend.setAutoUpdateModel(False)
        legend.setLegendFilterByMapEnabled(True)

It works... but a bit too much: my raster is not anymore visible in the legend, but the layer is set as invisible too (in the composer map and in the canvas)!
I thought I worked only in the legend model! May be I should have used another method like "removeFromLegend"? But I didn't find this anywhere.
More broadly, I think I have difficulties to understand boundaries between layer trees in the table of content and in then composer legend.
EDIT :
I read this thread :
Remove legend nodes of a categorized layer in a layout legend using PyQgis
an tried
if layertreenode.layer().type()==QgsMapLayer.RasterLayer:
    legend.model().removeRows(j,1)

It removes the legend item, but il also removes completely the raster from the QGIS project.


texts: Concatenating Attribut name for dynamic request

texts:
For my Workaround i'm using QGIS 3.12.0 and the expression editor.

The starting point is the following query:

round( (( "b_12052020" - "b_05052020" )/ "EWZ" )*100000)

The result is fine.

The feature database has one column for every day like b_01052020 and b_02052020 an so on.

The main goal is an automatic calculation from today to 7 days back.

I used the following code in an single expression and the result is fine.

format_date( (now() - to_interval( '7 days')),'ddMMyyyy')

I create with concat a string looks like the names of the columns.

concat( 'b_', format_date( (now() - to_interval( '7 days')),'ddMMyyyy'))

My Problem, i can't bring them together. for an example i show you a short version.

round( (( concat( 'b_', format_date( now(), 'ddMMyyyy')) - "b_05052020" )/ "EWZ" )*100000)

it brings me an error!

Also adding Quotation Marks to string brings no result.

round( (( concat( '"b_', format_date( now(), 'ddMMyyyy'),'"') - "b_05052020" )/ "EWZ" )*100000)

My question is, is it even possible to create dynamic table names for a query? is it only possible with the expression editor or do i need the funktion editor - and mayby how it works?
Thanks to all!


texts: Input field not showing up in Distance Matrix In QGIS 2.0?

texts:
I tried to perform a Distance Matrix analysis of building centroids. I have a two building centroid layers. The attributes in both the tables are name, osm_id, perimeter and area.

In the Distance Matrix Tool, I select this layer as input field then I have to select the unique field ID. My unique field ID is osm_id, but it does not show up in the list.

Does someone know why that could be?

I get all the other attributes.


texts:
So i have been trying to connect the following WFS resource to QGIS 3.4.4.

https://maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_Cadastre_WFS/MapServer

It connects without problems, lets me add the layer to the map without complaint but the layer is actually empty.
When I try to open the attribute table it says "Loading features for layer XY 0%" and never progresses.
This relates to the standard AddLayer/WFS tool. I have also tried to add it via the AddLayer/ArcGIS Feature Server tool. since the resource seems to lie on a ArcGIS Server, but it replied with en error message "Failed to retrieve service capabilities: Network error: Error transferring http://maps.six.nsw.gov.au/arcgis/services/public/NSW_Cadastre_WFS/MapServer/WFSServer?f=json - server replied: Bad Request"

The resource loads without problems in ArcGIS


texts: Remapping classified image using Descartes Labs

texts:
I am hoping to remap the values of a multi-value classified image to create a binary classification. In Google Earth Engine, there is a function remap() that takes the original values of pixels and remaps them to a new set, ex: [0, 1, 2, 3, 4, 5] to [0, 0, 0, 0, 1, 1].
Is there a comparable function in the Descartes Lab platform?
First, I identify the regions of California I am interested in and merge them:
import descarteslabs as dl
import numpy as np
from pprint import pprint
import json
import geojson
from functools import partial
#import pyproj
import shapely.geometry
import shapely.ops

places_client = dl.Places()
california = places_client.find(
    'california_redding'
    #'north-america_united-states_california'
)

places_client = dl.Places()
sac_geojson = places_client.shape(
    'north-america_united-states_california_sacramento-valley'
    #'north-america_united-states_california_san-joaquin-valley'
)

sj_geojson = places_client.shape(
    'north-america_united-states_california_san-joaquin-valley'
)

#pprint(new_mexico_geojson)

#sac_geo = sac_geojson.geometry
#sj_geo = sj_geojson.geometry

#sac_geo.geom_type

#geometry = sac_geo.union(sj_geo)

# pulling out the polygons
poly_sac = shapely.geometry.asShape(sac_geojson['geometry'])
poly_sj = shapely.geometry.asShape(sj_geojson['geometry'])

central_valley = poly_sac.union(poly_sj)

# using geojson module to convert from WKT back into GeoJSON format
geojson_central_valley = geojson.Feature(geometry=central_valley, properties={})

Next, I select the 2019 Cropland Data Layer from the catalog:
cropland_2019, ctx_19 = dl.scenes.search(geojson_out,
                    products='usda:cdl:v1',
                    start_datetime="2019-01-01",
                    end_datetime="2020-01-31",
                    limit=5
                   )

cld_19 = cropland_2019[0].ndarray('class', ctx_19)

cropland_2018, ctx_18 = dl.scenes.search(geojson_out,
                    products='usda:cdl:v1',
                    start_datetime="2018-01-01",
                    end_datetime="2019-01-31",
                    limit=5
                   )

cld_18 = cropland_2018[0].ndarray('class', ctx_18)

I am interested in remapping both the 2019 and 2018 cropland data layers for the central valley of California into a binary classification rather than a ~100 class classification, and then add together overlapping pixels of the two images together to create a single image.


texts: Raising raster to a power gives no data values

texts:
When performing the following operation in QGIS Raster calculator:

10 ^ (("20091201_20100101_chl_modis_malvinas_ns@1" * 0.0000581)-2)

I get no dada values for every cell where I should get < 1, and correct values for the rest.

I understand this has to be a silly thing that I am overlooking, but I could not make it work, nor find answers out there.

The same operation and dataset used to work OK with the Raster Calc plugin, but it appears it is not supported in QGIS 2.4.

Any other operation with the same layer works ok, rendering real numbers.


texts: Is the base from the first point of the vertex to the last one?


texts: Subtracting rasters after gradient field calculation in R fails?

texts:
I am attempting to write a code in R that imports 3 rasters with the same extent/projection, calculates the horizontal and vertical gradient fields (rise over run in cells), then subtract the outputs from each other (for example, after horizontal gradient fields are calculated for each of the three rasters, the horizontal gradient output for raster 1 is subtracted from the horizontal gradient output for raster 2). I’m using this code on multiple "triplets" of rasters. Sometimes this code works for some triplets, but not for others. This would suggest the problem is with the input rasters, and not the code, but when I inspect the rasters that cause the code to fail, they appear no different than other rasters that the code runs through successfully.

For rasters when the code does not work, I first get the warning message:

Warning message:
In .rasterFromRasterFile(grdfile, band = band, objecttype, ...) :
size of values file does not match the number of cells (given the
data type)

Then, when I try to subtract the vertical rasters from each other (this problem only happens for subtracting vertical gradient field rasters, not horizontal), I get this error message:

Error in setValues(r, methods::callGeneric(getValues(x))) :
length(values) is not equal to ncell(x), or to 1

It’s important to note I only get this error when trying to find the absolute value of the differences; if I just subtract the two rasters I get the warning, but not the error message.

Below is the code I have written for this.

library(raster)

#Load rasters
raster1<- raster(*input path to raster 1*)
raster2<- raster(*input path to raster 2*)
raster3<- raster(*input path to raster 3*)

#Calculate vertical (i) and horizontal (j) gradient fields of raster
#files. Warning message occurs for calculation of vertical gradient
#fields, but not horizontal.

zone1_igrad1<- focal(zone1_raster11, matrix(c(-1/2,0,1/2), ncol = 3))
zone1_jgrad1<- focal(zone1_raster11, matrix(c(-1/2,0,1/2), nrow = 3))
zone1_igrad2<- focal(zone1_raster21, matrix(c(-1/2,0,1/2), ncol = 3))
zone1_jgrad2<- focal(zone1_raster21, matrix(c(-1/2,0,1/2), nrow = 3))
zone1_igrad3<- focal(zone1_raster31, matrix(c(-1/2,0,1/2), ncol = 3))
zone1_jgrad3<- focal(zone1_raster31, matrix(c(-1/2,0,1/2), nrow = 3))

#Subtract Horizontal Gradient Field Outputs.
zone1_hori1<- abs(zone1_jgrad1-zone1_jgrad2)
zone1_hori2<- abs(zone1_jgrad1-zone1_jgrad3)
zone1_hori3<- abs(zone1_jgrad2-zone1_jgrad3)
zone1_hori_field<- zone1_hori1 + zone1_hori2 + zone1_hori3

#Subtract Vertical Gradient Field Outputs. (Error message occurs
#here)
zone1_vert1<- abs(zone1_igrad1-zone1_igrad2)
zone1_vert2<- abs(zone1_igrad1-zone1_igrad3)
zone1_vert3<- abs(zone1_igrad2-zone1_igrad3)
zone1_vert_field<- zone1_vert1 + zone1_vert2 + zone1_vert3


texts: How to write a vector layer with defined WKB type

texts:
I would like to save a selection of a polygon vector layer to a newly created line vector layer. Manually, this can be done by: Right-Click on vector layer > Save As > Geometry Type > Line String.

How would I get the same result with PyQGIS?

Until now this works, but creates a polygon vector layer (same as input):

QgsVectorFileWriter.writeAsVectorFormat(vlayer, lines, "utf-8", None, "ESRI Shapefile", 1)

From https://qgis.org/api/classQgsVectorFileWriter.html, I don't understand how or if it is possible to overwrite the output WKB type to MultiLineString. I tried several versions of:

QgsVectorFileWriter.writeAsVectorFormat(vlayer, lines, "utf-8", None, "ESRI Shapefile", 1, Qgis.WKBMultiLineString)

Does anyone know how to do this?


texts: Calculate the frequency of classes of neighbouring raster cells

texts:
I'm struggling with a problem when working with a raster dataset of biotope types:

For each cell I would like to find out the frequency of classes in the adjacent cells, e.g. from the 8 neighbouring cells 3 are of biotope type A, 4 of type B and 1 cell of type C.

I was recommended to try out 'Block Statistics' and work with a moving window, but so without any success.


texts:
I've been trying to use a ogr2ogr -sql call to create a ESRI.shp file with the joined attributes - but it keeps on failing due to

  "ERROR 1 SELECT from table inshape failed, no such table/feature"

My code is:

ogr2ogr -sql "select inshape.*, joincsv.* from inshape left join 'D:\testfolder\attributes.csv'.joincsv on inshape.ID = joincsv.ID" shape_join.shp D:\testfolder\europe_shapes.shp

The .shp and the .csv may have some ID's that not apear in the other table,
I wish to join them so that any .csv IDs that don't exist in the .shp file are droped, while the new colums of .shp IDs that don't exist in the .csv are filled with "NULL".

To be honest I don't quite understand the syntax:

Why is the .csv source defined within the sql statement while the
input .shp and shape_join.shp outside of the SQL statement?
Is shape_join.shp the output file/ path?


texts: Creating/digitizing roads using ArcMap editor?

texts:
This is probably very basic but I am not able to figure out. I am trying to digitize roads by drawing lines. I go to Editor>Star Editing, select the layer, select my layer under 'create features'. Then under 'construction tools', polygon is selected by default but there is no lines or road option. When I try to create a road, it creates a polygon.

How do I create a road/line?


texts: Intersection with closest point using PostGIS?

texts:
I'm trying to update the field of a point table with a value from a linestring one :

update object_elevation
set contour_elevation = c.contour
from contours c
where st_intersects(st_ClosestPoint(c.geom,object_elevation.geom),c.geom)

So what I want to do going from this query is update the contour_elevation field from the table object_elevation with the value stored in the field contour of the layer contours.

So what the query is supposed to do is  return the closespoint to the linestring and then get the value of the elevation from the line layer and update with it the point layer.

However the result returned is the same. All the features have the same value of the field after update and I know very well that shouldn't be like this.

What am I missing in this query?


texts: How to find all connected edges from gas device feature using C# API call to ArcGIS Pro SDK?

texts:
I am currently using connected trace to trace all the specified gas devices (using its Asset Group and Asset Type) to get output from the trace.

I want to find all the connected edge features from traced devices found in the trace output.

I am using ArcGIS Pro SDK to perform Utility Network-> Connected Trace to get the output as traced gas devices specified by a user (using its Asset Group and Asset Type).


texts: Accessing the values in a list in PyQGIS

texts:
I have got this code in PyQGIS to to save a selected feature into a variable.
class selectTool(QgsMapToolIdentifyFeature):

    def __init__(self, iface):
        self.iface = iface
        self.canvas = self.iface.mapCanvas()
        self.layer = self.iface.activeLayer()
        QgsMapToolIdentifyFeature.__init__(self, self.canvas, self.layer)
        self.iface.currentLayerChanged.connect(self.active_changed)

    def active_changed(self, layer):
        self.layer.removeSelection()
        if isinstance(layer, QgsVectorLayer) and layer.isSpatial():
            self.layer = layer
            self.setLayer(self.layer)

    def canvasPressEvent(self, event):
        found_feature = self.identify(event.x(), event.y(), [self.layer], QgsMapToolIdentify.LayerSelection)
        self.layer.selectByIds([f.mFeature.id() for f in found_feature])

        print (type(found_feature))
        print (found_feature)

    def deactivate(self):
        self.layer.removeSelection()

t = selectTool(iface)
iface.mapCanvas().setMapTool(t)

This is working so far. Just the print statement... print (found_feature) ... gives me:
[<qgis._gui.QgsMapToolIdentify.IdentifyResult object at 0x000001E2669BDDC8>]

My question is now: how can I print the features and values which are inside the created list? I guess I somehow have to make it accessible?
I am very new to PyQGIS and also Python.


texts: Delete duplicate label in GeoServer Style

texts:
I created complex label style in QGIS 3.8.1 and saved it to SLD file then uploaded it in GeoServer 2.15 as a new style, but I should add centroid function to avoid duplicate labels, I tried to add it but unfortunately I got error message:

line 43: cvc-complex-type.2.4.a: Invalid content was found starting with element 'se:Geometry'.
One of '{"http://www.opengis.net/se":Displacement,
"http://www.opengis.net/se":PerpendicularOffset,
"http://www.opengis.net/se":VendorOption}' is expected.

Update:

Download SLD file

But no duplicate labels, when i open the layer from "Layer preview" by open layer format.


texts: How to populate an item's tables list for a ArcGIS API for Python Map Service Item?

texts:
I am trying to obtain a specific table from a Map Service item, however the table list is empty.  The layers list correctly reports back 5 layers.

My code snippet below:

    search_result = gis.content.search(mapServiceName,'Map Service')

    if len(search_result) == 0:
        print ('Map Service: "%s" not found.' %  mapServiceName)
        return

    #take the first return from the list (should only be 1)
    item = search_result[0]

    if not item is None:
        print ('Found map service layer: %s (id=%s)' % (item.title, item.id))

        print ('There are %s tables in item : %s' % (len(item.tables), item.title))
        print ('There are %s layers in item : %s' % (len(item.layers), item.title))

Found map service layer: xxx (id=2358b6294549491da95820a679304bc0)
There are 0 tables in item : xxx
There are 5 layers in item : xxx

The Map Service itself has 2 tables and 5 layers:

I can access those 5 layers with no problem, but the "ArcGIS API for Python" doesn't list the tables.  How to populate that table list?


texts: Are there any cartography "rules" you should follow when labeling maps?

texts:
We've all seen poorly labeled maps before. In these maps, the wording is either illegible or unpleasant to view due to positioning, color choice or size, or the labels take up so much map space that your map appears "crowded". A poorly labeled map can signal unprofessionalism and averts attention from the message the map is attempting to convey.

What sort of efficient labeling strategies can we as GIS users utilize to make our maps clearer, more efficient and all around better for the end user? Are there any hard-and-fast "do not do this" rules to follow when it comes to labeling your maps?


texts:
I have a problem with OSM rendering stack (mapnik, renderd, mod_tile) + osm_carto style and planet.osm loaded in Postgres.

After the expiration period, tiles were updated and now, some of the country borders  and country labels disappeared.

I have manually re-rendered everything using render_list with force parameter for levels 1-5 and it still looks the same. Do you have any idea why did this happen? Maybe failed db updates?


texts: Using Get function jQuery calling Carto data on Leaflet?

texts:
I have problems with a get call for populate an empty geoJson. My purpose is the following:

I want to create a geometry. (Working)

Write into carto database. (Working)

Get the recent geometry and put it in a empty geojson.

Populate another geojson with the data recently got it.

What it happens is not all the reponses of the get call are data, some are empty. Then it cannot be written in the geojson. Others the get call comes with data but is not recognized o something similar.

At the first call the reply can come empty, until third or fourth attempt, when the call comes as expected, with data, but still the geometry is not written in the geojson. Is in the following atempt of geometry record when it starts to create the previous one. And then, sometimes the response can come empty.

//Function to get data to Carto
function cogeDatos(usuarioCarto, sql, geojsonLayer){
    $.getJSON("https://"+usuarioCarto+".carto.com/api/v2/sql?format=GEOJSON&q="+sql+"", function(data) {geojsonLayer.addData(data)});
}

escribeCarto(sqlN); //Function that writes the record in the carto Data Base is working nicely, so no need to show.
var geojsonTemp = geoJsonTemplate; //Temporal GeoJson based on an geojson object
var sqlTemp = "SELECT cartodb_id,the_geom, name from suc_zgz where name='"+nomCab+"';";//SQL query to get again the recently record
cogeDatos(userCarto, sqlTemp, geojsonTemp);//getting the new row
var nFeat = geojsonTemp.getLayers().length;//Finding the amount of new geometries
console.log('Numero elementos a introducir '+nFeat);
geojsonTemp.eachLayer(function (layer) {
                                        console.log(layer.feature.properties.name);//Trying to find the value of the new Data added to the temp geoJson
                                       });
geojsonSuc.addData(geojsonTemp.toGeoJSON());//adding the data to a previous geoJson layer.

The problem really comes with the get function cogeDatos, I don't know if the problem comes because of me, it's a kind of bug of Carto, jQuery.

EDIT: You can see the whole code working here. Now, some responses are randomly empty. So I can not add any geometry to the geojson when this happens. I have updated the fiddle, now less code it's on it. I hope all the problems are repaired.


texts: Creating polygons with different heights in Qgis2threejs

texts:
i'm trying to display a polygon in qgis-2threejs which works fine. Now i want to adjust the height on the sides of the polygon, so that i can vary the body of the polygon.

Basically i want to display different fractures inside of a mountain with the height of each polygon adjusted to the surface.

See here:


texts: running geoserver as background process in linux

texts:
I want to run geoserver as a background process:

sudo -b sh /usr/share/geoserver-2.8.3/bin/startup.sh

However, whenever I interact with the geoserver on the web surface, the process is suddenly re-appearing in the terminal an blocking further interaction. How do I prevent that?


texts: Getting “Cannot zoom to selected feature(s) : No extent could be determined” from PyQGIS

texts:
I am using QGIS 3.4.4.

I am trying to zoom on a specific place of my map using PyQGIS

In order to do so I've defined the layer I wanted to get the features :

couche_parcelles = QgsProject.instance().mapLayersByName('Parcelles')[0]

Then I defined the extent of my layer in order to zoom :

ext = couche_parcelles.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()
zoomRectangle = QgsRectangle(xmin,ymin,xmax,ymax)

Then I created the canvas in order to have a reference to my layer :

canvas = qgis.utils.iface.mapCanvas()
canvas.setExtent(zoomRectangle)
canvas.refresh()

Then I get the Id's I need to get the zoom on the specific point I'm aiming at :

identifiants_parcelles = [x['Idu'] for x in couche_parcelles.getFeatures()]

Finally, I write the block which is supposed to make me zoom on my section :

expr = "\"Id\"=36773"
selection = QgsExpression(expr)
vecteur = couche_parcelles.getFeatures(QgsFeatureRequest(selection))
elemts = [x for x in vecteur]
couche_parcelles.selectByIds(elemts)
canvas.zoomToSelected(couche_parcelles)
canvas.refresh()

It did work once without the need to specify the extent of my layer but then I've had the same error : Cannot zoom to selected feature(s): No extend could be determined

and yet I did determine my extent here so I don't understand why QGIS doesn't perform the zoom on the Id I'm giving him


texts:
var bandname = ee.String(ee.ImageCollection.first().bandNames())

gets a print of:

using

ee.Image.regexpRename(bandname, new_name)

or using

var bandname = ee.String(bandname)

results in :

print(band_name[0],'bandname')

returns

How do I call/convert the Bandname as/to String?


texts: Selecting multiple polygons based on common field in a point file (point does not fall within polygon) to perform service area analysis

texts:
I have 2 layers; a point layer defining poll locations (340 points), and a polygon layer defining poll divisions (areas assigned to a poll location) (2079 polygons). Each point relates to one or more polygon. They have a common point ID field. I would like to perform a service area analysis to determine if each point is within a 25 min drive of its associated polygon/polygons.  I would like to use a model to automate the process but I don’t know how to automate selecting multiple polygons for a single point in a related table.
Below is a selected study area with Polygon attributes displayed and related point features selected.
I am using ArcGIS 10.2.2 for Desktop.


texts: Calculating distances from points to centre of cells using QGIS?

texts:
I am using QGIS 1.8.0 (using older version because I need the Qmarxan plugin which is not available for 2.4.0).  I have a study area that has a vector laying of "planning units", or cells over the top (several thousand cells). Each planning unit has an individual ID and a couple of other fields in the attribute table. I then have a vector point layer with 8 points (ranger stations).

I need to calculate the distances from the centre of every cell to the nearest point, and have this as a new column in the attribute table.  Later, I will need to do the same calculation but from the centre of each cell to the nearest road.


texts: How can I hide custom model from mytoolbox in ArcGIS

texts:
I am created model in ArcGIS using model builder and python script. To build main model I have created some submodel and script in the same toolbox. I just want to show main model in the toolbox and other some model and script I want to hide so the user can not see them.

Please let me know how can achieve the same?

Thanks in advance.

Dharmesh


texts: Have editable label locations (dynamic) on a layer on which there is query

texts:
Is it possible to have dynamic labels, set by editable x y parameters in the attribute table, on a layer which also has a query applied, in QGIS?

i.e. I want to set the query and get all the peripheral data off the map canvas, then align the remaining labels with the "move label" tools in QGIS 1.8.

If I create a second shapefile just for labelling, later on when the master dataset has edits, I end up loosing the x y coords established for the labels.


texts: OpenLayers attribute substitution failing on KML file

texts:
I'm trying to read in mile markers for a marathon from a KML file and use the MILE attribute to set the label.  My layer code looks like this:

// add layer
mile_markers_layer = new OpenLayers.Layer.Vector('Mile Markers',
   {
    projection: new OpenLayers.Projection('EPSG:4326'),
    protocol: new OpenLayers.Protocol.HTTP({
        url: 'Mile_Markers.kml',
        format: new OpenLayers.Format.KML({
            extractAttributes: true
                    //extractStyles:  true
        })
    }),
    strategies: [new OpenLayers.Strategy.Fixed()],
            style: {
                'fillColor': '#888888',
                'fillOpacity': 1,
                'fontColor': '#000000',
                'fontFamily': 'arial, sans-serif',
                'fontSize': '.9em',
                'fontWeight': 'bold',
                'label': '${MILE}',
                'pointRadius': 8,
                'strokeColor': '#ff0000',
                'strokeWidth': 3
            }
});

I looked in Firebug and all the features of the mile_markers_layer have attribute of MILE with a value of "1" thru "26" (as expected).

Layer is displayed with a label of ${MILE} - by this I mean all 26 mile markers are displayed and each has a label of ${MILE}

I was expecting the value of MILE instead.  I'm a newbie, so I could be doing something extremely stupid.

The offending KML file is at:

https://docs.google.com/open?id=0BxhE8rRos0u1X2ExNldTdWwwWkk

Thanks for the suggestion, but every mile marker ended up with ${MILE.value} for a label.


texts:
Why does the MXD size increase over time when you save it?

Obviously you drop data in the map while doing work but then you delete it (temporary clipped shps, tiffs, excel/database).


texts:
I am looking for QGIS function to extract (or detect) edge from a raster, which was referred in here, How to digitize a raster automatically as answer by @MappaGnosis.

The name could have changed since then, or it may have been integrated into other tool... I cannot find it in my QGIS Processing Toolbox (2.14.10 and 2.18.2).

There are similar tools such as GRASS "i.zc Zero-crossing" and SAGA "Gaussian filter", but I am hoping to find this QGIS geoalgorithm.


texts: ArcPy: Error in For loop panToExtent + export to Tiff

texts:
I'm trying to trouble shoot a script, which has the objective of allowing a user to input a query for selected field, use the panToExtent function to pan the data frame at the specified scale to that feature, export as tiff, and then continue until the end of the features queried. This is the code I wrote:

import arcpy

#Set Overwrite Option
arcpy.env.overwriteOutput = True

mxd = arcpy.mapping.MapDocument("C:\\temp\\Graffiti_DM1.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
lyr= arcpy.mapping.ListLayers(mxd, "*", df)[0]

#USER INPUTS
lyr1="C:\temp\itscomplicated.shp"
query="Incident_Z" =="10006"
expdir="C:\temp\df"

for df in arcpy.mapping.ListDataFrames(query):
    df.rotation = 0
    df.scale = 1000
    outFile = expdir + df.name + ".tif"
    df.panToExtent(lyr.getSelectedExtent())
    arcpy.mapping.ExportToTIFF(mxd, outFile, df)
del mxd

Here is the traceback error:

Traceback (most recent call last):
  File "C:\temp\pantoextent.py", line 23, in <module>
    for df in arcpy.mapping.ListDataFrames(query):
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\mapping.py", line 1478, in ListDataFrames
    result = mixins.MapDocumentMixin(map_document).listDataFrames(wildcard)
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\mixins.py", line 727, in listDataFrames
    return list(reversed(list(self.dataFrames)))
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\mixins.py", line 694, in dataFrames
    return map(convertArcObjectToPythonObject, self.pageLayout.dataFrames)
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\mixins.py", line 678, in pageLayout
    return convertArcObjectToPythonObject(self._mxd._arc_object.pageLayout)
AttributeError: 'bool' object has no attribute '_arc_object'

Any help would be MUCH appreciated!


texts: What is this giving me? Trying to sum MODIS ET data by each month

texts:
I am currently evaluating water resource availability over the Bahamas using GEE (I'm new to this) and am running into some trouble. I am trying to sum the ET over the Bahamas for each month from MODIS and am getting some odd numbers. I want to make clear, I need the sum from each month in the time frame. I have done this for precipitation, but the numbers do not appear to be as "off" as the ET values. Below is my code:

//SETUP DATES
var startYear = 2001;
var endYear = 2019;
var startDate=ee.Date.fromYMD(startYear,01,01);
var endDate=ee.Date.fromYMD(endYear+1,12,31);
var months = ee.List.sequence(1,12);
//SET REGION
var dataset = ee.FeatureCollection("USDOS/LSIB/2013").filter(ee.Filter.inList('name', ['THE BAHAMAS']));
var region = dataset
var styling = {color: 'red', fillColor: '00000000'};
Map.addLayer(region.style(styling))

//SET VARIABLES
var sst = ee.Image('NOAA/AVHRR_Pathfinder_V52_L3/20120802025048').select('sea_surface_temperature')
  .rename('sst')
  .divide(100);
var P = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY").select('precipitation')
            .filterDate(startDate, endDate);
var ET = ee.ImageCollection('MODIS/006/MOD16A2').filterDate(startDate, endDate).select(['ET']);

//MONTH
//ET
var nMonths = ee.Number(endDate.difference(startDate,'month')).round();
var byMonth = ee.ImageCollection(
  ee.List.sequence(0,nMonths).map(function (n) {
    var ini = startDate.advance(n,'month');
    var end = ini.advance(1,'month');
    return ET.filterDate(ini,end)
                .select(0).sum()
                .set('system:time_start', ini);
}));

print(byMonth);

Map.addLayer(ee.Image(byMonth.first()),{min: 15, max: 35},'ET');

// plot full time series
print(
  ui.Chart.image.series({
    imageCollection: byMonth,
    region: region,
    reducer: ee.Reducer.mean(),
    scale: 1000
  }).setOptions({title: 'ET over time'})
);

I'm getting something like 200 mm of ET in a month....that can't be correct. That results in close to 3000 mm of ET in a year, and based on precipitation, that's impossible. I am trying to figure out if there is a calculation I am missing or something in the ET product I am unaware of. I have tried returning ET as a sum and a mean (see ET.filterDate(ini,end).select(0).sum()) and get not so great values for both. Am I missing something regarding what MODIS is giving me? The units are in kg/m2/8days which is mmH2O/8 days.


texts:
How do I generate geometric OGC filters for use in a WFS GetFeature request such as:

<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <ogc:Intersects>
        <ogc:PropertyName>geom</ogc:PropertyName>
        <gml:Box xmlns:gml="http://www.opengis.net/gml">
            <gml:coordinates decimal="." cs="," ts=" ">1298543.0367787,7107400.531374 1299312.4070614,7203464.9016567</gml:coordinates>
        </gml:Box>
    </ogc:Intersects>
</ogc:Filter>

The OpenLayers 2 supported geometric filters (e.g., Intersections). How do I write a filter in the OL3?
ol.format.ogc.filter seems does not support (geometric filters).


texts:
I am trying to remove anything in a map that has "DI" in the source path.

I have a test mxd with one df and 2 lyr's. The lyr data sources are:

  C:\Users\Corey.Fields\Python\Testing\DI_Remove\Fake_DI_Wells.shp
  S:\Departments\GIS\ARCUSERS\PERSONAL_FOLDERS\CFields\Test.gdb\Fake_DI_Test

My code is:

import arcpy
mxd = arcpy.mapping.MapDocument("C:\Users\Corey.Fields\Python\Testing\DI_Remove\Untitled.mxd")
for df in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
        if lyr.dataSource == "*DI*":
            arcpy.mapping.RemoveLayer(df, lyr)
mxd.save()
del mxd

The script is running but not removing the layers.

What do I have wrong?


texts: QgsMapToolEmitPoint.canvasReleaseEvent not firing with SHIFT modifier

texts:
I've created a plugin that implements a Map Tool. It captures mouse click events, which works fine. Now  I want to capture Shift+Mouseclick. This doesn't work at all, while Alt+Mouseclick and Ctrl+Mouseclick do work. What's causing this?

The way I'm capturing the events is as follows:

class GeometryTool(QgsMapToolEmitPoint):
    def canvasReleaseEvent(self, event):
        QgsMessageLog.logMessage("Click: " + str(event.button()))


texts:
I'm about creating a python plugin for my project , I already compiled it but I'm now facing a real problem while the installation , it seems to me that  there is issue about the version but I'm not sure ... this is the error O get :

Impossible de charger l'extension dataexporter provoque une erreur lors de l'appel de sa méthode classFactory()

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISBR~1/apps/qgis/./python\qgis\utils.py", line 208, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/soukaina/.qgis2/python/plugins\dataexporter__init__.py", line 34, in classFactory
    from .dataexporter import dataexporter
ImportError: cannot import name dataexporter

Version de Python :
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

Version de QGIS :
2.6.0-Brighton Brighton, exported

Chemin vers Python : ['C:/PROGRA~1/QGISBR~1/apps/qgis/./python', u'C:/Users/soukaina/.qgis2/python', u'C:/Users/soukaina/.qgis2/python/plugins', 'C:/PROGRA~1/QGISBR~1/apps/qgis/./python/plugins', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\bin\python27.zip', 'C:\PROGRA~1\QGISBR~1\apps\Python27\DLLs', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\plat-win', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISBR~1\bin', 'C:\PROGRA~1\QGISBR~1\apps\Python27', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', 'C:\PROGRA~1\QGISBR~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg']


texts: Point otbApplication for python to the relevant OTB DLL's

texts:
I installed the latest OTB and Monteverdi on a Win 7 machine. Via osgeo4w I downloaded otbApplication for python, without any dependencies.

Python, gdal, numpy scipy and so forth was custom build from source and using it with the osgeo python and related dependencies "break" the custom build python and modules.

I copied the three otbApplication files from the osgeo4w install location to the custom python install location (C:\Python27\Lib). When importing otbApplication, an Import error is raised stating

>>> import otbApplication

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
    import otbApplication
  File "C:\Program Files\Monteverdi-3.0\lib\orfeotoolbox\python\otbApplication.py", line 28, in <module>
    _otbApplication = swig_import_helper()
  File "C:\Program Files\Monteverdi-3.0\lib\orfeotoolbox\python\otbApplication.py", line 24, in swig_import_helper
    _mod = imp.load_module('_otbApplication', fp, pathname, description)
ImportError: DLL load failed: The specified module could not be found.
>>>

My question is:
How do I go about pointing the otbApplication for python to the relevant OTB DLL's in a windows 7 environment?


texts: Which conditions does the geoJSON file read from the server in?

texts:
I've added an alert to the function myReadFunction in order to detect function calls.

As according to my observations, the function readMyFunction is called sometimes while I have been navigating on the map, not every time. The function readMyFunction has not been called when a feature is being viewable or not during map navigation. Could you please explain in which conditions the geoJSON file is read from the server?

 mVectorLayer = new OpenLayers.Layer.Vector("Overlay", {
        strategies: [new OpenLayers.Strategy.BBOX()],
        protocol: new OpenLayers.Protocol.HTTP({
            url: 'http://localhost:56786/jlist.geojson',
            format: new OpenLayers.Format.GeoJSON({
                'read': myReadFunction,
                'internalProjection': map.baseLayer.projection,
                'externalProjection': new OpenLayers.Projection("EPSG:4326")
            })
        }),
        projection: new OpenLayers.Projection("EPSG:900913")
    });

Thanks a lot for your helps and clarifications


texts: How to avoid QGIS layout from including transparent layers?

texts:
I have a number of shapefile layers loaded with 0% opacity so they cannot be seen on the map canvas, only 1 layer is visible with 100% opacity. But when I export this as a pdf from the layout, it shows all layers as being visible.
Is there a setting somewhere to fix this? I'm using QGIS 3.16.


texts:
I am working with a colleague who is learning QGIS and has a question about the R plugin. His library path (.libPaths()) has two folders:

"C:/Users/RDCERWAW/Documents/R/win-library/3.0"
"C:/Program Files/R/R-3.0.3/library"

When he downloads packages, they are downloaded to his personal library because the library in his Program Files folder is not writeable and he is unable to change permissions (dang Windows!). At any rate, everything works fine within R not run under QGIS, but when he is using R within QGIS, R is only trying to access packages in C:/Program Files/R/R-3.0.3/library. Since he can't write to this folder, the packages that he needs to run are not available.


texts: Trouble with opening/using i.atcorr in GRASS

texts:
I'm relatively new to GRASS. I'm running version 6.4.3 on a Windows 7 machine.  I'm trying to process Landsat 5 imagery from multiple states, so I was hoping to approximately follow the workflow here:
Is there a way to automate the preprocessing of LANDSAT data?

I have no problem running i.landsat.toar, but when I try to use i.atcorr I get errors:

When I start to type into the command console, usually I type the command, then a space, the enter and it opens the GUI window.

When I type i.atcorr, as soon as I press the space, this displays:

  Traceback (most recent call last):
   File "C:\Program Files\GRASS GIS
  6.4.3\etc\wxpython\gui_core\prompt.py", line 1133, in
  OnKeyPressed

  self.cmdDesc = gtask.parse_interface(GetRealCmd(cmd))
   File "C:\Program Files\GRASS GIS
  6.4.3\etc\python\grass\script\task.py", line 470, in
  parse_interface

  tree = etree.fromstring(get_interface_description(name))
   File "C:\Program Files\GRASS GIS
  6.4.3\etc\python\grass\script\task.py", line 448, in
  get_interface_description

  "\n\nDetails: %(det)s") % { 'cmd' : cmd, 'det' :
  decode(cmderr) }
  grass.script.core
  .
  ScriptError
  :
  Unable to fetch interface description for command
  'i.atcorr'.
  Details:

[Note that no details appear after it says "Details:"]

Then if I press "Enter", an error pops up with the same text:

  "Unable to fetch interface description for command
  'i.atcorr'.
  Details: "

I thought perhaps the problem was that I didn't actually have the function installed on my computer, so I went through g.extension and get this:

  g.extension extension=i.atcorr svnurl=http://svn.osgeo.org/grass/grass-addons/grass6
  which: wget: unknown command
  WARNING: The GRASS program path: [C:/Program Files/GRASS GIS 6.4.3] contains a space. The >build which follows will likely fail.
  Fetching  from GRASS-Addons SVN (be patient)...
  svn: E170000: URL 'http://svn.osgeo.org/grass/grass-
  addons/grass6/imagery/i.atcorr' doesn't exist
  ERROR: GRASS Add-on i.atcorr not found in repository or no network connection or another >problem"

Google tells me that I'm not alone in this "contains a space" error, but I couldn't find any clear information on how to solve this or work around it.

Any insight or advice would be much appreciated.  I am not sure if I can circumvent the problem by just using the command console, but as I looked through the wiki on the topic (http://grasswiki.osgeo.org/wiki/Atmospheric_correction) it wasn't clear to me if I was supposed to build my own icdn text file?  That's not something that came when I downloaded the data from USGS. I was hoping the GUI version of the function would make this more clear to me but since I can't get it to open, I'm still confused.

Thanks in advance.


texts:
I have a shapefile containing many points.

Some of these points fall within an area for which I have raster data.

I would like to load the values of the raster which underlie each point into the shapefile such that they appear as a column in its DBF.

Is there a way to do this with, say, GDAL?


texts: I have posted a question about the use of sDNA in ArcGIS. It is said the problem is caused by non-ascii characters. I tired to modify the shapefile and kept the Object ID and Shape length only. However, the analysis still failed. Here is the attribute table of the data and the steps I tried.


texts: Get all Nodes from a List of adjacent LineStrings in PostGIS

texts:
I use PostGis and pgRouting.
I have a table with geometries of the adjacent edges returned by a shortest path query. I want to concatenate these edges and retrieve a list of all nodes these edges consist of. How do I do this in SQL?

I think the solution would be a mix of LineMerge and GeomUnion, but is there a smart way to do it?

Thanks.


texts: Group sublayers in WMS service to a single layer

texts:
I created a WMS service in QGIS that consists of several sub layers. It is a base map of cadastral data that is stored in different tables to get the final look. In the QGIS project for the WMS I grouped the sublayers.
If I add the WMS to another QGIS project, I can still see all sublayers. Is there a way to make only the group layer visible / exposed via WMS?


texts:
I have been trying to find a marker animation plugin for OpenLayers to sync lat, lon, and time data along a track like this (in Leaflet):

https://github.com/hallahan/LeafletPlayback

http://leafletplayback.theoutpost.io/examples/example_1.html

Is it possible to animate a marker in OpenLayers, given the way it draws on the canvas?


texts:
I'm new to GDAL python API, when I tried to rotate a TRMM image, I got a weird result compared to the output of ENVI. I just want to rotate the image with 270 degrees, the result is also flipped vertically.

the output of ENVI desktop

the result with using gdal

import math
from osgeo import gdal
import numpy as np
from osgeo import osr

# open the file
precip = gdal.Open('HDF4_SDS:UNKNOWN:"/home/wanghl/wxq/3B43.19980101.7.HDF":0')

# get the data of the precipitation dataset
preciparr = precip.ReadAsArray()

# get geotransform
transform = precip.GetGeoTransform()

#rotate

rotation = (math.pi/180) * 270

# set geotif driver
driver = gdal.GetDriverByName('GTiff')

# get x,y dimensions of the map
xs = precip.RasterXSize
ys = precip.RasterYSize

new_geotransform = (transform[0],
                    transform[1],
                    -math.sin(rotation) * transform[1],
                    transform[3],
                    math.sin(rotation) * transform[5],
                    transform[5])
# set output name
outname = "/home/wanghl/wxq/3B43.19980101.7.tif"

# set projection
target = osr.SpatialReference()

# write dataset to disk
outputDataset = driver.Create(outname, xs, ys, 1, gdal.GDT_Float32)
outputDataset.SetGeoTransform(new_geotransform)
outputDataset.SetProjection(target.ExportToWkt())
outputDataset.GetRasterBand(1).WriteArray(preciparr)
outputDataset.GetRasterBand(1).SetNoDataValue(-9999)
outputDataset = None


texts: How to change Datasource of a Raster in sde using ArcObjects?

texts:
I'm building a C# application that lets the user change their DataSources of all the layers in their MXD.  I have it working on FeatureLayers and LYR files but I'm stuck on how to change Raster layers from one sde database to another sde database.  Below is all I have so far with Rasters.  To change FeatureLayers I used IMapAdmin2 but theure is no option to change Rasters with that interface.  Any ideas or additional code?  Many thanks

UID uid_ras = new UIDClass();
uid_ras.Value = "{D02371C7-35F7-11D2-B1F2-00C04F8EDEFF}"; //IRasterLayer
IEnumLayer raslayers = pMap.get_Layers(uid_ras, true);
raslayers.Reset();

pRas = raslayers.Next() as IRasterLayer; //get the first layer
while(pRas != null)   //Now change each Raster's DataSource to the new Datasource
{
    //I'm unsure what to do from here..
}


texts: Restrict extent and Zoom in OpenLayers3 upto 1 mile

texts:
I need to restrict map to a certain extent plus one mile. How to do that?

here is my code

    map.setView(view);

    map.getView().fit(bbox, map.getSize());

    view = new ol.View({
        projection: 'EPSG:27700',
        center: center,
        extent: bbox,
        resolutions: {/* list of resolutions*/},
        resolution: {/* resolution value*/}
    });

    map.setView(view);

The bbox is bounding box for the extent.


texts: Is it possible to use arcpy.mapping.UpdateLayer on a raster layer?

texts:
Is it possible to set a raster symbology to unique values using arcpy.mapping.UpdateLayer.   Using 10.3 I am receiving an AssertionError

import sys, os, arcpy
sourceLayer = r'D:\Chad\GIS\TEMP\raster_unq_val.lyr'
mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
lyr = arcpy.mapping.ListLayers(mxd, "MapunitRaster_fy16_10m", df)[0]

print lyr.name
print lyr.symbologyType
arcpy.mapping.UpdateLayer(df, lyr, sourceLayer, True)
desc = arcpy.Describe(lyr)
flds = [x.name for x in desc.fields]
for f in flds:
    if f.find('rating') <> -1:
        targetField = f

values = list()

with arcpy.da.SearchCursor(lyr, targetField) as rows:
    for row in rows:
        aVal = round(row[0], 2)
        if not aVal in values:
            values.append(aVal)
values.sort()

lyr.symbology.valueField = targetField
lyr.symbology.addAllValues()
lyr.symbology.classValues = values
lyr.symbology.classDescriptions = values
lyr.symbology.ShowOtherValues = False

arcpy.RefreshActiveView()
arcpy.RefreshTOC()

del values

with the error:

  Runtime error  Traceback (most recent call last):   File "",
  line 9, in    File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
      return fn(*args, **kw)   File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1889, in
  UpdateLayer
      assert isinstance(source_layer, Layer) AssertionError


texts: SVG control not visible on top of Leaflet.js

texts:
I'm using Leaflet 0.7.1 and want to draw a radial menu (like openstreetmap's iD editor) on top using d3.
Some examples I found use Leaflet's overlayPane to append the svg element to:

var svgContainer= d3.select(map.getPanes().overlayPane).append("svg");

On mouseclick, I add the menu and reposition it to the screen xy coords:

map.on('contextmenu', function(e) {
    console.log('contextmenu()');
    var tooltip;
    var center = [e.layerPoint.x, e.layerPoint.y];
    var menu = svgContainer.append('g')
        .attr("class", 'leaflet-zoom-hide')
        .attr('class', 'radial-menu')
        .attr('transform', 'translate(' + center + ')')
        .attr('opacity', 0);
        menu.transition()
            .attr('opacity', 1);

        menu.append('path')
            .attr('class', 'radial-menu-background')
            .attr('d', 'M' + r * Math.sin(a0) + ',' +
                             r * Math.cos(a0) +
                      ' A' + r + ',' + r + ' 0 ' + 0 + ',0 ' +
                             (r * Math.sin(a1) + 1e-3) + ',' +
                             (r * Math.cos(a1) + 1e-3)) // Force positive-length path (#1305)
            .attr('stroke-width', 50)
            .attr('stroke-linecap', 'round');
    });

Somehow, the SVG paths get drawn (this is visible in the Chrome Inspector) but are behind the map object.
When I 'drag' the SVG element just below the body tag in Inspector, I can see the circle paths.

Any ideas on how I can draw SVG 'menu' elements on top of Leaflet?

Thanks!!

See this fiddle for a demo. Right click or hold to add an invisible element.

Ps. I also asked this question on StackOverflow. Sorry for crossposting.


texts: How to create a polygon from the union of the areas formed by a circle running in a line formed by two geographic points?

texts:
How can I create a polygon like in the right side of the picture from a line formed by two geographic points like in the left side of the picture in PostGIS?

I was thinking in calculating the area by two ways:
1) The polygon is the union of the area of all infinite circles with an radius x running inside the line formed by the two geographic points
2) The area of a rectangle centered in the line with a x*2 width united with the two circles of radius x, one centered in the top and the other in the bottom geographic points.
The points are geographic points (like -72.1235 42.3521) and I like the radius of the circle to be informed in meters.


texts:
I have several Rasterstacks (b1_mask, b2_mask,...) of several layers representing only a determined pixel-intensity range. (So point "A" in some layers is "NA" - in others a "value" of that certain range). My goal is to plot (boxplot) 100 random points of this raster stack - to see the variety of that random point in the different layers of my raster stack. I am working in R and my approach would be the following:

setting random points:

library(dismo)

set.seed(1)
random<-randomPoints(b1_mask, n=100)

the next step is difficult for me: As I have 6 spectral bands (b1_mask - b7_mask) my approach was to extract the point values of the 6 bands and save them into a dataframe from where they are plotted as a boxplot. I am not good in for loops so please help me with this part...

for(i in (1:nrow(random))){

  d<-data.frame(extract(b1_mask, random[i,]),
             extract(b2_mask, random[i,]),
             extract(b3_mask, random[i,]),
             extract(b4_mask, random[[i,]]),
             extract(b5_mask, random[[i,]]),
             extract(b7_mask, random[[i,]]),
  )
  par(mfrow=c(4,4)) #plot 16 figures in one plot

  boxplot(d, las=2, main="Reflectance of bands at Point"random[[i]], ylab="Reflectance" xlab ="band 1-7")
}


texts: What spatial analysis did the US Supreme Court review in the Perry v. Perez opinion?

texts:
The Supreme Court recently issued an opinion on Perry v. Perez  which seems like it must have involved evidence based on spatial analysis.  The L.A. Times reports:

  The 9-0 decision set aside a new map of congressional districts drawn
  by a special federal court in San Antonio that gave Latinos and
  Democrats a good chance to win three or possibly four new seats in the
  House of Representatives.

Does anyone know what type of spatial analysis - if any - was used to create the evidence presented?

Update:
Thanks to mkennedy for interpreting this, and also for pointing me in the direction of preclearance.  I see that the Federal Register Vol 76 No. 27 Part III goes into what I was wondering about under the "Retrogressive Effects".  Fortunately, this is very googleable.


texts: Find closed lines inside SpatialLines objects in R

texts:
I need to find which lines are closed lines and which aren't inside a SpatialLines object in R. A closed line mean that the first point is the same as the last point. Is there any tool or function in R to find/check this? I give a reproducible example below.

How can I query the SpatialLines object called spLines to know that the first line (blue) is closed and the second line (green) opened?

# Libraries
library("sp")

# Points
points.x1 <- c(-54.00000, -54.43489, -54.78560, -54.00000, -53.57810, -53.22097, -54.00000)
points.y1 <- c(-40.74859, -40.00000, -38.00000, -37.49485, -38.00000, -40.00000, -40.74859)

points.x2 <- c(-52.53557, -52.00000, -50.00000, -48.00000, -46.40190)
points.y2 <- c(-41.00000, -40.60742, -40.08149, -40.82503, -39.00000)

# Lines
line1 <- Line(coords = cbind(points.x1, points.y1))
line2 <- Line(coords = cbind(points.x2, points.y2))

# SpatialLines
spLines <- SpatialLines(LinesList = list(Lines(slinelist = list(line1, line2), ID = "A")))

# Plot
plot(spLines, cex = 0, xlab = "Longitude", ylab = "Latitude", main = "Plot")
lines(spLines@lines[[1]]@Lines[[1]], col = "#3465A4", lwd = 4)
lines(spLines@lines[[1]]@Lines[[2]], col = "#73D216", lwd = 4)
box()

Note: this is only a simple example with two lines, the answer should work with a SpatialLines object with many lines.


texts: Expression to run two geometry functions in expression builder

texts:
I have feature polygons that sit within and across larger tenancy polygons, I need to transfer the tenancy reference from the tenancy polygons to the feature polygons. To do this I am using AutoFields Expression Builder to automatically fill in the feature polygons reference field as a new polygon is created.

My problem is I can either run the expression to return intersections or within, I can not work out how to do both at the same time. My current expressions are

geomintersects( 'NAT_CAP_Test_Tenancies',  'REFERENCE' )

and

geomwithin( 'NAT_CAP_Test_Tenancies',  'REFERENCE' )

is there a way to combine both?


texts:
There is a very strange problem with my QGIS 3.0 to show some font styles.

For example you can see that font Gill Sans Ultra Bold has 4 styles and QGIS 2.18.19 can see them:

But QGIS 3.0.3 fails to do that:

I've tried to reinstall fonts, reinstall QGIS 3.0.3 but the issue isn't fixed.

(Update) It fails to read even Ms Shell Dlg 2:

Any ideas how to solve this?


texts: How to calculate forest loss in google earth engine

texts:
I am using the method outlined in a previous question - Extracting land cover type in Google Earth Engine? - to identify habitat type and area with a defined polygon. I would now like to take this information from different years to calculate forest cover loss over time. How is this best achieved?
I am using google earth engine for this process.


texts: Calculation of a new field with conditions

texts:
I am struggling with the calculation of a new field that is based on another field.

I have one field with numbers that appear twice ("Field_A") and I need to calculate a new field ("Field_B") that has the value -90 when one number appears the first time and the value 90 when the same number appears the second time.

I am working with ArcGIS 10.6. I guess I need a Python code block.


texts: Spatial joining (nearest) between features that fulfill query?

texts:
I have two point feature classes, one with locations wildfires were discovered, and one where lightning strikes were detected. Both feature classes contain temporal information. For each fire, I want to search for lighting strikes within a 5 km buffer that were recorded on the day that the fire was discovered or any of the three days prior.

I believe the below process is fairly straightforward, but I am having trouble coming up with a proper iterative solution.

Draw a buffer of 5 km
Clip lightning to buffer
Select lightning by querying date range
Use Near to join closest lightning strike attributes to fire

How can I most effectively perform this process iteratively over each point in the fire feature class?


texts: Deleting graphics from Go To XY and Add Labeled Point?

texts: How do I delete a green point with coordinates that I added to my map using Go to x,y / Add Labelled Point?


texts: How to best get Qt Designer for QGIS plugin development on Mac?

texts:
I'm using QGIS 2.6.1 on Mac OSX 10.9.4, trying to create a plugin as described in the PyQGIS Programmer's Guide. I have the plugin 'skeleton' (main widget with header and ok/cancel buttons) running in QGIS, and I want to edit the plugin with Qt Designer.

I haven't installed Qt on its own and I don't really understand if, by having compiled and installed PyQt and QGIS, I already have Qt Designer installed. I know the Qt Designer classes are available in PyQGIS, and when I run:

/usr/sbin/system_profiler SPFrameworksDataType

I find:

QTKit:

  Version: 7.7.3
  Obtained from: Apple
  Last Modified: 1/26/15, 1:23 PM
  Kind: Intel
  64-Bit (Intel): Yes
  Signed by: Software Signing, Apple Code Signing Certification Authority, Apple Root CA
  Get Info String: QTKit 7.7.3, Copyright 2003-2012, Apple Inc.
  Location: /System/Library/Frameworks/QTKit.framework
  Private: No

I really don't want to confuse things by installing a separate version of Qt if I don't have to- but I can't for the life of me figure out how to open Qt Designer or any of the Qt modules without first separately downloading and installing Qt.

Should I download and install Qt or can I already use Qt Designer? If I already have Qt Designer, how can I access it (preferably from command line)?


texts: Google Earth Engine API, convert Javascript to Python

texts:
I have the below code in Javascript using Google Earth Engine Code Editor and it generates a chart and I want this in Python if someone help to convert it.
var roi = ee.Geometry.Polygon([66.169021,33.698883,66.165845,33.698669,66.166102,33.696526,66.164471,33.695598,66.168677,33.695384,66.170308,33.697812]);

// Import the Landsat 8 TOA image collection.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');

var withNDVI = l8.map(function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
});

// Create a chart.
var chart = ui.Chart.image.doySeriesByYear({
  imageCollection: withNDVI.filter(ee.Filter.date('2018-01-01', '2021-01-01')).select('NDVI'),
  bandName:'NDVI',
  region: roi
}).setOptions({title: 'NDVI over time'});

// Display the chart in the console.
print(chart);


texts: Drawing lines between all possible combinations of points using ArcGIS Desktop?

texts:
I have a shapefile with several points (A, B, C, D, E,...) and I want to draw a line between all possible combinations of points (A-B, A-C, A-D, B-C, B-D,...).

As I know that it is quite simple in ArcGIS I don't want to use a Python-script or something similar. The "Points To Line" tool just draws me one line that is connecting all the points (A-B-C-D-...).

How can I do it using ArcGIS Desktop (all extensions available) ?


texts: Querying for extent of layer in arcpy.mp?

texts:
We're migrating automated mapping scripts from ArcGIS Desktop to ArcGIS Pro — Python 2.7 to 3, arcpy.mapping module to arcpy.mp module. All these scripts iterate through a loop to select various features of interest, and then focus the dataframe extent on that selection. For example:

gas_main_layer = arcpy.mapping.ListLayers(mxd, "GasMain", data_frame)[0]
data_frame.extent = gas_main_layer.getSelectedExtent()
# adjust scale slightly so view is larger than polygon extent
data_frame.scale *= 1.2

I easily figured out how to transition ListLayers to the new mapping module's workflow, but I can't find a comparable method on the Layer class to get the selected extent and then reset the zoom/extent of its Map. What is the way to get the same functionality out of the new module?


texts: Getting hazed figure when exporting a high resolution image in QGIS

texts:
I have a raster map made in QGIS. I am trying to export that map in .tiff format in high resolution (Somewhere around I can get a file of 300 MB).

So I am exporting the image with a resolution of 960 dpi (using the option Export Map to Image...

But all I am getting is a hazed figure like this:

What could be the possible solution for this issue?

At 600 dpi, I am getting an image like this on exporting:


texts: Create sequence numbers based on location of line features using Arcobjects

texts:
I want to add a Sequence number property for line features (Pipe), which will hold a value that will determine where the line is in respect to other lines that are attached to the from point and to point.

Here is a picture that will provide a better explanation.

So based on the picture above, AB1 sequence would be 0, AB2 sequence would be 1, and AB3 sequence would be 2

This is the class presenting the lines

  internal Pipe(IFeature pipe, Ipoint fromPnt, IPoint toPoint){
      Shape=pipe.ShapeCopy;
      PipeOid=pipe.Oid;
      FromPoint=fromPnt;
      ToPoint=toPoint;
    }

    public int PipeOid;
    public int SequenceNumber;
    public IPoint FromPoint;
    public IPoint ToPoint;
    public Polyline Shape

What i have tried so far.
Looping over a list of Pipes, and using QueryNormal from the middle of each line,and passing the length in a positive and a negative value, so I can get results in both directions, in the picture it would be from top and bottom of line. Then I would check if the normal if it touches other Pipes that are connected to the nodes, and based on those results, I would get my sequences.

foreach(pipe in commonPipes){
var pipeShap=pipe.Shape;
 ILine normalLine;
 pipeShap.QueryNormal(esriSegmentExtension.esriNoExtension, 0.5,true,10,normalLine);
 // test for if normalline touches anything,
 pipeShap.QueryNormal(esriSegmentExtension.esriNoExtension, 0.5,true,-10,normalLine);
 // test if it touches anything
 }

I found it rather long to do so, and that the positive and negative lengths are sometimes reversed. Also the length of the normal could fall short of touching the links below or on top.

I tried creating a point in the middle of each line, and then using ProximityOperator, or just plain X,Y but that is also long since the line could be horizontal or vertical.

Any ideas in the Arcobjects that might help me solve this problem ?


texts: OpenLayers 6 ol.source.ImageWMS does not use/follow provided resolutions

texts:
Even though layer source resolutions are as follows

var source = new ol.source.ImageWMS({
    url: this.url,
    params: parameters,
    ratio: 1,
    resolutions: [396.8757937515876, 132.2919312505292, 52.916772500211672, 26.458386250105836, 10.5833545000423344, 3.9687579375158754]
});

Images are being requested for intermediate resolutions, for example when resolution is 19.xxx as reported by map.getView().getResolution()

FYI: map resolutions are different. Providing that range only ..., 66.1459656252646, 26.458386250105836, 13.229193125052918, ...


texts:
I've been reading the "Whats New in ArcGIS" for the ArcGIS 10.3 update and came across the Python and ArcPy portion. It states that:

  ArcGIS 10.3 has been upgraded to include Python 2.7.8.

Now this is fine, I have zero issues with any 2.x version of python. My concern is that Esri will eventually push ArcGIS to Python 3x. I hope they don't in the immediate future, but at the same time Python 3x will need to be implemented at some time.

I do a ton of coding in python for all sorts of GIS and non-GIS applications. I always use the version current with my ArcGIS release. If they move to 3x, I'll be spending a solid month doing 2x to 3x upgrading for my applications.

I'm just wondering if anyone has any insider information out there concerning this topic.

Does Esri have any plans in the near future to move to Python 3 or will they continue on with Python 2x?


texts: How to resolve this "IndentationError: unexpected indent" error?

texts:
Used this code from http://www.qgistutorials.com/en/docs/building_a_python_plugin.html

layers = self.iface.legendInterface().layers() layer_list = [] for layer in layers:
    layer_list.append(layer.name())
    self.dlg.comboBox.addItems(layer_list)

and getting this error after reloading

IndentationError: unexpected indent
Traceback (most recent call last):
  File "C:/OSGeo4W/apps/qgis/./python\qgis\utils.py", line 333, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/vedan/.qgis2/python/plugins\SaveAttributes\__init__.py", line 34, in classFactory
    from .save_attributes import SaveAttributes
  File "C:/OSGeo4W/apps/qgis/./python\qgis\utils.py", line 607, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/vedan/.qgis2/python/plugins\SaveAttributes\save_attributes.py", line 185
    layers = self.iface.legendInterface().layers()
   ^

Using QGis 2.18.3


texts: Filtering noise out from bare-earth DEM with terrace fields?

texts: I have a 1m bare-earth DEM, derived from LiDAR data, representing a terraced landscape. For the sake of any subsequent use of the DEM (e.g., least-cost path, cost surface analysis, and the like) I am wondering if it is sound/necessary to 'blur' the DEM in order to de-emphasize any possible 'noise' produced by the terrace fields (e.g., sharp changes in slope values due to the terrace walls). If the answer is yes, what type of approach is better: something like a low-pass filter (I am not familiar with it, and how it can be performed in ArcGIS 10.1) or should I go back to the LiDAR data and create a bare-earth DEM from scratch with a larger cell size?


texts: Globe not displaying in my cesium application

texts:
The below code doesn't display the default Cesium globe. When I run my file, its just a black background. Why is it so?

Loading...


texts: "L.Proj is undefined" error when using Proj4Leaflet in Angular 5

texts:
I am facing the error "ERROR TypeError: Cannot read property 'CRS' of undefined" caused by L.Proj being undefined.

Here what I did :

I added Leaflet, Proj4 and Proj4leaflet packages using npm commands.
I added @types/leafelt, @types/proj4 and @types/proj4leaflet.
I included required js paths in .angular.json.

"scripts": [
          "../node_modules/jquery/dist/jquery.min.js",
          "../node_modules/bootstrap/dist/js/bootstrap.js",
          "../node_modules/tether/dist/js/tether.min.js",
          "../node_modules/leaflet/dist/leaflet.js",
          "../node_modules/proj4/dist/proj4.js",
          "../node_modules/proj4leaflet/src/proj4leaflet.js"
          ]

In my Typescript file I included the following imports :

import * as L from "leaflet";
import * as proj4 from "proj4";
import * as proj4leaflet from "proj4leaflet";

Finally here is an extract of the code :

var crs = new L.Proj.CRS('EPSG:3006',
    '+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs', {
        resolutions: [
            8192, 4096, 2048, 1024, 512, 256, 128
        ],
        origin: [0, 0],
        bounds: L.bounds([218128.7031, 6126002.9379], [1083427.2970, 7692850.9468])
    });

//32.3340497, -6.3590815
this.mymap = L.map("mymap", {
    crs: crs,
    attributionControl: false
}).setView([32.3340497, -6.3590815], 12);

L.tileLayer("http://api.geosition.com/tile/osm-bright-3006/{z}/{x}/{y}.png", {
    minZoom: 7,
    maxZoom: 18
}).addTo(this.mymap);

Could you please give me any help.


texts: Using PostgreSQL stored procedure with GeoServer?

texts:
I was wondering if it is possible to call PostgreSQL stored functions from GeoServer?

I tried to use GeoServer's WPS to do map algebra, but found that the operations it offered was somewhat limited to things such as gs:AddCoverages. This makes me wonder if it is possible to use PostgreSQL/PostGIS functions and stored functions directly to do arbitrary algebra, e.g.:

SELECT ST_MapAlgebra(t1.rast,t2.rast
        ,'([rast1]+0.001)/([rast2]+0.001)') AS rast
FROM img1 t1, img2 t2;

I googled around, and didn't find much, except a related post here about using stored procedure of Microsoft SQL server.


texts: Unable to merge features in KML with QGIS?


texts: Python not appearing as option for parser in label expression dialog of ArcGIS Desktop 10.0?

texts:
I have ArcMap 10.0 and python is not appearing as a parser option in the label expression dialog box, only VBScript and JScript. Python is installed and the python window works fine.

Any idea how I can get this function?


texts: Is there a "Disperse Markers" tool for QGIS?

texts:
I am using ArcGIS, but either I would like to use QGIS and maybe want to switch completely later.
A have a point overlapping problem, as a result of geocoding of addresses.

I have not found any tool in QGIS which could replace the "Disperse Markers" tool from ArcGIS, because if I understand well, the "Points displacement" tool in QGIS can displace the overlapping polygons only in "circle way" which has finally a very artificial effect.

Please check out these two short articles about Disperse Markers to see what I mean. Please if you know any possible solution in QGIS for this, let me know.

http://blogs.esri.com/esri/arcgis/2008/08/01/how-to-scatter-marker-symbols-for-stacked-or-clustered-point-features/


texts: Leaflet geoJson in front of markers in a FeatureGroup

texts:
I am loading some markers and json data to a leaflet map.

I need to get my json in front of the markers. I have tried placing the markers and json in seperate FeatureGroups and then tried bringToFront() and bringToBack() to no avail

var markerHolder = new L.FeatureGroup().addTo(map).bringToBack();
var jsonHolder = new L.FeatureGroup().addTo(map).bringToFront();

I also tried using this panes work around but when using a Feature group in place of the tileLayer as in the example below, an error was thrown at getContainer()

var topPane = map._createPane('leaflet-top-pane', map.getPanes().mapPane);
    var topLayer = L.mapbox.tileLayer('lxbarth.map-vtt23b1i').addTo(map);
    topPane.appendChild(topLayer.getContainer());
    topLayer.setZIndex(9);

Here is a fiddle illustrating my problem - could anyone demonstrate how to get my json features in front of the markers?


texts:
I am unable to chart a histogram for my region for each image in an image collection.

// Create clipping function
var clipper = function(image){
  return image.clip(geometry);
};

var modis = ee.ImageCollection('MODIS/006/MOD09GQ').map(clipper);
var months = ee.List.sequence(1,12);

// Group by month, and then reduce within groups by mean();
// the result is an ImageCollection with one image for each
// month.

var byMonth = ee.ImageCollection.fromImages(
      months.map(function(m) {
        return modis.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .filterBounds(geometry)
                    .median()
                    .set('month', m);
}));

print("Monthly Images Median",byMonth);

//chart histogram of a single band
print(ui.Chart.image.histogram({
  image: byMonth.first().select(['sur_refl_b01']),
  region: geometry.geometry().bounds,
  // maxBuckets: 10,
  scale: 250
}).setOptions({title: 'modis median band 01'}));

Error:

Unrecognized argument type to convert to a FeatureCollection: function standin(_0) {
    return newF.apply(this, arguments);
  }

Even though I am charting one month at a time, I would ideally like 12 histograms for 12 images in the collection, and a single chart with all 12 histograms drawn.

EDIT I tried the following and I get another error saying unexpected argument type

print(ui.Chart.image.doySeriesByRegion({
  image: byMonth,
  bandName: 'sur_refl_b01',
  regions: geometry,
  // reducer: ee.Reducer.first(),
  scale: 250
}).setOptions({title: 'monthly median sur_refl_b01'}));

Error:

Unexpected arguments to function doySeriesByRegion(): image


texts: CartoDB adding custom legend to sublayer

texts:
Here's what I'm trying to do:

I have a map with selectors that load data from the same table (using SQL) based on what you selected. This all works great (code below). What I want is for a legend to be visible and for it to only show the values included in the current layer. For example, if the layer has only two possible categories, then I want the legend to show only those two categories.

Bonus points for being able to customize the legend (category names) as I create the layer.

I have no idea where to include anything about legends in this code:

function createSelector(layer) {
    var sql = new cartodb.SQL({ user: 'myusername' });
    var $options = $('#layer_selector li');
    $options.click(function(e) {
      // get the area of the selected layer
      var $li = $(e.target);
      var attribute = $li.attr('data');
      // deselect all and select the clicked one
      $options.removeClass('selected');
      $li.addClass('selected');
      // create query based on data from the layer
      var query = "select * from attribute_table";
      if(attribute == 'gender') {
        query = "select *  from attribute_table where "+ attribute + " > '0'";

        new_css = '#attribute_table{marker-width: 4;}#attribute_table[gender=2]{marker-fill: #A6CEE3;}#attribute_table[gender=1] {marker-fill: #1F78B4;}'

      }
        else if(attribute == 'marital') {
       query = "select *  from attribute_table where "+ attribute + " > '0'";

        new_css = '#attribute_table{marker-width: 4;}#attribute_table[marital=1] {marker-fill: #7B00B4;}#attribute_table[marital=2] {marker-fill: #136400;}'

      }
      else if(attribute == 'religion') {
        query = "select *  from attribute_table where "+ attribute + " > '0'";

        new_css = '#attribute_table{marker-width: 4;}#attribute_table[religion=1] {marker-fill: #A6CEE3;}#attribute_table[religion=3] {marker-fill: #1F78B4;}#attribute_table[religion=4] {marker-fill: #B2DF8A;}#attribute_table[religion=5] {marker-fill: #33A02C;}#attribute_table[religion=2] {marker-fill: #FB9A99;}#attribute_table[religion=7] {marker-fill: #E31A1C;}#attribute_table[religion=8] {marker-fill: #FDBF6F;}#attribute_table[religion=6] {marker-fill: #FF7F00;}'

      }
      else if(attribute == 'all') {
        query = "select *  from attribute_table";
        new_css = '#attribute_table{marker-width: 4;}#attribute_table{marker-fill: #2e5387;}'

      }
      // change the query in the layer to update the map
      layer.setSQL(query);
      layer.setCartoCSS(new_css);

    });
  }
  function main() {
    cartodb.createVis('map', 'https://myusername.cartodb.com/api/v2/viz/f6b34e86-b382-11e4-b544-0e853d047bba/viz.json', {
      tiles_loader: true,
      center_lat: 39.8282,
      center_lon: -98.5,
      zoom: 5,
      legends:true,
      fullscreen:true,

    })

    .done(function(vis, layers) {
      var subLayer = layers[1].getSubLayer(0);
      createSelector(subLayer);
    })
    .error(function(err) {
      console.log(err);
    });
  }
  window.onload = main;


texts: RuntimeError: Field is not editable when using InsertCursor on ArcSDE geodatabase with ArcFM installed?

texts:
We are running ArcGIS Desktop 10.1 SP1, and also have ArcFM installed.

I have this simple code:

workspace = "C:\\my.sde"
arcpy.env.workspace = workspace

edit = arcpy.da.Editor(arcpy.env.workspace)
edit.startEditing(False, True)
edit.startOperation()

with arcpy.da.InsertCursor("C:\\my.sde\\dataset\\feature", ("Comments")) as icur:
    try:
        icur.insertRow(["Test comment"])

    except Exception as e:
        print e

edit.stopOperation()
edit.stopEditing(True)

I receive the very specific (/s) error:

  SystemError: error return without exception set

I've tried the same code on a different feature class in the same geodatabase and I get:

  RuntimeError: Field is not editable.

But I can't think of any reason why it's not editable. I can edit the data in ArcMap just fine. Any idea what could be wrong?


texts: Setting vertical offset for SLPK layer in ArcGIS Runtime .NET sdk?

texts:
I recently started development of an app using ArcGIS Runtime .NET SDK (WPF). I have an issue related to the the elevation offset for scene layer package (SLPK).

When I load my SLPK to ArcGIS Online it overlays above the ground, but in ArcGIS Online I can set the elevation offset that corrects the vertical position of the model.

Now if I use the portal item or ArcGIS Online service in .NET WPF app it works fine because my model has been published with an offset.

But I want to load my SLPK file local path, and I am unable to find any such option to set the offset value for ground overlay.

Uri _serviceUri = new System.Uri(@"D:\Projects\test\lod75.slpk");
ArcGISSceneLayer sceneLayer = new ArcGISSceneLayer(_serviceUri);
myScene.OperationalLayers.Add(sceneLayer);


texts: How can I join a spreadsheet to a point layer?

texts:
I have a layer with point data for the grave positions in our cemetery.  We are gathering data from each headstone. If I enter the data to a spreadsheet can I later associate the spreadsheet with the point data layer.

Will I have to go back to the point data layer and add all the columns from the spreadsheet as attributes to make the point layer searchable by data from the headstones?

Would a better solution be to use a PostgreSQL to make a database to associate with the point layer?

Am I going about this the right way?


texts: ArcPy SummarizeWithin not availabe in ArcMap with Advanced license?

texts: When running "arcpy.SummarizeWithin_analysis()" in PyCharm, I get an error saying the module isn't available. I've done "import arcpy" and I have ArcGIS Desktop 10.6 Advanced. It also doesn't show up as a tool in my Statistics toolbox in ArcMap. According to the webpage, it should be available for my license: http://pro.arcgis.com/en/pro-app/tool-reference/analysis/summarize-within.htm#L_


texts:
I have a table with coordinates of some postcodes and I want to display them in ArcMap. The points should be on the street centrelines. When I bring in the csv file and right click to display the points, they are being located on the properties instead of the street centrelines. I have nearly 5000 records in the table so I have to find a automated way to deal with this.

Does anyone knows what tool to run to bring the point onto the centrelines?


texts: Getting distance between route and points around that route?

texts:
I have a route between 2 points by means of OSM routing service.

I want to get the distance between route and points around. How can I achieve this?

I created this route by the following example https://jsfiddle.net/jonataswalker/079xha47/:

var points = [],
    msg_el = document.getElementById('msg'),
    url_osrm_nearest = '//router.project-osrm.org/nearest/v1/driving/',
    url_osrm_route = '//router.project-osrm.org/route/v1/driving/',
    icon_url = '//cdn.rawgit.com/openlayers/ol3/master/examples/data/icon.png',
    vectorSource = new ol.source.Vector(),
    vectorLayer = new ol.layer.Vector({
      source: vectorSource
    }),
    styles = {
      route: new ol.style.Style({
        stroke: new ol.style.Stroke({
          width: 6, color: [40, 40, 40, 0.8]
        })
      }),
      icon: new ol.style.Style({
        image: new ol.style.Icon({
          anchor: [0.5, 1],
          src: icon_url
        })
      })
    };

console.clear();

var map = new ol.Map({
  target: 'map',
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
    vectorLayer
  ],
  view: new ol.View({
    center: [-5685003, -3504484],
    zoom: 11
  })
});

map.on('click', function(evt){
  utils.getNearest(evt.coordinate).then(function(coord_street){
    var last_point = points[points.length - 1];
    var points_length = points.push(coord_street);

    utils.createFeature(coord_street);

    if (points_length < 2) {
      msg_el.innerHTML = 'Click to add another point';
      return;
    }

    //get the route
    var point1 = last_point.join();
    var point2 = coord_street.join();

    fetch(url_osrm_route + point1 + ';' + point2).then(function(r) {
      return r.json();
    }).then(function(json) {
      if(json.code !== 'Ok') {
        msg_el.innerHTML = 'No route found.';
        return;
      }
      msg_el.innerHTML = 'Route added';
      //points.length = 0;
      utils.createRoute(json.routes[0].geometry);
    });
  });
});

var utils = {
  getNearest: function(coord){
    var coord4326 = utils.to4326(coord);
    return new Promise(function(resolve, reject) {
      //make sure the coord is on street
      fetch(url_osrm_nearest + coord4326.join()).then(function(response) {
        // Convert to JSON
        return response.json();
      }).then(function(json) {
        if (json.code === 'Ok') resolve(json.waypoints[0].location);
        else reject();
      });
    });
  },
  createFeature: function(coord) {
    var feature = new ol.Feature({
      type: 'place',
      geometry: new ol.geom.Point(ol.proj.fromLonLat(coord))
    });
    feature.setStyle(styles.icon);
    vectorSource.addFeature(feature);
  },
  createRoute: function(polyline) {
    // route is ol.geom.LineString
    var route = new ol.format.Polyline({
      factor: 1e5
    }).readGeometry(polyline, {
      dataProjection: 'EPSG:4326',
      featureProjection: 'EPSG:3857'
    });
    var feature = new ol.Feature({
      type: 'route',
      geometry: route
    });
    feature.setStyle(styles.route);
    vectorSource.addFeature(feature);
  },
  to4326: function(coord) {
    return ol.proj.transform([
      parseFloat(coord[0]), parseFloat(coord[1])
    ], 'EPSG:3857', 'EPSG:4326');
  }
};


texts: Adding a table and styling it in arcgis online pop up

texts:
I am researching adding a table and styling the table in a popup in arcgis online.  Could someone provide me a html, css, and javascript examples?  I am fairly new at these languages, but I am hoping to have a relatively good customized pop up with a table and styling before the end of the week.

The table in the pop up will have two columns.  The first column will have the feature name and the second column will have the value associated with that feature class.  The text will need to be in Arial font.  The table itself will need to have alternating color rows: white and light blue.  Finally, the text alignment for the first column needs to be aligned to the right and the second column text aligned to the left.  Eventually this pop up will include hyperlinks and images to reports collected in the field using ArcGIS Worforce, ArcGIS Survey 123, and ArcGIS collector.

Here is some sample code I wrote and applied to the customized pop up.

ObjectID:
{OBJECTID}

FIPS Code
{FIPS}

County:
{County}

State:
{STATE}

Total Population:
{Total_Population}

Population Under 15:
{Population_Under_15}

Unfortunately the CSS rule I tried to apply to the left column did not work.  Any suggestions?


texts: Splitting a GeoTIFF raster file into tiles in meters

texts:
I am new to GIS in general and I have started working with GDAL and Rasterio. What I am looking for is a method in GDAL or Rasterio that mimics the ArcGIS's SplitRaster function when the parameter tile_size is in meters (i.e. "unit" flag set to metres). I have seen multiple similar questions asked but I can only see solutions that would split/tile a large GeoTIFF by pixels.
In summary,
My input: A large GeoTIFF file (Planetscope)
Expected output: Bunch of GeoTIFF files (preserving geo-references) that are tiled of size 100m x 100m.
I have tried gdal_translate and techniques described in this question but I need tiling dimensions to be in meters not pixels.
Edit (Added full gdalinfo of the raster below):
Driver: GTiff/GeoTIFF
Files: Chicago.tif
Size is 10134, 9211
Coordinate System is:
PROJCS["NAD83 / UTM zone 16N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","26916"]]
Origin = (433874.999999512860086,4646969.999883539974689)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  433875.000, 4646970.000) ( 87d47'53.12"W, 41d58'19.02"N)
Lower Left  (  433875.000, 4619337.000) ( 87d47'42.00"W, 41d43'23.11"N)
Upper Right (  464277.000, 4646970.000) ( 87d25'52.21"W, 41d58'26.10"N)
Lower Right (  464277.000, 4619337.000) ( 87d25'46.20"W, 41d43'30.13"N)
Center      (  449076.000, 4633153.500) ( 87d36'48.39"W, 41d50'55.12"N)
Band 1 Block=10134x1 Type=Float32, ColorInterp=Gray


texts: Add WMTS (GRB) to Arcgis Silverlight application

texts:
I want to simply add the WMTS of the GRB (agiv belgium) to my ArcGIS silverlight application (in xaml). I have found and succesfully tested the sample code on the arcgis developers website, but when I use the GRB WMTS is gives an error "Tiling scheme not supported at this time". Does this have something to do with the CRS?

The GRB WMTS = http://grb.agiv.be/geodiensten/raadpleegdiensten/geocache/wmts?service=WMTS&request=Getcapabilities

    <esri:Map WrapAround="True" Name="MyMap"  Margin="0,0,-592,-347">
        <esri:WmtsLayer ServiceMode="KVP"
              Url="http://grb.agiv.be/geodiensten/raadpleegdiensten/geocache/wmts?"/>
    </esri:Map>


texts: How to insert georeferenced PNG into correct location in AutoCAD Map3D?

texts: I have a PNG of building plans that I will need to digitize in AutoCAD Map3D. I have georeferenced the PNG in ArcGIS, which has produced a world file (.pgwx) in the same directory in which the image is stored. When I attempt to bring this image into AutoCAD using either Insert Image or Data Connect, I am prompted to specify an insertion point. It is clearly not recognizing the resource file. How can I open a georeferenced PNG in AutoCAD Map3D in the correct location?