texts: ArcPy - List layers and list fields in one code

texts:
I am trying to combine ListDataframes, ListLayers and ListFields in one ArcPy code.

I can do it one by one (first list layers in one code, and then list fields in second code), but I am trying to do in one code. So, basically I want to list all data frames, layers, all layers with their projection and get field info on those layers (field name, field type). Code by code is not a problem, but combining two codes is.

So far, I have a piece of code, but I do not get results I want. I do get field name and type, but no layer names and dataframes. I know I also have to include Describe tool for spatial reference, but not sure where. This is my code so far. Works, but not as expected.

import arcpy
mxd = arcpy.mapping.MapDocument(r"CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
lyrList = arcpy.mapping.ListLayers(mxd, "", df)
for lyr in lyrList:
     if lyr.supports("dataSource"):
         for field in arcpy.ListFields(lyr.dataSource):
             print (field.type, field.name, field.aliasName)
for lyr in arcpy.mapping.ListLayers(mxd, df):
    name = lyr.name
    print name
    del mxd

As I said, this is a result of piece of my code, And I get this:

00_BASE
SUBGROUP
shore_polygon
shore_polyline
sea

On this list, I would like to look like this:

00_BASE (This is layer name)
SUBGROUP (This is subgroup name)
shore_polygon (FID, Shape, ID) (ObjectID, Shape, Long)
shore_polyline (FID, Shape, ID) (ObjectID, Shape, Short)
sea (FID, Shape, Name) (ObjectID, Shape, Text)

EDIT 1:
This code works for printing layer name, field name, and field type. Still can not get to print group and subgroup layer names.

import arcpy
mxd = arcpy.mapping.MapDocument(r"CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "Main Map")
lyrList = arcpy.mapping.ListLayers(mxd, "", df)
for lyr in lyrList:
    if lyr.supports("dataSource"):
        for field in arcpy.ListFields(lyr.dataSource):
            print (lyr.name, field.name, field.type)
del mxd

And result looks like this:

 (u'shore_polygon', u'FID', u'OID')
 (u'shore_polygon', u'Shape', u'Geometry')
 (u'shore_polygon', u'Id', u'Integer')
 (u'shore_polyline', u'FID', u'OID')

Still want group and subgroup layer names, to look like this:

 GROUP LAYER NAME
 SUBGROUP NAME
 (u'shore_polygon', u'FID', u'OID')
 (u'shore_polygon', u'Shape', u'Geometry')
 (u'shore_polygon', u'Id', u'Integer')
 (u'shore_polyline', u'FID', u'OID')


texts:
map.getView().calculateExtent(map.getSize())

I am able to get extent using the above but I want the extent on every zoom/pan of map.. do we have any functions that can do it? i looked at zoom interactions and zoom to extent but none of those do what i want.


texts: Print Openlayers WMS map with fixed scaling

texts:
I need to print a map with a fixed scale (1:500 and 1:1000) on A4. The online map is done with Openlayers, the print request is sent to QGIS Server using WMS getPrint function.

I mange to get an approximate scale of around 1:500 (1:594 horizontally and 1:554 vertically), but I can't figure out how to fix my formulas to get an exact solution. The formula I use to calculate the extent in meters is:

extent_width = width*printResolution/screenScale
extent_height = height*printResolution/screenScale

To calculate printResolution in m/px I use an approach described in this post:

  function printResolution(scale) {
    var unit = map.getView().getProjection().getUnits();
    var inchesPerMetre = 39.3700787;
    var dpi = 120;

    return scale / (ol.proj.Units.METERS_PER_UNIT[unit] * inchesPerMetre * dpi);
  }

To calculate the actual screenScale I use the approach described in this post:

  function screenScale() {
    var unit = map.getView().getProjection().getUnits();
    var resolution = map.getView().getResolution();
    var inchesPerMetre = 39.3700787;
    var dpi = 96;

    return resolution * ol.proj.Units.METERS_PER_UNIT[unit] * inchesPerMetre * dpi;
  }

And finally I put it together calculation the size of the print box on screen using:

  var w = Number(size[0])*printResolution(scale)/screenScale()*18000;
  var h = Number(size[1])*printResolution(scale)/screenScale()*18000;

  var bounds = map.getView().calculateExtent([w,h]);
  var printBox = [
    [bounds[0], bounds[1]],
    [bounds[0], bounds[3]],
    [bounds[2], bounds[3]],
    [bounds[2], bounds[1]],
    [bounds[0], bounds[1]]
  ];

So I'm aware of errors in the calculation (as I have to multiply by 18000 to get an approximate scaling), but I can't figure out what is wrong.

Here is the complete code in a jsfiddle showing the described behavior.

PS: This is a follow up to this question posted some time ago.


texts: Wrong projection reference in coordinate system Modifying original coordinate system in ArcGIS Desktop?

texts:
Working in Germany, where Gauss_Kruger is very often used (as well as ETRS1989_zone 32N), and trying to find out about WKID when I realised that on my machine, the projection reference of the DHDN_3_Degree_Gauss_Zone_3 is set on Transverse Mercator, whereas for Gauss_Kruger 1,2,4 or 5, it's set on Gauss_Kruger, which makes sense.

I have compared with the Parameters on other machines, where they are correct (Gauss_Kruger). Somehow, I have changed them and I need to modify, and not by setting a new custom Coordinate System but by modifying the file within ArcGIS or even, deleting it. Now I read things there is no Coordinate System Folder but " Instead, a virtual directory structure is built directly from the coordinate system data."    (http://support.esri.com/en/knowledgebase/techarticles/detail/39488)

How can I do this?


texts: How to eliminate lines that are added when "join multiple line" plugin is used?

texts: I'm new in QGIS. I have three road shapes. If I merge them, I cannot use "network analysis" (shortest path) because the merged lines are not continuous. I tried to use "Join Multiple Lines" and, actually, it worked. However, when I used it, some lines were added and I don't how to eliminate them.


texts:
I have 15 years' rasters along time serials. Each year contain 46 rasters.The file such as N2000001.tif, N2000009.tif,......, N2000361.tif

I want to stack same year rasters as single file.Such as N2000.tif, which contians 46 bands.

Run code as follows:

import numpy as np
from osgeo import gdal
import os

inpath=''
outpath=''

driver=gdal.GetDriverByName('GTiff')
driver.Register()
file0=
data0=gdal.Open(file0,gdal.GA_ReadOnly)  # base file for driver.CreatCopy()

 # *get the all tif rasters filenames*
filenames = [filename for filename in sorted(os.listdir(inpath))
            if os.path.splitext(filename)[1]=='.tif']
# *sort all rasters name*
single=sorted(set(filenames),key=filenames.index)

for year in range(2000,2016):
     i=1
      #buid each year stack file
     outfile=os.path.join(outpath,'N'+str(year)+'.tif')
     outds=driver.CreateCopy(outfile,data0)
     for name in single:
          if str(year) in name:        # *run each year*
              file=os.path.join(inpath,name)
              print name+'is processing'

              data=gdal.Open(file,gdal.GA_ReadOnly)
              band=data.GetRasterBand(1)
              arr=band.ReadAsArray()
              outds.GetRasterBand(i).WriteArray(arr)

              i+=1                 # *the band sequence number*
      outds=None

system hint error:

  File "<ipython-input-5-a57257d55372>", line 35, in <module>
    outds.GetRasterBand(i).WriteArray(arr)

AttributeError: 'NoneType' object has no attribute 'WriteArray'


texts:
I have a point layer of students, and a polygon layer of district boundaries.

How can I select those students who are not in the district boundaries, a Select By "Not-in" Location, so to speak.


texts: Fastest way to join many points to many polygons in python

texts:
I need to attribute a geographical area (based on shapefile polygons) to 8 million GPS points. There are 405 polygons. I currently have the code below, and a quick extrapolation tells me it will take 84 hours to compute this.

I am new to GIS, but intuitively, I feel that there has to be a smarter way around this than, for each point, randomly testing each polygon with a point-in-polygon function until there is  match. For example, simply putting the points and polygons in 2 groups (such as "East" and "West") simply based on their easternmost/westernmost longitude may divide by two the number of operations to perform by almost 2.

Are there indeed such much efficient algorithms? Are any of them easily implementable in python?

from shapely.geometry import shape, Point, Polygon
import pandas as pd
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

plt.figure(figsize=(16,16))
m = Basemap(projection='merc',llcrnrlat=15,urcrnrlat=55,llcrnrlon=70,urcrnrlon=140,lat_ts=20,resolution='c')
m.readshapefile(china_shapefile, 'prefectures', drawbounds=True)

print len(m.prefectures)

#get data
china_shapefile="C:/Users/.../CHN_adm2_SIMPLIFIED/CHN_adm2"
points_csv="C:/Users/.../sample_map_data.csv"

points=pd.read_csv(points_csv)
points=points.dropna()

lat=points["Latitude"].tolist()
lon=points["Longitude"].tolist()
coordinates=zip(lat,lon)

for coordinate in coordinates:
    for info,shape in zip(m.prefectures_info,m.prefectures):
        poly = Polygon(shape)
        #print info["NAME_1"],info["NAME_2"]
        xpt,ypt = m(coordinate[1],coordinate[0])
        point1 = Point(xpt,ypt)
        if point1.within(poly):
            print coordinate,info["NAME_1"],info["NAME_2"]#,info["NAME_3"]
            break


texts: Truncate Table tool with geoprocessing service?

texts:
A model I built with ModelBuilder works as expected in Desktop yet generates numerous error messages when published as a geoprocessing (GP) service in ArcGIS Server 10.1 SP1. I've decided to start simple then and am trying to get just a single tool to work.

Here is my model. It has one parameter, a table or a feature class to truncate. When running in ArcGIS Desktop 10.1 SP1, user can specify what object should be truncated. When publishing as a model after getting a successful result, I use default settings for all settings (and Info for Message Level).

When running the GP service task in the same ArcGIS Desktop, I get this window in order to supply an object for truncating.

After pointing to a dataset (stored in a 10.1 file geodatabase on a local disk drive), the GP task gets running. The task does not truncate table and produces an error:
000187 : Only supports Geodatabase tables and feature classes.

My guess is that it has something to do with the Input mode of the task (User defined value) which makes it impossible to supply to my task a feature class from a file geodatabase in the file system. It seems as model uses the empty FeatureSet object instead of a feature class I supply when running the GP service task.

Is there a way to make the Truncate Table GP tool to use another feature class (not the same used when running the model first in Desktop before publishing) when published as a GP service?

The model for truncating a table works fine when there are no parameters exposed and when running the GP task of a published service, it runs over the same feature class I've used when running the model first to generate a result.

So, the GP service does work with feature classes yet only with the same one used when running the model before publishing.

Is it possible to supply a feature class to a GP task that differs from the feature class used to run the model and generate the result?

NB: the feature class is not copied to the arcgisserver folder, it is being processed properly right in the geodatabase on the disk, so this part works properly.


texts: What is the Formula to get an Image's Bounding Box from it's GDAL GeoTransform Array?

texts:
I have a process which needs to read in an Image's GeoTransform Array, and convert it to the Bounding Box

i.e. If the Input is:

<GeoTransform>
   7.8087700868360430e+05,  9.9974951285737870e+00,  0.0000000000000000e+00,
   2.0276558410121971e+06,  0.0000000000000000e+00, -1.0000024189695811e+01
</GeoTransform>

I need the Output as

Upper Left:780877.009, 2027655.841
Lower Left:780877.009, 2010425.799
Upper Right:799692.295, 2027655.841
Lower Right:799692.295, 2010425.799

I see that you can do this using GDAL's Python bindings, but I need to do this without any dependencies.

Is there a formula that could be applied to get the expected Bounding Box?


texts:
QGIS is not starting because of deleted file libgeos_c.so.1,How can I fix it
It was originally running but some how I have deleted this file for some porpose installation of postgis for postgress but can't start it now.
It shows error as follows on ubuntu 14.04 machine

  /usr/bin/qgis.bin: error while loading shared libraries: libgeos_c.so.1: cannot open shared object file: No such file or directory


texts: Running the ogr2ogr utilities in a batch file. Problems with command project

texts:
I have made a  batch file containing a list of two commands. The first (-clipsrc) works properly, but the second not (-s_srs EPSG:4326 -t_srs EPSG:3857).

Could anyone tell me where is the fail?.

This is the batch file:

cd /d c:\data\PhiladelphiaBaseLayers
set ogr2ogrPath="c:\program files\QGIS Dufour\bin\ogr2ogr.exe"
for %%X in (.shp) do %ogr2ogrPath% -skipfailures -clipsrc c:\data\PhiladelphiaBaseLayers\clipFeature\city_limits.shp c:\data\PhiladelphiaBaseLayers\clipped\%%X c:\data\PhiladelphiaBaseLayers\%%X
for %%X in (.shp) do %ogr2ogrPath% -skipfailures -s_srs EPSG:4326 -t_srs EPSG:3857 c:\data\PhiladelphiaBaseLayers\clippedAndProjected\%%X c:\data\PhiladelphiaBaseLayers\clipped\%%X


texts: I would like to get a list of the values of a GeoTIFF. Is that possible using 'gdalinfo' command? If not, is there any other software that can show that?


texts: QGIS: Is it possible to merge two stacked raster layers keeping color rendering and transparency?

texts:
Extremely new to QGIS and was hoping to merge 2 layers and maintain an "overlay" keeping layer transparency and color rendering?

I seem to loose the top layers layer properties when merging and it reverts to it's original look.

Before "Merge":

After "Merge":

Using QGIS 3.2.2


texts: Aggregate table with dates and average geographic coordinates


texts: I am trying to make a 3D representation of the map below. I have 3 polygons which represent contour lines. It starts at 0ft for the outside polygon, 10ft in the middle and 20ft for the very center polygon. The end result would look something like a weird shaped pancake that is 20ft thick within the middle polygon and shrinks down to 0ft at the very edges. As an end result I would like to be able to compute the volume of the end shape. I am using ArcGIS10.3.1. I do not think I have the 3D Analyst extension so that may be a problem.


texts: What is a proxy entity (AcDb Zombie Entity) and how are they introduced into a file?

texts:
I am doing some transformation and translation of a few DWGs using FME desktop and I received this error,

  AutoCAD proxy entity type 'AcDbZombieEntity' not supported.

Some preliminary research revealed that 'AcDb Zombie Entity' is another term for 'proxy entity' (link) and explains that the Export to AutoCAD command should be used to rid your data of the zombie entities (link). This helps me proceed with my task but does not help me understand the AutoCAD object data model any better.

What exactly is a proxy/zombie entity? How are proxy entities introduced into a file? The creator of the DWG does not know what a proxy entity is or how they inadvertently introduced them into the file.


texts: For Loop - Map Algebra (ArcGIS) using Python

texts:
ORIGINAL QUESTION

Im relatively new to python and these fourms but Im am having trouble utilizing for-loops for spatial analyses within ArcGIS. It does not appear to recognize the variable, RasterImage!

Here is a part of the code im having trouble with...

import arcpy, os, sys, string
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")

arcpy.env.workspace = r'O:\DATA\test'
InputFolder = arcpy.env.workspace

Temp4 = InputFolder + '\\' + 'temp4'

RasterList = arcpy.ListRasters()
for RasterImage in RasterList:
    print (RasterImage)
    Temp4 = FocalStatistics('RasterImage', NbrRectangle (9, 9, 'CELL'), "RANGE", "DATA")
    print ('Focal Statistics Done')

I receive the following error

"Traceback (most recent call last):
File "O:\test.py", line 67, in <module>
Temp4 = FocalStatistics('RasterImage', NbrRectangle (9, 9, 'CELL'), "RANGE", "DATA")
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 4796,   in FocalStatistics
ignore_nodata)
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in   swapper
result = wrapper(*args, **kwargs)
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 4790,    in wrapper
ignore_nodata)
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_base.py", line 474, in <lambda>
return lambda *args: val(*gp_fixargs(args))
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000865: Input raster: RasterImage does not exist.
Failed to execute (FocalStatistics)."

It does not seem to recognize my for in loop of the RasterImage - any ideas? The print statment after the loop shows that it is recognizing the file (in python) but that the spatial analyses used in arcpy are not.

I thought it might have to do something with the workspace used so I thought that adding

RasterImage = InputFolder + '\\' + RasterImage

might help locate the full path of the raster datasets but to no avail! Im lost for ideas and the only help files that utilize spatial analyses on the ArcGIS resource page do not utilize loops in their scripts. ¨

UPDATE

Issue has been resolved by implementing the following:

The original question displaying the error 'RasterImage' does not exist was cause in part due to the '' marks and the spaces around the NbrRectangle. I have implemented part of my new code using the suggestions by nmpeterson whereby the geoprocessing is called through the function

def focalStats(raster):
   Temp4 = FocalStatistics(raster, NbrRectangle(9, 9, 'CELL'), "RANGE", "DATA")
   return Temp4

While methodologies in which several temporary files rely on one another (i.e. a geoprocess within another geoprocess) i've used the .save function into a temporary directory as shown by om_henners before deleting all those contents afterwards.


texts: I am creating a map in QGIS and I am using a hex bin layer. To create this layer I used the MMQGIS plugin and created a grid layer. However, the size of my hexagons that I created are too small, and I was wondering how I can change their size to make them bigger.


texts: Select the bands which not start with 'B' using regex in GEE

texts:
The image has the band : ['B1', 'B2', ..., 'B11', 'NDVI', 'NDWI'].
The goal is to select the bands which not start with 'B' using regex. In another words, ['NDVI', NDWI'].
What I know is that using image.select("B.*") to select bands start with 'B'. But I do not know to reverse the selection.


texts:
I have feature class (Name: V1) in personal geodatabase with columns:

Height [float] and

Id [int] representing group number

I need to do SQL query that will select eight highest 'Height' values in every group.
I know that syntax below is not working but this is what I want to accomplish:

Reason I choose personal geodatabase instead of file one is that I can't create query below (I don't know why, because in personal geodatabase I can):


texts: I am exporting my map with qgis2web and I would like to export the opacity of my layers as well. It does not work, somebody has an idea?


texts:
I am trying to get an overlay (intersection) of two collections of polygons
in geopandas:

splits = gp.overlay(nhds, cts, how='intersection')

but result is very weired (see the picture below).

I opened an issue on geopandas github, but wondered if someone saw similar behavior or can detect what is going on (I am inclined to blame spatial index) here from the image


texts: How to create memory layer and populate fields?

texts:
Seems like a basic operation using PyQGIS but can't see what I am missing. I want to create a memory layer with the exact same features and attributes from a shapefile. I have looked several posts such as:

Copying attributes from point layer to memory linestring layer PYQGIS
pyQGIS manually copy all features with attributes from a layer into a memory layer

The following script creates a memory layer, gets the correct fields and correct number of features but it doesn't populate the fields with data:

input = "C:/Users/Me/Desktop//example.shp"
layer = QgsVectorLayer(example,"line","ogr")

temp = QgsVectorLayer("LineString?crs=epsg:4326", "result", "memory")
temp_data = temp.dataProvider()
temp.startEditing()

layer_fields = layer.dataProvider().fields().toList()
attr = layer.dataProvider().fields().toList()
temp_data.addAttributes(attr)
temp.updateFields()

feat = QgsFeature()
for elem in layer.getFeatures():
    feat.setGeometry(elem.geometry())
    feat.setAttributes(attr)
    temp.addFeatures([feat])
    temp.updateExtents()

temp.commitChanges()
QgsMapLayerRegistry.instance().addMapLayer(temp)


texts: Calculate number of pixels and more variables by classes in QGIS

texts:
I have one band raster that has 3 possible values: 1,2 or 3.
I want to calculate for this raster different variables, for example how many pixels we have for each values and also to do some band math for them.
Another challenge is that the raster I have is one but it contains data about different polygons.

My end goal is to have kind of tabular data based on the different polygons that will contain data about how many pixels I have in each polygon, how many pixels are in each class, etc.


texts: Self split geopandas Linestring geodataframe in a fast way without loosing all attributes

texts:
I have geodataframe of many LineStrings. The LineStrings intersect but are not split at those intersections. My current solution to achieve this is from here:
network = gp.read_file(filenameNetwork)
newNetwork = gp.GeoDataFrame()
for splittedGeom in network.geometry.unary_union:
    part = gp.GeoDataFrame([[splittedGeom]], columns=['geometry'])
    newNetwork = newNetwork.append(part)

However I loose all columns with this solution. I also tried this, but it takes way too long:
from shapely import ops

streets = streets.reset_index(drop=True)
streets = streets[['geometry', 'costs']]
headers = list(streets.columns)

index = 0
newStreets = gp.GeoDataFrame( columns=['geometry'])
for line in range(len(streets)-1):
    print(line, len(streets))
    linegeom = streets.at[line, 'geometry']
    isNotSplitted = True
    for line2 in range(len(streets)):
        if line2 == line:
            continue
        linegeom2 = streets.at[line2, 'geometry']
        if linegeom2.crosses(linegeom):
            try:
                linegeomsplitted = ops.split(linegeom, linegeom2)
            except:
                continue
            isNotSplitted = False
            for split in range((len(list(linegeomsplitted.geoms)))):
                splittedline = (list(linegeomsplitted.geoms))[split]
                for head in headers:
                    if head == 'geometry':
                        headValue = splittedline
                    else:
                        headValue = streets.at[line, head]
                    newStreets.at[index, head] = headValue
                index += 1
    if isNotSplitted:
        for head in headers:
            headValue = streets.at[line, head]
            newStreets.at[index, head] = headValue
        index += 1
streets = newStreets
streets = streets.drop_duplicates(subset=['geometry'],
                                  keep='first')

Any suggestions?


texts: Enhance performance of a join query for mapnik using postgis

texts:
I use mapnik to render map with datasources from a postgis database. One query is a join on two table, one with a geometry (a polygon) and the other one with some fields I want to use for rendering. Unfortunately, this query is very slow...

Here is the slow query:

SELECT table1.id, table1.geometry, table2.fk_table1, table2.field FROM table1
INNER JOIN table2 ON (table1.id = table2.fk_table1)
WHERE table1.geometry && SetSRID('BOX3D(x1 y1, x2 y2)'::box3d,900913);

Here is the query plan I got from an EXPLAIN:

Hash Join  (cost=118896.94..262210.34 rows=27942 width=1099)
    Hash Cond: (table1.id = table2.fk_table1)
    ->  Seq Scan on table2  (cost=0.00..73740.89 rows=4786689 width=16)
    ->  Hash  (cost=113876.63..113876.63 rows=33865 width=1083)
        ->  Bitmap Heap Scan on table1  (cost=1112.15..113876.63 rows=33865 width=1083)

This query can be faster if the bbox is first applied on table1.geometry. However, mapnik does not help me there. Indeed, mapnik builds a query from a table/subquery by adding a WHERE clause on it

  SELECT AsBinary("geometry") AS geom
  FROM subquery
  WHERE "geometry" && SetSRID('BOX3D(x1 y1,x2 y2)'::box3d,900913)

With the following subquery:

( SELECT table1.id, table1.geometry, table2.fk_table1, table2.field FROM table1
INNER JOIN table2 ON (table1.id = table2.fk_table1) ) as subquery

Basically, I need to enhance the speed of this query knowing that mapnik will add a where clause on the subquery. Thx for your help.


texts: Using ModelBuilder to Iterate, Join & Export Feature Classes?

texts:
I am trying to make what I thought was a very simple batch processing tool, but I can't get the output to work properly. I would like the model to iterate through a geodatabase, join each feature class with the counties feature class, export that join (Copy Features) to a new feature class, remove the join, and repeat with the next for a total of 240 iterations.

After looking at the geoprocessing results, I believe the issue lies somewhere in the Copy Features tool. The model results in one feature class with the name of the first iteration's output, "AGWT", and the values of the last iteration, "YTVI", in the table. So it appears to run through each of the 240 feature classes, but keeps the name of the first output and thus overwrites all subsequent outputs.

The Copy Features Output looks like this: C:\ ...\Exports.gdb\%Name%c

Any ideas/questions/recommendations?


texts: Add MapImage to MapImageLayer in ArcGIS API for Javascript 4?

texts:
I am looking for a way to add a georeferenced image file to a map in the ESRI 4.x Javascript api.  In the 3.x version this was possible through the addImage method on MapImageLayer.  It looked like:

mapImageLayer.addImage(new MapImage({
    extent : extentForDataUrl,
     href : someDataUrlString
}));

In reading the docs for the 4.2 MapImage they state “The MapImage class can also be used to create a geo-referenced image file for use with the MapImageLayer” but I cannot find how this is done.  The addImage method on MapImageLayer appears to be gone and its constructor requires a portal item id or a service url.

The images that I want to add to the map reflect radar/meteorology data that is updated every couple of minutes.


texts: Geographically Weighted Regression - Problem to run model with floating explanatory variables

texts:
Can anybody tell me please why I cannot run the GWR model in ArcGIS 10.3 with a floating explanatory variable?

When I convert the floating explanatory variable in an integer variable, then, the model runs.But, when I use floating values, I got following error message - 040038 : Results cannot be computed because of severe model design problems. And at the description of the error in help tool box is - there is a multicollinearity problem. But, I used only one explanatory variable. So, multicollinearity or redundancy problem cannot exist.


texts: How force bing logo to bottom quickmap service

texts:
I'm trying to compose a map on QGIS 2.14.6 using Bing Satellite basemap with quickmap service plugin.

However, contrary to other basemaps like OSM Mapnik, I am not able to force showing bing logo to bottom of a map.

Going through Bing Satellite layer's properties I have switched on the box "place the credit on bottom right corner".

Does anyone knows how to force bing credit to bottom quickmap service, so I don't have legal problems using this service?


texts: Shapefile to ArcGIS JavaScript JSON format?

texts:
Possible Duplicate:
  How to convert SHP into GPX or JSON?

We're in the proces of adding a new feature to our Web Application where we employ the use of ESRI's ArcGIS Javascript API alongside other ESRI server products.

How would one go about converting an ESRI Shapefile format into ArcGIS Javascript API compatible JSON?

Near as I can tell, there are many examples of doing this very "thing" but converting to GeoJSON. AFAIK ArcGIS Javascript API does not understand or work with GeoJSON.


texts: Check if ArcGIS Server 10.0 Service is running and start if not

texts:
I am using AGSSOM.exe to start and stop ArcGIS Server 10.0 MapServices. I would like to slightly adapt this script to first check if the services are running, and then to start them if they have stopped.

Does anyone know if this is possible and how?


texts:
I have create a grid in QGIS. But I don't have coordinates, I would like to sketch on the basis of Google satelitte, that uses Pseudo Mercator projection. I tried GRASS mkgrid, but something went wrong.

  ERROR: North must be larger than South, and other warnings.

I read many forums about Pseudo Mercator mistake, but I cannot find the solution. What are my options for useable method?

I tried this settings for fishnet

The create grid not solving for me, because necessary define number of rows and columns, therefore I need to use mkgrid


texts: Confusion about "alpha" band for MODIS Aqua

texts: I'm working with modis:09:CREFL data for a harmful algal blooms paper, and I am trying to figure out how to call the chlorophyll-a band, which I believe is the "alpha" band, but I'm not sure because it doesn't have a clear description. Is it possible I have to derive it myself? I doubt it because NASA's SeaDAS software has chl-a as a built-in option.


texts: How to generate arcs from attribute data?

texts:
What I am trying to achieve is an ellipse-like shape from a digitised point location and cardinal measurements.

For example, point location at 531500,181500 with the attributes North - 4m, South - 6m, East - 5m, West - 5m.

I am guessing the best approach is to calculate the start and end point of 4 arcs in postgis but i'm not sure where to start.


texts: Map library supporting linestring curve rendering to produce flow map

texts:
Is there some function/plugin for any javascript map library (OL, leaflet, ArcGIS API for JS) to render WFS or Geojson LineStrings not as simple line but as curves?

I want to make mobility flow maps of region sequences (centroid of areas) and drawing everything just as LineStrings is just messy and confusing.

I would like to be able to create not just overlapping LineStrings, but curves and possibly have styling options like adding an arrow representing a direction of flow (middle, per interval, etc.).

Another possibility is to draw "subway" like map/diagram where LineStrings are shifted next to each other.

Also, there could be some function to automatically sum the "edge" values between same nodes and represent that part only once with weighted LineString.


texts:
I have a dataset with lat/longs mostly in the UK and want to add a column indicating parliamentary constituency. The whole dataset is ~50k rows by 30 columns. A short sample is:

    lat      long Age city_id

1  50.78839 -3.648326  26    1721
2  51.47254 -2.500891  19    1724
3  52.76901 -1.202122  27    1727
4  53.57782 -1.002811  27    1730
5  54.66645 -5.913355  45    1733
6  51.06004 -1.309250  27    1736
7  51.50796 -2.606150  23    1724
8  51.47989 -1.442205  44    1739
9  53.72372 -1.862787  55    1742
10 50.72398 -1.788323  31    1745

I've downloaded shapefiles for Westminster parliamentary constituencies from the ONS: https://data.gov.uk/dataset/westminster-parliamentary-constituencies.

I've then used over in the sp package in R to check which points fall within which polygons. The code looks like:

WMCs<-readOGR(file.choose())
data_SPDF<-SpatialPointsDataFrame(data[,c("long","lat")],data,CRS("+proj=longlat +datum=WGS84"))
data_SPDF<-spTransform(data_SPDF,CRS(proj4string(WMCs)))
data_over_WMCs<-sp::over(data_SPDF,WMCs)
data$WMC<-data_over_WMCs$NAME

This works for most of the data, but around 8k points just get NA. For those first ten lines:

        lat      long Age city_id              WMC
1  50.78839 -3.648326  26    1721             <NA>
2  51.47254 -2.500891  19    1724        Kingswood
3  52.76901 -1.202122  27    1727     Loughborough
4  53.57782 -1.002811  27    1730             <NA>
5  54.66645 -5.913355  45    1733             <NA>
6  51.06004 -1.309250  27    1736             <NA>
7  51.50796 -2.606150  23    1724             <NA>
8  51.47989 -1.442205  44    1739          Newbury
9  53.72372 -1.862787  55    1742          Halifax
10 50.72398 -1.788323  31    1745 Bournemouth East

These points do fall within the UK and within parliamentary constituencies. The ones that fail don't generally fall at the boundaries between constituencies - they're mostly inside, and seem randomly interspersed with points that succeed. Identical code with different shapefiles has allowed me to successfully add a column with MLSOAs.

Any ideas what could be going wrong?


texts: Displaying 2 rasters side by side in ArcGIS Pro Notebook

texts:
I have drone imagery that I am performing various processes through a notebook. One thing I would like to do is display the RGB imagery next to a single band plant height raster for visual comparison using Render, but it is only displaying the second raster.
Is it possible to do this?
Here is my code:
name = "PLAZ"
intgdb = "G:\\...\\IntermediateData.gdb\\"
height = Raster("G:\\...\\name + "_vegHeight.tif")
heightS = arcpy.management.Shift(height, intgdb + name + "_shift", 775, 0)
arcpy.ia.Render(ortho, {'bands': [3, 2, 1]})
arcpy.ia.Render(heightS, {'bands': [1]})

I was told that because the two rasters occupied the same XY space I had to shift one of them, so I tried that in line 4 but no difference. It seems to only be rendering the last command.
Is there a way to put them side by side to take advantage of the big empty space?
I've also included a screenshot.

And here is the result of the Shift, so that tool worked (and just to make sure my question is clear, this is basically what I would like to display in my notebook).


texts: How to move vectors under geom.Polygon.circular layer

texts:
The vectors are hidden by the geom.Polygon.circular, how do I move the vectors under layer of Circular?

Circular is added like ol.Feature (localization)
Vectors (points)

The Z-index for the vectors element does not change its position for the circular feature.

How do I solve this?

I need solution strictly for OSM OpenLayers and JavaScript


texts:
Am trying to clip multiple rasters but getting a runtime error, what could be the problem with my code

import arcpy, os

arcpy.env.workspace = r"E:\Geos4s\GeoS4S Module Geospatial Support for Arid Ecosystem Management\prec_30s_esri\converted"
arcpy.env.overwriteOutput = True

# Clip region
poly = r"E:\Kenya GIS Resources\Kenya New Districts1\Kenya_CountyID.shp"
out_folder =  r"E:\Geos4s\GeoS4S Module Geospatial Support for Arid Ecosystem Management\prec_30s_esri\clipped"
rasters = arcpy.ListRasters()
for raster in rasters:
    out_raster =(out_folder,raster)
    arcpy.Clip_management(raster,'#',out_raster,poly,'#','ClippingGeometry')
    print 'Clipped: %s' %raster
print arcpy.GetMessages()


texts: ArcMap.exe processes stay open after closing ArcMap?

texts:
I realized that this is happening months ago when I was not able to alter table structure in an instance of ArcMap after another instance was closed.
When, for example, I remove or add some fields to a feature class in ArcMap, save and close the document and open ArcCatalog and attempt to delete that particular Feature Class I get the Delete Failed: Failed to delete selected object(s) error.  This usually happens when an instance of ArcMap that contains layers which I am trying to delete is opened so the workaround is to fire up the Windows Task Manager and kill the ArcMap.exe processes that are for some reason still opened.

Anybody else having this issue?

Now running SP3 and the issue has not been solved.

I am running a very similar setup on a comparable machine and I do not have this issue there.


texts: Does Model-Dependent Realism mean the Map is the Territory after all?

texts:
Stephen Hawking's Pipe.

In his book The Grand Design, Stephen Hawking says:

  There is no picture- or
  theory-independent concept of reality.
  Instead we will adopt a view that we
  will call model-dependent realism: the
  idea that a physical theory or world
  picture is a model (generally of a
  mathematical nature) and a set of
  rules that connect the elements of the
  model to observations. This provides a
  framework with which to interpret
  modern science.

Does this mean I no longer need a disclaimer on my apps saying the Map is Not the Territory?


texts: How configure pop-ups in ArcGIS Online for open there attachments?

texts:
Is it possible to open attachments by dint of pop-ups (as files jpg, doc)? According to this Help – it’s possible:

“Pop-ups can be configured to show attachments as links if the map includes a features layer that has attachments in its data. You will only see an option to show attachments if the layer has them. When the pop-up is displayed, the attachment links are listed at the bottom of the window. Clicking one of the links opens the attachment in a new window”

But they don’t explain how to browse theme with attachment to ArcGIS Online. How I can do it, if shp theme can’t have attachment, but to ArcGIS Online from ArcGIS it’s possible to browse only shp file as a zip folder.?


texts: KML to shapefile will not convert attributes embedded in the description field

texts:
I have a point shapefile on which I used the "Layer to KML" tool on ArcMap. This tool successfully converted the shapefile into a KML containing all the attributes embedded in the "Placemark Description" field in Google Earth which serves my purposes. However when I try to go the opposite direction, namely from KML back to layer or shapefile using "KML to Layer" tool in ArcMap, all the attributes are lost and I cannot see them in the resulting Shapefile.

Is there any way to somehow expose these embedded attributes?

Or maybe export the shapefile such that the attributes are not embedded in the description in the first place?

I have third parties who need to edit my shapefiles in Google Earth and send them back to me to convert back to shapefiles. The back and forth from Arc to GE is quite paramount to my situation.


texts: How to prevent a geoprocessing task from being run twice?

texts:
If I run twice a geoprocessing task (using the REST API), what does happen?

Does the geoprocessing service execute the two tasks simultaneously, ...?

How can I configure the geoprocessing service from executing two tasks at the same time?


texts:
I'm trying to use ArcGISTiledMapServiceLayer to show a cached set of image tiles on our local server. There is good help on the arcgis developer's site but I can't seem to find good examples. None of the samples they list actually call this method.

I'm using esri javascript api version 3.7


texts:
I use GeoServer with Oracle DB. The number of records in the database is fewer than 4000, but a feature request takes 5-6 seconds. How can I speed up this?

P.S. Native JAI is false in my GeoServer

UPDATE:

I have view which joins two tables. int view sdo.geometry is polygon with 12 points.

When I open it from layer preview it takes from 5 to 6 seconds.

UPDATE 2

I removed view and joins and it works great. I had to copy columns from one table to another for reading directly from table all fields. If someone knows solution when using views or joins please share :)


texts: Get coordinates from polygon points on MapBox draw?

texts:
This is the code I am using for polygon to draw and I am using MapBox GL draw I am trying to get longitude and latitude of every point from polygon.

var draw = new MapboxDraw({
  displayControlsDefault: false,
  controls: {
     polygon: true,
     trash: true,
  }
});

map.addControl(draw);

map.on('draw.create', function(e) {

    if (lngLat.lng && lngLat.lat == 'polygon') {
        console.log(draw.getAll());
    }
});


texts:
Preface

  For a while, I was wondering about applying the looping/iterable
  functions in QGIS's Virtual Layer. There were several cases when the
  application of a loop would be beneficial and useful, e.g. Intersect
  multiple layers at once
  QGIS
  or Incorrect result while adding points along geometry in
  QGIS.

There is one polyline layer that represents flows (attributes: "ORIGINE", "DEST", "FLUX", and "DIST_KM"). Its attribute table looks as following.

My desire is to obtain some information based on categories and corresponding rules for each category.

By any means, this query will give what I actually need

SELECT COUNT(*) AS anzahl, SUM(FLUX) AS summe
FROM "flows_workday_GK5"
WHERE "DIST_KM" > 0 AND "DIST_KM" < 1
UNION ALL
SELECT COUNT(*) AS anzahl, SUM(FLUX) AS summe
FROM "flows_workday_GK5"
WHERE "DIST_KM" >= 1 AND "DIST_KM" < 2
UNION ALL
SELECT COUNT(*) AS anzahl, SUM(FLUX) AS summe
FROM "flows_workday_GK5"
WHERE "DIST_KM" >= 2 AND "DIST_KM" < 3
UNION ALL
SELECT COUNT(*) AS anzahl, SUM(FLUX) AS summe
FROM "flows_workday_GK5"
WHERE "DIST_KM" >= 3 AND "DIST_KM" < 4
UNION ALL
SELECT COUNT(*) AS anzahl, SUM(FLUX) AS summe
FROM "flows_workday_GK5"
WHERE "DIST_KM" >= 4

Output

In the above example, I showed only five categories, but what if I will need to work with 10, 100, 1000 categories.

Since Virtual Layer supports SQLite library it might be smart to embed a RECURSIVE expression in the afore query. I was inspired by this article: basic recursive query on sqlite3?

My desire is to achieve a table as following

category | anzahl | summe
   1     |   96   | 56996

These are several efforts were done, but somehow I am getting errors.

I. When I try with this query

WITH RECURSIVE cat(x) AS (
     SELECT 0
     UNION ALL
     SELECT (x + 1) AS category, COUNT(*) AS anzahl, SUM(FLUX) AS summe
     FROM "flows_workday_GK5", cat
     WHERE "DIST_KM" > x AND "DIST_KM" < x + 1
     AND x < 4
     UNION ALL
     SELECT (x + 1) AS category, COUNT(*) AS anzahl, SUM(FLUX) AS summe
     FROM "flows_workday_GK5", cat
     WHERE "DIST_KM" > x
     AND x = 4
)

SELECT *
FROM cat

I end up with this error

II. Even if I simplify the query like

WITH RECURSIVE cat(x) AS (
     SELECT 0
     UNION ALL
     SELECT (x + 1) AS category, COUNT(*) AS anzahl, SUM(FLUX) AS summe
     FROM "flows_workday_GK5", cat
     WHERE "DIST_KM" > x AND "DIST_KM" < x + 1
     AND x < 4
)

SELECT *
FROM cat

In this case, I end up with the following error

So, my question is: How shall I adjust the query to be able to get the result by means of RECURSIVE function?

The efforts that were made by @M Bain are absolutely brilliant and lead to the desired output. As well as I got to know CAST-expression. Nevertheless, there are somehow several pitfalls that may lead to the deficient result.

The last 'category' might be incomplete because it won't include all the remained records, see image below. The idea is to aggregate all the rest values in the last 'category'. In my case, the 5th 'category' has to have a '2105282'-value.

ROUND(DIST_KM + 0.5, 0) owns an inevitable 'rounding effect'. 'category' number at certain point will jump over some essential numbers, see below.

Last but not least, I am truly passionate to learn/understand the RECURSIVE technique.


texts: QGIS network analysis, points on graph cannot be reached

texts:
Using QGIS 2.18.14, I try something around with qgis.networkanalysis module as described here: https://docs.qgis.org/2.18/en/docs/pyqgis_developer_cookbook/network_analysis.html

director = QgsLineVectorLayerDirector(network, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = network.crs()
builder = QgsGraphBuilder(crs)

tiedPoints = director.makeGraph(builder, [start_geom] + target_list)
# target_list is a list of QgsPoints()
# start_point is a QgsPoint()
graph = builder.graph()

Everything seems to work so far, tied point are generated correctly, but when I try to find paths from start_point to the targets by

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idSource, 0)

some of the tied_points cannot be reached in the graph (i.e. tree[id_target] =-1 for certain point).

My first thought was that there is something wrong with my vector layer but it seems topologically correct (no gaps or duplicate consecutive point e.g.), no matter if LineStrings or MultiLineStrings.

Since I thought that each tied point should be reachable in the graph, what can be the reason for the 'non-reachability' of some points on such? Could this be a problem of LineString directions in the network layer?


texts: Anyway to remove points outlying shape in ggplot2?

texts:
I am to map a plot of an area and then plot points over it (given by Easting and Northing). However, I found out that there are certain mistakes in the data entered for Easting and Northing resulting in some outlying points, and ggplot (which I am using) tends to reshape the map in order to fit the outlying points.

Is there any quick way to fix this issue? I do not know of any except to identify and correct the points manually. I am using shapefiles from the Ordnance survey.

Here is the screen-shot. The image does not appear to be skewed. However, upon comparison with other images, the differences become apparent.

**

EDIT

**

The code I'm using is as follows.

countyRegion<- readShapePoly(file.choose())
norfolkCounty<- countyRegion[countyRegion$NAME=="Norfolk County",]
#convert shape file into data that can be plotted on graph
gpclibPermit()
norfolk<- fortify(norfolkCounty,region="NAME")
Norfolk<- merge(norfolk, norfolkCounty@data,by.x="id",by.y="NAME")

#eit data
jul2012<- read.table(file.choose(),header=TRUE,as.is=TRUE,blank.lines.skip=FALSE,sep=",")
jul2012$Crime.type<- factor(jul2012$Crime.type)
jul2012<- jul2012[jul2012$Crime.type!="",]
jul2012$Crime.type<- factor(jul2012$Crime.type)
levels(jul2012$Crime.type)<- c("Anti-social behaviour","Theft/burglary","Criminal damage/arson","Drugs","other crime","Theft/burglary","Public disorder and weapons","Theft/burglary","Vehicle crime","Violence and Sexual Offences")

#plot
jul12Map<- ggplot(data=Norfolk,aes(long,lat)) + geom_polygon() + geom_point(data=jul2012,aes(Easting,Northing,colour=as.factor(Crime.type)),alpha=0.6) + scale_colour_brewer(palette="Set1",name="Category") + theme_bw() + scale_x_continuous(breaks=0,labels="") + scale_y_continuous(breaks=0,labels="") + theme(axis.ticks=element_blank(),panel.grid.major=element_blank(),panel.grid.minor=element_blank())

The  shp file and the crime stats file are linked, as well as the dbf, shx, and prj files.


texts: How to loop Grass command in Python and put variable in output file

texts:
I want to reclassify the same raster several times and then save the file with the value loop variable at the end ...
I cannot get it working.

    import sys
    import os
    import atexit

    import grass.script as grass

    def cleanup():
        pass

    def main():
        for i in range (2,4):
        out = "reclasstest" + str(i)
        rulesvar = "C:\Grass_data\datareclassfile") + str(i)
        grass.run_command("r.reclass",
                          overwrite = True,
                          input = "glg@PERMANENT",
                          output = out,
                          rules = rulesvar)

        return 0

    if __name__ == "__main__":
        options, flags = grass.parser()
        atexit.register(cleanup)
        sys.exit(main())


texts: How to create least-cost path between two polygons with GRASS?

texts:
I am completely a beginner in GRASS but I would like to use it to have the cumulative cost of moving between two polygons. More exactly, I would like to get the smallest cost from any point situated at the edge of one polygon (polygon A) to any point situated at the edge of another polygon (polygon B).

Until now, I used the function CostDistance and CostPath of ArcGIS like this:

CostDistance("start_polygon", "costs.tif", "", "outCostLink.tif")
CostPath("end_polygon", "outCostDist.tif", "outCostLink.tif", "EACH_ZONE","FID")

I would like to do the same thing with GRASS. I think that the function r.cost allows to do this. But the problem is that I don't understand how to specify in parameters the two polygons in GRASS ?

I have :

a cost raster where each cell has a cost value.
a shapefile for the first polygon
a shapefile for the second polygon

Have you got an example of how to use r.cost with two polygons ?

Thanks very much for your help.


texts: Using ArcPy and PyGRASS tools in same Python script?

texts:
I´m trying to write a Python script, that includes different ArcPy / ArcGIS commands and want to combine it with GRASS geoprocessing tools.

Unfortunately importing the GRASS libraries doesn´t work.

import grass.script as grass

ends in an error:

Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
ImportError: No module named script

I know that GRASS GIS installs it's own Python version. Isn´t there a way to combine the ArcGIS and the GRASS GIS installation of Python? I tried copying the pygrass module (or what I thought it might be it) to the C:\Python26\ArcGIS10.0\Lib\site-packages\ folder. I´m getting a different error message, but it´s still not working.

Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\Python26\ArcGIS10.0\lib\site-packages\grass\script\__init__.py", line 1, in <module>
    from core   import *
  File "C:\Python26\ArcGIS10.0\lib\site-packages\grass\script\core.py", line 38, in <module>
    gettext.install('grasslibs', os.path.join(os.getenv("GISBASE"), 'locale'))
  File "C:\Python26\ArcGIS10.0\lib\ntpath.py", line 96, in join
    assert len(path) > 0
TypeError: object of type 'NoneType' has no len()

I´m using ArcGIS 10 and GRASS 6.4.2 on a Windows 7 64 Bit machine.


texts:
I have created a view by joining a feature layer and 2-3 other tables in ArcSDE database using help from
James Richard Post.  The view gets created fine in both SQL Server 2008 and also in SDE.  In SQL Server it shows records and geometry but when it comes to ArcGIS it does not show table or geometry in ArcGIS Desktop 10.1.
.
It shows only the geometry in ArcCatalog but not the table in preview table.

Below is my code:

SELECT SPE.[OBJECTID] ,SPE.[BaseId], SPE.[PlayId] , B.[BasinName], SP.[ShaleName], AGS.[AgeS],  SP.[Ref] , SPE.[Shape] , SPE.[Shape_STArea__], SPE.[Shape_STLength__]
    FROM [dbo].[Play] SPE inner join [dbo].[Base] B on SPE.[BaseId]=B.[OBJECTID]
        inner join [dbo].[shale] SP on SPE.[Shale]=SP.[OBJECTID]
        inner join [dbo].[Age] AGS on SPE.[AgeId]=AGS.[OBJECTID]


texts: Esri Terrain-Tools v1.1 Sky Tool Cannot Save File

texts:
I am trying to use Esri Terrain-Tools v1.1 with ArcMap 10.5 to create a hillshade from either a 90 m SRTM DEM or 30 m DEM. I have tried saving this to a file geodatabase (FGDB) and .tif. In both cases the error reported is cannot save to "output format" GRID/TIFF
Traceback (most recent call last):
  File "C:\GIS_local\TOOLS\TerrainTools_v1_1\Scripts\SkyModel.py", line 81, in <module>
    skymodel(*args)
  File "C:\GIS_local\TOOLS\TerrainTools_v1_1\Scripts\SkyModel.py", line 63, in skymodel
    prevRaster.save(outRaster)
RuntimeError: ERROR 010240: Could not save raster dataset to C:\GIS_local\hillshade\sky1.tif with output format TIFF.

Failed to execute (SkyModel).
Failed at Mon Feb 26 17:18:37 2018 (Elapsed Time: 3 hours 46 minutes 44 seconds)

The output file path is to a local internal drive, scratch filepath is to local default GDB, and there should be more than adequate space to save the file. This seems to have been commonly reported in the past in comments on the ESRI website, but no answer was provided.
Does anyone have thoughts on what may be causing this error?


texts:
I use QGIS 2.18.0 on Windows 10.

I wanted to add a layer from PostGIS using the DB Manager. After choosing the "unique column" and the "geometry column" the DB Manager did not load the layer. However, without any error message or hint.
The reason why it did not load was that there were columns with the same name (duplicate column names).

Does DB Manager usually tell me about duplicate column names?

(or is this "feature" not implemented?)
(or is my QGIS version too old?)


texts: How do I convert EPSG:26986 coordinate pairs to EPSG:4326 from the MassGIS web service, the Linux command-line, or PHP?

texts:
I'm a GIS newbie, normally a PHP web applications developer.

I am trying to map all boundaries of the MA state representative districts.  So I just need to convert polygon coordinates from the MassGIS web service's output into Google-friendly lon/lat coordinates (the WGS84 SRS).  In GIS parlance I need to execute a "re-projection."

Below is the link to the XML (6MB, careful).  The default SRS is EPSG:26986.  You can add the GET parameter &srsName=26986 to the end of it and then toy with that number but I've tried plugging in the WGS84 equivalents in EPSG, which folks are saying are 4326 or 3875, and apparently MassGIS does not have those output formats available.

Here is the response XML:
http://giswebservices.massgis.state.ma.us/geoserver/wfs?request=GetFeature&version=1.0.0&typeName=massgis:GISDATA.HOUSE2012_POLY

Then I started playing with the Linux tool gdaltransform which seems promising, but when I plugged in over 100 of the original (EPSG:26986) pairs, it only spits back about four pairs (in the correct EPSG:4326 format).  So something is wrong there.  If tried dumping the coordinate pairs in a text file of source data, but then the gdaltransform complains that it is not a supported file format.

Obviously the easiest thing would be if someone has some experience with the MassGIS web service and can tell me the right parameter to get these coordinates in the format I need them.  But another fine solution would be to help me find a simple, quick conversion script or utility (PHP or Linux) that I could integrate into my web app behind the scenes.


texts: Extract intersections of 2 shapefiles with field subset in QGIS

texts:
I have 2 polygon shapefiles: A (with integer field "a_id") and B (integer field "b_id").

My goal is to extract intersections of polygons in layer A with polygons in layer B where attribute "a_id" is equal to attribute "b_id" + 55.

I'm trying to implement it with the QGIS graphical modeler, but "Select by expression" algorithm supports only one layer for running expressions, and "Extract by location" has no attribute filter.


texts: Blank Polygon issue using QGIS Contour Plugin

texts:
I am using QGIS with Contour plugin, want to generate filled contour polygons with user selected contour levels.
However there are some errors showing blank polygons.
Is it bug or is there any way to fix it?
Parameter is here
{ 'ContourInterval' : 1, 'ContourLevels' : '0.0 0.2 0.4 0.6 0.8 1.0 1.5 2.0 3.0 4.0 5.0 6 7 8 9 10 12 14 16 18 20 25 30 35 40 50 60 70 80 90 100', 'ContourMethod' : 4, 'ContourType' : 1, 'DuplicatePointTolerance' : 0, 'ExtendOption' : 0, 'InputField' : '\"RN-60m\"', 'InputLayer' : 'file:///[path].csv?encoding=UTF-8&delimiter=,&xField=LON&yField=LAT&crs=epsg:4326', 'LabelDecimalPlaces' : -1, 'LabelTrimZeros' : False, 'LabelUnits' : '', 'MaxContourValue' : None, 'MinContourValue' : 0, 'NContour' : 400, 'OutputLayer' : 'C:/[path]/out.geojson' }


texts:
I have about 40 QGIS2 projects. If I open them with QGIS3 labels and svg icons are not displayed. I don't know if there is any other feature that is not working when upgrading to QGIS3, but I suppose there must be more of them.

For labels, for example, I can see no configuration in QGIS3, as if labels config was not defined in the layer. Seems that QGIS3 does not understand that part of the layer configuration saved in QGIS2.

Is there any tool/procedure to upgrade a QGIS2 project to QGIS3?


texts: Why is my geoserver css style not working when I add multiple filters?

texts:
I am using Geoserver and the CSS Style plugin. The following code works fine and shows an 'L' when the feature is 'Left' and the '*' filter catches anything else and shows it as an X.

* {
 mark: symbol('ttf://DejaVu Sans#0x2A02');
 mark-size: 24;
}

[left_right = 'Left'] {
  mark: symbol('ttf://DejaVu Sans#0x004C');
  mark-size: 18;
}

However, when I add a filter for 'Right' then 'R' and 'L' show up but then the catch-all * filter stops working.

* {
 mark: symbol('ttf://DejaVu Sans#0x2A02');
 mark-size: 24;
}

[left_right = 'Left'] {
  mark: symbol('ttf://DejaVu Sans#0x004C');
  mark-size: 18;
}

[left_right = 'Right'] {
  mark: symbol('ttf://DejaVu Sans#0x0052');
  mark-size: 18;
}

I have changed the * rule to the lines below and the symbol still doesn't show up.

[left_right <> 'Right'] [left_right <> 'Left']
[left_right <> 'Right' or left_right <> 'Left']
[left_right = '']

I have also left the * filter alone and added [left_right = ''] as a fourth filter and it doesn't show up either.

I guess I am missing a basic concept of styling rules but I can't figure out what it is. It looks simple to me - I should be able to change my symbol based on the feature attribute. I am not symbolizing unrelated things so I don't see why I would need separate styles.

I have not tried changing the @mode rules since I do not think this should generate an overly complicated SLD style.

What can I change or am I doing wrong?

Full version of the code minus the filter for 'Right':

* {
 mark: symbol('ttf://DejaVu Sans#0x2A02');
 mark-size: 24;
}

[left_right = 'Left'] {
  mark: symbol('ttf://DejaVu Sans#0x004C');
  mark-size: 18;
}

:symbol {
 fill: black;
}

[@scale > 2500] {
  mark-size: 6;
}


texts: How to install QGIS 2.2 instead of 2.4 to get a working OpenLayers plugin again?

texts:
I use ubuntu . I don't know how to return installed qgis 2.2
I try uninstall qgis 2.4 and try find how to install old version qgis 2.2
Please .. tell me how to install old version qgis 2.2 for ubunut 12.4


texts: Multiple values for condition in IIF statement in Arcade

texts:
I am trying to create an immediate calculation attribute rule for a feature class in ArcGIS Pro (that I hope carries over when published as a feature service on AGOL).
I understand the basic syntax: IIF(condition, trueValue, falseValue). So if I wanted to create a rule that would take a value in my fictitious FoodName field and auto-calculate my FoodType field, I could do some thing like...
var food = $feature.FoodName
IIF(food == 'Carrot', 'Vegetable', 'Fruit')

This would, whenever the value 'Carrot' appears in the FoodName field, calculate the FoodType field to 'Vegetable.' My question: How do you pass a list of values to the condition in these statements? This has been surprisingly not straight forward to find online.
For instance, I would like the values 'Carrot', 'Celery', 'Parsnip', and 'Broccoli' to all to auto calculate to 'Vegetable.' Is it possible to pass a list to the condition, rather than writing multiple IIF statements? What is the preferred syntax?


texts: Style points using color gradients with QGIS

texts:
I need to style a shapefile with a specific color gradient, but considering absolute values.
It's point shapefile with values between 0.19 and 0.70 (considering they're in an absolute scale between -1 and 1).

I'm using QGIS 2.14.3.

I've build a color map using Matlab (.svg format based on the cpt-city) but the problem is that the color steps inside the .svg file are divided with percentages.
I can't give the exact value the style/color I need, unless I do that manually, one by one.

I did the trick when using raster data: I've build color ramp with the exact steps and the correspondent color.


texts: Setting up indexes for PostGIS distance queries

texts:
I'm building an application that is supposed to query and return every Record in a table that is X kilometers away from PointX. Records and PointX's positions are determined from (long/lat) information provided by Google Geocode API.
I'm new to PostGIS. After a quick research, I've found this question. The answer seems to be along the lines of:
SELECT *
FROM your_table
WHERE ST_Distance_Sphere(the_geom, ST_MakePoint(your_lon,your_lat)) <= radius_mi * 1609.34

The problem is that even though I'm only getting started at GIS, when I look at the above query, I can't possibly imagine how this can use an index. There are 2 function calls. I imagine the table being scanned for every Record.
Does PostGIS have any index type capable of making the above query performant and, if not, what would be the recommended approach doing what I need?


texts:
When I print as PDF a work done with a shape point and OpenLayers plugin I notice that the point are shifted

In the dual screen image you can see the PDF exported image (left) and the QGIS working one (right). You can notice that the point are a little shifted in the exported image and I am not able to understand the reason.

In the image you can see the same yellow street and the point on the street (right image) and shifted in the exported image

I use QGIS 2.6.1


texts: QGIS2.4.0: Cannot import qgis module in python "from qgis._core import * \n DLL Load Error"

texts:
Having some trouble importing qgis into python. I'm very new to python so it could be a simple mistake. The path seems to be correct. It gets this far:

    File "<stdin>", line 1, in <module>
    File "qgis\__init__.py", line 36, in <module>
    from qgis.core import QgsFeature, QgsGeometry
    File "qgis\core\__init__.py", line 1, in <module>
    from qgis._core import *
    ImportError: DLL load failed: The specified module could not be found.


texts: Field Calculator: Split address into house number, street name

texts:
I have an ADDRESS field in a table with values like:

1 MAIN ST
20 SIDE AVE W
300 JOHN DOE JR RD S

I want to use the field calculator to split the addresses and insert the parts into HOUSE_NUMBER and STREET_NAME fields.
Example:

300 into HOUSE_NUMBER
JOHN DOE JR RD S into STREET_NAME

I also want the tool to check if the value that will be put into HOUSE_NUMBER is a valid integer (and ignore it if it isn't).

How can I do this?


texts: PostGIS query to split a polygon by n other polygons

texts:
I'm trying to write a PostGIS query to do the following:

Get polygons from a table that overlap geometry A above (in this case, B and C).
Get the intersecting area.
Split the intersecting area by the overlapping polygons (in this case, returning D and E).

In this case there are 2 overlapping polygons, but I need to support an arbitrary number.

I'm stuck on how to get started with this. I know I can use ST_Intersects to get the intersecting polygons and ST_Intersection to get the overlapping area, but I'm not sure the splitting will work. It looks like ST_Split only accepts two geometries Can anyone think of how I could get this to work with n polygons?


texts: QGIS Points are not retaining original attributes

texts: I am trying to map birth and death locations for several hundred people in the US and Europe (for a total of 532 points in one layer). Twice now I have discovered that somewhere during this process (of manually going to each location and placing a point and filling in the attribute table) the original attribute information has been switched. For example, attributes given to a point in England are now associated with a point in North Carolina. I don't understand why this is happening. Why would attribute information move around like that? It makes no sense. Please help. I've been working on this for months. I just might lose my mind if I have to do this all over again.


texts: Why does ArcGIS 10 Convert Time Field geoprocessing tool return  when using custom format?

texts: I have dates in my attribute table that are stored as strings - they look like 02/22/2012 (mm/dd/yyyy) - and I would like to convert them to Date type. When I try to input a custom date format in the Convert Time Field geoprocessing tool it returns  in all the records and I am not sure what I am doing wrong. I've tried MM/dd/yyyy and MMddyyyy - is there something I'm missing about setting up a custom date format? Or, is there another way I can do this?


texts: How to draw all possible lines between a set of points in QGIS?

texts:
I have a points layer in QGIS defined by a csv file in the following format:

group,name,lat,long
1,A,…,…
1,B,…,…
1,C,…,…
2,A,…,…
2,D,…,…
3,E,…,…
3,F,…,…

(lat and long are coordinates in WGS84)

How can I draw all of the possible lines between all of the points that are in the same group? So the above data would yield a triangle with vertices at A, B and C, a line with end points A and D, and another line with end points E and F.

I have tried a number of solutions suggested on these boards (PointsToPaths plugin, Points2One plugin, Connect Points plugin, GRASS v.to.lines, SAGA Convert lines to points) but all of these only draw the minimum number of lines to connect the points. This works fine when only two points need to be connected, but fails when there are more than two points. In the example above, group 1 comes out as a line from A to B to C, without the direct line from A to C. Sometimes there are up to eight points in the group, which should yield 28 lines.


texts: Joining table to shapefile without creating new shapefile in ArcPy?

texts:
I am trying to join a csv file to a shapefile using a unique field and a python script. Then I want to append one of the fields from the csv file to the shapefile and remove the join.

Is there a way to do this without creating a new "output" shapefile? The shapefile that I am working with is a very large shapefile in a .gdb and copying all of its information to a new shapefile, deleting the old shapefile, and moving the copied information back to the original location would be non-ideal.


texts: Export/Import a layer as QLR in QGIS 3.4 using PyQGIS

texts:
I am trying to find a way to export a layer in a QLR file in PyQGIS. I want to do the same as right-click and export to layer definition file.

I know that I have to use QgsLayerDefinition.exportLayerDefinition and it needs a QgsLayerTreeNode parameter. But I don't know the steps to take a QgsVectorLayer and put it in a QgsLayerTreeNode.

I also need to do the reverse, load a QLR file and show the layer in QGIS.


texts: Conditional Labels, Multiple Labels ArcGIS

texts:
I'm attempting to write a code to show 6 fields that all have some Null values. The fields will have accompanying text labels that are not a field, and each one needs to be on a different line.

The 6 fields are [TF] [INV] [INV_N] [INV_S] [INV_E] [INV_W]
The desired output would be
"TF [TF]"
"INV [INV]"
etc.

I have tried to do research and can't seem to find anything with this many values. The vb-script code I attempted to write is:

    Function FindLabel ( [TF], [INV], [INV_N] , [INV_S] , [INV_E] , [INV_W]  )
     FindLabel = "TF " & [TF] & vbCrLf
      If IsNull( [INV] ) Then
     FindLabel = ( [INV_N] )
      If isNull ( [INV_N] ) Then
     FindLabel = ( [INV_S] )
      If isNull ( [INV_S] ) Then
    FindLabel =( [INV_E] )
     If isNull ( [INV_W] ) Then
    Else
     FindLabel = "INV " & [INV] &vbCrLf
    Else
     FindLabel = "INV N " & [INV_N] & vbCrLf
    Else
     FindLabel = "INV S " & [INV_S] & vbCrLf
    Else
     FindLabel = "INV E " & [INV_E] & vbCrLf
    Else
     FindLabel = "INV W " & [INV_W]
    End If
    End If
    End If
    End If
    End If
    End If
   End Function

My coding experience is 1 Java class my freshman year of college so I'm pretty lost when it comes to this.

If the field has a Null I would not want the text to be shown nor for there to be an empty space


texts: Multipatch to polygon/polyline conversion in ArcGIS Pro?

texts:
I have a file captured from Bentley software already converted to .shp format. I would like to have polygon or polyline to be able to simply edit it but I don't know how to convert Multipatch to polyline/polygon without loosing Z-values. I tried Multipatch Footprint but as the result I have 2D layer. As well I tried to generate points using Feature Vertices to Points and it works, but my goal is to create TIN model and here I don't want to include all points (some of them are below desired surface).

Does ArcGIS Pro have another tool which could help me?


texts: How to align ASTER GDEM with Google Maps?

texts:
I'm using ASTER GDEM V2 to pick heights at lon-lat coordinates obtained from Google Maps but I found that at some locations - particularly San Francisco: https://goo.gl/maps/w7p1K - the DEM seem to be horizontally shifted by about 100m in NNW direction which is more than declared ASTER GDEM precision.

I'm using QGis 2.0.1 under Ubuntu.

Reproduce:
Open freshly downloaded/unmodified ASTGTM2_N37W123_dem in QGis. The tile should open in it's native GCS EPSG:4326. Note coordinates of outline of San Francisco Treasure Island and display those in Google Maps.

The question is: Am I missing something about coordinate systems and CS references involved? As far as I know Google's lon-lat is also EPSG:4326.

Or is that a bug in ASTER GDEM's georeferencing?

Further analysis:

The shift is not homogeneous. It is definitely smaller or negligible in Czech Republic (ASTGTM2_N50E014_dem) but the area lacks such strong features as coastlines to tell for sure.

In San Francisco I compared ASTER to SRTM and SRTM does not seem to exhibit described shift.


texts: Importing cvs file with coordinates into QGIS

texts:
I got a problem. I got some coordinates in Excel. They are written with comma 14459091,1010 ex.

When I try to import them as csv-file into QGIS, it will only create a table and no points. When I look in the properties it says that my "x" and "y" columns are strings.

I use a semicolon to separate the fields.


texts: Can't upgrade file geodatabase from 10.0 to 10.5

texts:
I want to examine / use ArcGIS Desktop 10.5.1 on some data that was all created in 10.0 in a file geodatabase.  However, when I try to "Upgrade Geodatabase" (either through the geoprocessing tool by right-clicking Properties -> Upgrade Geodatabase), 10.5.1 thinks that the geodatabase is "up to date" at 10.0.

I know everything's not fine because trying to manipulate the data every so often throws me this error:

Any suggestions?


texts:
I have a google Layer map on top of which I have added a raster layer of Image and then mapped image coordinates from Google layer.

I made transparent (75%) image map layer and now mapping and creating polygons on top of image layer (Added new shapefile layer)
I drew few polygons and then saved it.

When I reopened the file, clicked on toggle editing, I could see my polygons, but I am not able to add/edit/delete existing/new polygons in that shapefile layer. (All icons are disabled)

Is this expected behavior?

My old question on what exactly am I doing.
How to create/obtain district data for Singapore?


texts:
It appears that I have corrupted / overwrote the index planet_osm_nodes for my Open Street Map installation.  I would like to recreate it, but am not sure of the access method (b-tree, GIST...) which osm2pgsql initially created it as.  I am using postgres.

I have tried this:

  create index planet_osm_nodes on public.planet_osm_nodes USING gist (id)

But get errors

Can anyone let me know how best to recreate this index?


texts:
I'm relatively new to programming and ArcObjects. I'm trying to simply calculate the length of a line (test.shp) using ArcObjects and Java. I've tried the code below with no luck. This is simple in ArcMap because one just needs to use 'calculate geometry' on the field within the attribute table. But I need code to do this.

Any suggestions would be greatly appreciated!

FeatureClass fc1 = null;
ShapefileWorkspaceFactory shapefileWorkspaceFactory = new ShapefileWorkspaceFactory(); Workspace shapefileWorkspace = new Workspace(shapefileWorkspaceFactory.openFromFile("c:/PFRAS/TEMP", 0));
fc1 = new FeatureClass(shapefileWorkspace.openFeatureClass("test.shp"));

Field fe = new Field();
fe.setName("Dist");
fe.setType(esriFieldType.esriFieldTypeDouble);
fe.setScale(2);
fe.setPrecision(16);
fc1.addField(fe);

IFeature feature = fc1.getFeature(0);
Line line = new Line();

//Assume that the underlying feature is a Line.
IGeometry geom = feature.getShape();
    if (geom instanceof Line)
    line = (Line)geom;
    double length = line.getLength();

int fieldIndex = fc1.findField("Dist");

FeatureCursor fCursor = new FeatureCursor(fc1.IFeatureClass_update(null, false));
        Feature feature2 = (Feature) fCursor.nextFeature();
        while (feature2 != null)
        {
         feature2.setValue(fieldIndex, length);
         fCursor.updateFeature(feature2);
         feature2 = (Feature) fCursor.nextFeature();
        }
        fCursor.flush();
        Cleaner.release(fCursor);


texts: change projection in an ascii (x,y,z) file

texts: I am having an ASCII file (lat,long,depth) in the Mars projection. I need this file to be converted into WGS-84 projection. How can I do the conversion for the ASCII files.


texts: How to cluster points into clusters of a maximum diameter in PostGIS?

texts: I'm working with transit stops in GTFS. I want to group together stops (N~8000 for a medium sized agency in the USA) which are proximately related in order to reduce the number of dimensions but also group together stops that are a reasonable distance from each other: e.g. transfers. K-means requires an a priori specification of the number of clusters, whereas I would like the cluster diameter to be specified by the user. A bottom-up hierarchical clustering algorithm which stops at a specified diameter is what I need.


texts: QGIS is not able to find GRASS. Choose the path of installation folder of GRASS

texts:
I've downloaded and then correctly installed QGIS 1.7.4 for windows (Win 7 64bit is my version) in C:\Programs\Quantum GIS.

But everytime I start QGIS this popup appear: "QGIS is not able to find the installation of GRASS. Choose the path of the GRASS installation folder".

I've chosen every installation folder in Quantum GIS but the program not recognise it!
So if I don't choose the folder, I can't use GRASS.
In the installation folder there are a file called "grass64.bat", but not "grass64.exe". Could it be the problem?

I've installed QGIS also in other folder and D:\Programs\... but without results.


texts: Changing Color and Opening Label on Hover

texts:
I'm initializing GeoJSON markers in Leaflet, and I'm trying to change the color and open a label on mouseover. Upon mouseout, I want to change back to the default color and close the label. However, the layer only handles one event listener -- in the code below, only the mouseover event is fired. When the mouseover listener is commented out, only then will the mouseout event be active. What am I missing in my code to handle more than one listener?

// load geojson data with default variables shown
var mapLayer = L.geoJSON(data[0]['features'], {
    pointToLayer: function (feature, latlng) {
        var marker = L.marker(latlng, {icon: defaultMarker});
        return marker;
    },
    onEachFeature: onEachFeature
}).addTo(map);

// on each marker, do the following
function onEachFeature(feature,layer) {
    layer.bindTooltip(feature.properties.label);
    layer.on('mouseover', function(e){
        layer.setIcon(highlightMarker);
        layer.openTooltip();
    });
    layer.on('mouseout', function(e){
        console.log("mouseout");
        layer.setIcon(defaultMarker);
        layer.closeTooltip();
    });
}


texts: Is there a way to change environment settings of a dataset without creating a new one?

texts:
I am trying to change the environment settings of an existing dataset. I found a workaround by first changing the environment settings in the SDE-database and then create a new dataset with now correct settings. Now I copy&paste my features from the old to the new dataset. But I guess copying is rather dangerous because some database information might be lost? So I would like to avoid it.
The values I would like to change are:
XY-Tolerance
XY-Resolution
Z-Tolerance
Z-Resolution
Z-min.Domain
M-Tolerance
M-Resolution
I am aware that with changing these values the features will grow/shrink/move a tiny bit.Any help is appreciated!


texts: Point to line distance using dist2Line and 'sf' library

texts:
I need to calculate the distance from a set of points to the closest line. I followed the approach outlined here.

dist_plsscentroid_rail <- geosphere::dist2Line(p = st_coordinates(plss_centroidWgs84.sf), line = st_coordinates(railWgs84.sf)[,1:2])

The input line file looks like the image on the left. The figure on the right shows a subset of the returned distances (in blue) and the line geometry (in red). As you can see, the geometry of the lines is different after calculating the distances. In particular, several additional segments are included. Any ideas on what might be happening here?

Thanks for your help!


texts: Append element to an array in Google Earth Engine

texts:
I am trying to write a for loop to append elements dynamically to an array. However, the dimensions of the array change on every loop run, throwing off error.

var arr = 1
for (var i=0;i<4;++i){
  var newval = somefunction()
  arr = ee.Array.cat([newvar,arr],0)
}
print(arr)

Error:
Array.cat: Arrays must have same dimensionality. Expected 0, found 1 for array 1.

EDIT:  The new value to be appended is coming out of a for loop, so I can't use a list.map or iterate function


texts: PostGIS connect points to 1 closing/bounding line

texts:
I have this points represented as geom:

And when I apply: ST_MakeLine(points.geom) I got this:

When I need something like closing/bounding line (and even better the whole area as 1 geom):

This is how I get the points for example:

SELECT DISTINCT (ST_DumpPoints(buildings.geom)).geom  --AS geom
FROM buildings
WHERE ST_DWithin(buildings.geom, ST_GeomFromText('POINT(lon lat)', 4326), radius)

EDIT 2:

Applying ST_ConcaveHull(ST_Collect(vertices.geom), 0.99, false) resulted almost the solution:

The reason is as documented here: "he basic approach is that it first creates a convexhull of the geometry and then uses the ST_ClosestPoint function introduced in PostGIS 1.5 to cave in the hull to transform it into a concave hull. That's the basic idea". We can see that the points not connected properly are those where the issue with the closest point is Misleading...


texts: Add multiple fields to attribute table by python, with a prompt window for parameters

texts:
I know how to add field to arcmap shapefile using Python:

import arcpy
arcpy.env.workspace = "H:\pythonexercise\NZL_rds\NZL_roads.shp"
arcpy.AddField_management("NZL_roads", "LENGTH", "TEXT","","","50")

But how to make the python code prompt a window for "field", "data type" and "length"? Also, what if I want add many fields at once like below:

field "SPEEDLIMIT" data type as "num" and length "3".
field "BREADTH" datatype as "num" and length "4"
field "LANE" datatype as "num" and length "4"

It should be dynamic, not fixed like my code.


texts:
const mapnik = require('mapnik');
const sphericalmercator = require('@mapbox/sphericalmercator');

mapnik.register_default_fonts();
mapnik.register_default_input_plugins();

const TILE_HEIGHT = 256;
const TILE_WIDTH  = 256;

const SRS_MERC = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over";

const mercator = new sphericalmercator({
    size: TILE_WIDTH
});

var map = new mapnik.Map(TILE_WIDTH, TILE_HEIGHT, SRS_MERC);
map.background = new mapnik.Color(255, 255, 255, 0);
var style = `
<Map>
  <Style name="census_wards">
    <Rule>
      <PolygonSymbolizer fill="#00ffff" fill-opacity="0.75"  />
      <LineSymbolizer stroke="#000000" stroke-width="0.5" stroke-opacity="1.0"  />
    </Rule>
  </Style>
</Map>
`;
map.fromStringSync(style);
var postgis = {
  'type': 'postgis',
  'dbname': 'test_db',
  'host': '127.0.0.1',
  'user': 'test_map',
  'password': 'TestUser',
  'table': '(SELECT * FROM z4.census_wards) AS census_wards',
  'geometry_field': 'the_geom_webmercator'
};

var layer = new mapnik.Layer('census_wards', SRS_MERC);
layer.datasource = new mapnik.Datasource(postgis);
layer.styles = ['census_wards'];
map.add_layer(layer);

var box = mercator.bbox(2940, 1846, 12, false, '900913');
map.zoomToBox(box);
console.log('box = ' + box);

var im = new mapnik.Image(TILE_WIDTH, TILE_HEIGHT);
console.log('rendering');
map.render(im, function(err, im) {
    console.log('inside render');
    if (err) throw err;
});
process.exit(0);

Control never reaches the "console.log('rendering');" line. Couple of days ago something like this was working. I am assuming that I am missing a very simple but important part due to my lack of experience with node and mapnik.

Apart from the console.log output, there is no response from the program. No errors, no other output.

Any pointers on how to get this working?


texts:
I'm struggling to understand VSIMEM. We're saving it in memory, similar to a "MEM" or "Memory" dataset, but it has a path under vsimem. If we're saving it to a path, then how is it in memory?
I've read the documentation as well as searched around, but am struggling to understand the concept. It also seems to introduce "Unlinking" after we're done, whereas with "MEM" or "Memory" datasets, it would operate like Shapefiles/TIFF datasets and we would just set it to None after we're done with it. Is Unlinking the equivalent of setting to None?


texts: Zoom to centroid of polygon from form-input in Leaflet

texts:
I want to  write the name of the polygon in a form and then click a button to zoom or to flyto to the centroid of the polygon. I've writen a request in PHP that give the centroid of the polygon.

function localiser5() {
  var val6 = document.getElementById('val6').value;
  console.log(val6);
  var x;
  var y;
  var z;
  $.ajax({
      type: 'POST',
      url: 'php/survoler.php',

      data: {
        val6 : val6
      },
      success:function(data){

     map.flyTo(new L.LatLng(x, y),12);
      },
      error:function(response){
        console.log(response);
      }
    })

}

PHP:

$db = pg_connect("host=localhost port=5432 dbname=d3topo user=postgres password=khadija") or die('could not connect :'.pg_last_error());

if(!empty($_POST['val6'])){
  $nom_table='points_rattachement';
  $champ_geometrique = 'geom';
  $val6 = $_POST['val6'];

  $sql = "SELECT st_x(st_centroid(geom)) as x, st_y(st_centroid(geom)) as y
          FROM plan_situation where nom_plan='".$val6."'";

  $rs = pg_query($db, $sql);
  if (!$rs) {
    echo "Erreur SQL survenue!.\n";
    exit;
  }
  $testRow = pg_fetch_row($rs);
  echo $testRow[0];

}

HTML:

<input  type="text" placeholder="Nom plan" name="val6" id="val6" />
<button type="submit" onClick="localiser5()" >FlyTo Polygon </button>


texts:
is there a way to order the returned results of a query task?

this:

query.where = "COD > 0";

works fine, but this:

query.where = "COD > 0 order by COD";

throws:

dojo.io.script error undefined

Please note that our server has standardized queries enabled and we cannot change that.


texts: Is ST_GeomFromText better than providing direct geometry?

texts:
I have been working with postgis recently,and in my query if I use ST_GeomFromText it execute faster than running a sub-query to get geom.

I thought ST_GeomFromText will be more expensive but after running many tests every time I got the result faster, my question Is there any explanation behind this? because for me getting the geom directly in sub-query is better than getting geom as text then added as GeomFromText.

Thanks


texts: Converting DNs to radiance and reflectance for Planet Labs Imagery

texts: I am trying to convert pixel DNs to reflectance from a Planet BGRN Analytic Scene (from 2016, I can't tell what level it is processed to).  However, I can't find info on the exo-atmospheric irradiance of the 4 bands.  Is this available somewhere?


texts:
My goal is to add Callouts (Manhattan type) to labels in a point layer with PyQGIS.
Here is what I got so far, which basically renders the labels without the callouts:
font = QgsTextFormat()
font.setFont(QFont('MS Shell Dlg 2'))
font.setSize(12)
p_layer = QgsPalLayerSettings()
p_layer.setFormat(font)
p_layer.fieldName = 'label'
p_layer.enabled = True
labels = QgsVectorLayerSimpleLabeling(p_layer)
layer.setLabeling(labels)
layer.setLabelsEnabled(True)
layer.triggerRepaint()

What I have tried:
I read in the documentation that the QgsPalLayerSettings class has a method named setCallout, which if I understand correctly, takes an object of QgsCallout class as argument.
When I try to start setting/using the QgsCallout() class, I get the following message in the QGIS Python console:

TypeError: qgis._core.QgsCallout represents a C++ abstract class and cannot be instantiated

According to this answer in StackOverflow:

"The error means there are some methods of the class that aren't implemented. ..."

Question: How to add Callouts to labels with PyQGIS? I'm using QGIS version 3.14.1


texts: Set Input Directory Path for ArcMap Geoprocessing Tool

texts:
How do I set the input file directory in a arcmap geoprocessing tool?

I want to be able to input a file from a specific directory without having to browse for the directory every time I want to run a geoprocessing tool.


texts:
I'm trying to find a plug-in on leaflet in which the user will enter the address and the html (aka the plug-in) will geocode the location and will return a latitude and longitude. So I have  try the Leaflet-Control Geocoder,but I am not able to the lat,long.

Does anyone have an example?


texts: Applying less than and greater than threshold in image segmentation in Google Earth Engine

texts:
I am trying to segment the images based on the threshold values. For instance, I want to put a threshold range between > 800 and < 1200, but with my existing code in Google Earth Engine, I can put the only one function at a time that is greater than, below is my code, how to define greater than and less than at the same time.
 // Threshold the thermal band to set hot pixels as value 1, mask all else.
var hotspots = s2a.gt(3500) // i want put here both grater than and less than function
  .selfMask()
  .rename('hotspots');

// Display the thermal hotspots on the Map.
Map.addLayer(hotspots, {palette: 'FF0000'}, 'Hotspots');


texts: GeoJSON properties not displaying in popup-leaflet GeoDjango

texts:
I want to show the markers in different colours based on column(rastervalues) values in leaflet. My geojson data is

{'geojson': "[{'type': 'FeatureCollection', 'features': [{'type': 'Feature', 'ge
ometry': {'type': 'Point', 'coordinates': [74.886779, 33.571307]}, 'properties':
 {'rastvalues': 9, 'id': 1, 'species': 'Poa pratensis L.'}}, {'type': 'Feature',
 'geometry': {'type': 'Point', 'coordinates': [75.015556, 33.588333]}, 'properti
es': {'rastvalues': 9, 'id': 2, 'species': 'Poa pratensis L.'}}, {'type': 'Featu
re', 'geometry': {'type': 'Point', 'coordinates': [75.06231, 33.592838]}, 'prope
rties': {'rastvalues': 9, 'id': 3, 'species': 'Phleum alpinium L.'}}, {'type': '
Feature', 'geometry': {'type': 'Point', 'coordinates': [75.196944, 33.550556]},
'properties': {'rastvalues': 12, 'id': 4, 'species': 'Pinus roxburghii Sargent'}
}, {'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates': [75.224722,
33.540278]}, 'properties': {'rastvalues': 12, 'id': 5, 'species': 'Pinus roxburg
hii Sargent'}}, {'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates':
 [75.103889, 34.485833]}, 'properties': {'rastvalues': 3, 'id': 6, 'species': 'P
inus roxburghii Sargent'}}, {'type': 'Feature', 'geometry': {'type': 'Point', 'c
oordinates': [75.143611, 34.536667]}, 'properties': {'rastvalues': 3, 'id': 7, '
species': 'Phleum alpinium L.'}}, {'type': 'Feature', 'geometry': {'type': 'Poin
t', 'coordinates': [74.329145, 33.818276]}, 'properties': {'rastvalues': 9, 'id'
: 8, 'species': 'Pinus roxburghii Sargent'}}]}]"}

My javacript code is

function main_map_init (map, options) {
             var promise  = $.getJSON('{% url "datapot" %}');
            // Download GeoJSON via Ajax
            promise.then(function(data) {
            var allbusinesses = L.geoJson(data);
            var cafes = L.geoJson(data, {
            filter: function(feature, layer) {
                return feature.properties.rastvalues == "3";
                },
                pointToLayer: function(feature, latlng) {
                return L.marker(latlng, {

                }).on('mouseover', function() {
                    this.bindPopup(feature.properties.species).openPopup();
                });
            }
        });

             var others = L.geoJson(data, {
            filter: function(feature, layer) {
                return feature.properties.rastvalues != "3";
            },
            pointToLayer: function(feature, latlng) {
                return L.marker(latlng, {
                }).on('mouseover', function() {
                    this.bindPopup(feature.properties.species).openPopup();
                });
            }
        });
         map.fitBounds(allbusinesses.getBounds(), {
            padding: [50, 50]
        });
        cafes.addTo(map)
        others.addTo(map)

        });
        }
    </script>

In chrome debugger, for the line filter: function(feature, layer) {
I'm getting the following value

feature = Object {type: "Feature", id: 1,properties: Object, geometry: Object}, layer = undefined

Properties namely species and rastvalues are not appearing.

How to solve this?

solution is

def extract_raster_points(request):
    conn = psycopg2.connect(dbname="geodjango",host='localhost',user='postgres',          password='postgres', port=5433)
    cur=conn.cursor()
    dict_cur = conn.cursor(cursor_factory=psycopg2.extras.DictCursor)
    res=dict_cur.execute("""SELECT row_to_json(fc) FROM ( SELECT      'FeatureCollection' As type, array_to_json(array_agg(f)) As features FROM (SELECT 'Feature' As type, ST_AsGeoJSON(lg.geom)::json As geometry ,row_to_json((SELECT l FROM (SELECT id, species,rastvalues) As l )) As properties FROM pft_account As lg   ) As f )  As fc;""")

    points=dict_cur.fetchone()
    #datapot = json.loads(str(points))
    datapot = json.loads(json.dumps(points))
    datapotg=datapot[0]
    print(datapotg)
    return JsonResponse(datapotg,safe=False)


texts:
https://chocolatey.org/packages/QGIS/2.6.1#individual.
I used the above link to download QGIS 2.6.1 but in this QGIS 2.6.1 Desktop is no where to be seen after installation.
Please suggest.


texts:
I have install on a Debian 9 a qgis-server, python-qgis, apache2 and libapache2-mod-fcgid. Apache service is runnning.

I have a folder uder usr/lib/cgi-bin/ but when I request http://localhost//cgi-bin/test/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&MAP=projet_test.qgs

I have edited the conf file to enable fcgi.

Everything works on qgis-server 2.14 but not on qgis-server 3.

Any ideas?


texts: Search plot by total area of two adjacent plots

texts:
I've got some QGIS data where each plot is described with some attributes e.g.  Shape_Area which is the field describing the area of plot in m2.
Normally I search for plot using the expression e.g. Shape_Area >= 1000 AND Shape_Area <= 1100, however this time I need to find two connected plots given their total area.

What's the best way to approach this problem?


texts: arcGISMapImageLayer works in WPF, not in Xamarin.Forms

texts:
I am trying to get an overlay working in Xamarin Forms (all code in shared folder) with Esri.ArcGISRuntime SDK v100.1.0 for .NET, testing with Android as start-up project.  I load an ESRI basemap, then add a busmap/transitmap (from our own GIS server) as a layer on top of that. When I monitor the layerviewstatus, it goes from “loading” to “error”.   I do get the basemap to display, but nothing is overlayed.

If I try the sample layers hosted by ESRI, like the census overlay, it works fine.

I can get the same code to display the layers in WPF, but not in Xamarin Forms.


texts: How to make a layer or group invisible in the Layers Panel?

texts:
I have a number of layers in the Layers Panel inside a group but I don't want them shown (there's quite a lot so it looks untidy). These layers are also joined to several others.

Is it possible to make these layers invisible to the user (i.e. not shown in the Layers Panel, it is not about setting the visibility of the layer) but can still be joined to other layers? Or perhaps make the group invisible?

Closest I found was this post but hoping there's a more direct approach:

Showing only some layers in QGIS legend?


texts:
I need to clip a large and very detailed shapefile in QGIS showing environmental data similar to a temperature map.

Is there something like the very simple crop tool like in Adobe Photoshop where you can drag a rectangle over an area of interest and crop the shapefile to this area only?

I do have an additional shapefile of neighborhoods/districts that I could also use for cropping but I do not know how to specify "extract a shp/kml only within the borders of district X, row xxxxx in the attribute table of districts.shp"

How can I extract the data primitively within a rectangle of choice or per district?

I need smaller pieces from the large map that has been provided by an environmental geo-institute to be implemented into an mobile Application (Android) to view the quite heavy data in a lighter way.

I made some screenshots of the data and posted them here: http://merglindev.wordpress.com/2012/10/22/extracting-districts-from-big-shapefile-data/.

You can see the districts and the data layers. I need to have every district as a separate file / separately export every dataset per district to KML.


texts: How to build and pass remap for floating point raster without getting Reclassify ERROR 010255?

texts:
I am receiving the following message when executing my python script:

  ERROR 010255: String remap must be used with a string field.

My input raster is floating point and I would like to reclassify it in python using gp.Reclassify_sa.

I build my remap up as a string, but this doesn't seem to work.

How do I build and pass my remap for a floating point raster?

The relevant code extract is as follows:

remap1 = "0 1.25 100000"
remap2 = "1.25 10 200000"
full_remap = remap1 + ";" + remap2
gp.Reclassify_sa(in_raster, "VALUE", full_remap, out_raster, "NODATA")


texts: convert REGCAN95 coordinates to OpenStreetMap coordinates (WGS84-Form in lat/long)

texts:
I have some coordinates in datum REGCAN95, and I want to know how to convert them to something usable in OpenStreetMap.

I have found that REGCAN95 uses the EPSG codes as follows:

4082 REGCAN95 for UTM 27
4083 REGCAN95 for UTM 28

For example I have the  coordinate 28BR205306 or 28BR067530, which I want to find in OpenStreetMap or other webmaps. How can I convert them?


texts: Selecting only portion of shapefile using ArcMap?

texts:
I have two layers: One of the Misissippi River Basin  and one of the USA counties within the 48 states.

I'm having no trouble selecting the counties within the Basin, but I also need to select the portions of the counties only partially within the basin. I know how to select the entirety of the counties that are partially within, but I want to select only the portions that are within the basin shapefile so I can calculate the percentage of those counties that are within the MRB.

Edit for more info: Take an example that shows a selected county - the counties layer is red and part of the county selected contains a portion of the Mississippi river basin shapefile (green). Rather than select the whole county when I select by location (I've selected those that intersect the MRB shapefile, those that are wholly within it, and those only partially within it - I still get entire counties) I would like to select only the portion of each county that contains the MRB shapefile shown in green.

Furthermore, I would like to retain the borders of each county - My end goal is to extract county-level data from a raster, but only the portions of counties within the River Basin.


texts:
I want to UnderStand what is cache Tiled Services and how can it improve performance.I am facing a performance issue in my application.I created a cache services and loaded the map without adding feature layer it performed well

  var opLayer1 = new esri.layers.ArcGISTiledMapServiceLayer("http://10.32.9.106:6080/arcgis/rest/services/KX_Mapping/World_Countries_WMS_Cached/MapServer");
    map.addLayer(opLayer1);

but when i added feature layer again performance got hit , how can i imporve performance of this SO actually using cache services dint help or DId i do some mistake please help.

adding feature layer

featureLayer = new FeatureLayer("http://10.32.9.106:6080/arcgis/rest/services/KX_Mapping/World_Countries_VD_Look_v2/MapServer/0",{
            mode: esri.layers.FeatureLayer.MODE_SNAPSHOT,
             outFields: ["*"],
             opacity: 0.75
          });


texts: True long int fields in a File Geodatabase

texts:
I have been wondering about the 'long ints' in FGDBs. They are 32 bit, which is just a regular int in most places. Is it possible to actually have a field be a 64 bit int?

I have some bigint's (64 bit) in a PostGIS database that I am translating into a FGDB with ogr2ogr and they all max out in long ints (2,147,483,647) in the FGDB. I know I can cast them as text and it will get the numbers in there, but that makes it harder to deal with them as numbers. I played with the doubles a little bit, but they start to go into scientific notation once they are too large, which is also not useful.