texts: How to only ask for a input if checkmark has been checked? (arcgis script tool)

texts:
I am currently working on a python script, that i want to convert to a tool.
I have managed to set it up in Arcgis as a tool, but it does not behave exactly as i want it to:

What i want it to do is:
if the first checkmark box is not checked, the next input parameters need to be greyed out (not available for the user to specify), and if the first box is checked, parameters will be enabled, indicated by the purple arrows.
and then, if the first checkmark box is checked, (which then enables the second checkmark box), it should enable the last parameter only if the 2. checkmark is checked. Basically creating a "dependency" hierarchy of some sort.

I have tested out the tool-validation, but with very little luck and was hoping someone had some more experience with customizing a tool like this! right now my validation script looks like this:

def updateParameters(self):
    if not self.params[3].hasBeenValidated:
        self.params[4].enabled = True
        self.params[5].enabled = True
        if not self.params[5].hasBeenValidated:
            self.params[6].enabled = True
        else:
            self.params[6].enabled = False
    else:
        self.params[4].enabled = False
        self.params[5].enabled = False
return

(edit: the above script does not behave properly: when you open the tool, all input parameter fields are NOT greyed out, but when you click on the first checkmark, the last field [what is the offset above floor level?] is greyed out. then if you click the last check mark, both the last checkmark and the field [Height field] is greyed out. and if you uncheck the first checkmark, the field [Height field] and the last checkmark is un-greyed out.   )


texts:
I have a polygon shapefile I need to extract the line segments where the polygon features intersect each other and also the segment where they don't.

Can we automate this process or can scripting like ArcPy be used?


texts: Which affine format to use in Rasterstats zonal_stats?

texts:
I am working with the Rasterstats package to run zonal statistics on two raster datasets. The documentation calls for the following:

import rasterio
from rasterstats import zonal_stats

with rasterio.open('tests/data/slope.tif') as src:
    affine = src.affine
    array = src.read(1)
zs = zonal_stats('tests/data/polygons.shp', array, affine=affine)

However, src.affine is no longer supported in Rasterio. Instead, the documentation recommends using src.transform, which yields something like:

Affine(1.0, 0.0, 670625.0016573452,
       0.0, -1.0, 5282947.795216726)

When I do something like the following, I get an error:

from rasterstats import zonal_stats

dataset = rasterio.open('/path/to/image.tif')
affine = dataset.transform * (0,0) # upper left corner
zs = zonal_stats(segments, ndvi, affine=affine, stats="mean")

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

What affine format does the Rasterstats zonal_stats function require?


texts:
I am new to R for GIS and I am trying to get universal kriging working in Autokrige in R but I am running into problems. Ordinary kriging works with my dataset so I think the problem is to do with the predictor grid which is an ASC file. The data I used is here. https://www.dropbox.com/sh/lkg8qtp2eowmml8/AABWYWegF-FgsbIT1i-5y6VEa?dl=0

library(sp)
xdata <- read.table(file="C:\\Users\\Robert\\Dropbox\\Dplaces.txt",header=TRUE,na.strings = "NA", sep="\t")

p <- xdata
p$g <- p$Asin.gather
coordinates(p) <- ~ x + y

##Set the coordinate system
proj4string(p) <-CRS("+proj=utm +zone=10 +datum=WGS84")

##Load first temperature GeoTiff raser  from worldclim
library(raster)
f <- "C:\\Users\\Robert\\Documents\\QGIS Database\\bcmidbi_10m\\bcmidbi11.tif"
       r <- raster(f)
bcmidbi11  <- aggregate(r, fact=2)  ## By default aggregates using mean, but    see fun=

##Load second temperature GeoTiff raser from worldclim
       f <- "C:\\Users\\Robert\\Documents\\QGIS Database\\bcmidbi_10m\\bcmidbi10.tif"
       r <- raster(f)
bcmidbi10  <- aggregate(r, fact=2)  ## By default aggregates using mean, but see fun=

   ##calculate global effective temperature from rasters which is my predicator variable
x <- bcmidbi11  /100
y <- bcmidbi10  /100
a <-   18*y
b <-   x
c <-   y-x
d <-   c+8
e <-   a-b
MET_map <-   e-d
MET<-read.asciigrid("org_matter.asc")
proj4string(MET) <-CRS("+proj=utm +zone=10 +datum=WGS84")

ra <- aggregate(MET_map, fact=2)  ## By default aggregates using mean, but see fun=
writeRaster(ra, "C:\\Users\\Robert\\Documents\\outfile.asc", format="ascii")

Error in .startAsciiWriting(x, filename, ...) :
x has unequal horizontal and vertical resolutions. Such data cannot be stored in arc-ascii format

When this I run into some conversion problem. I tried to experimentally continue with the Meuse dataset, using the Meuse 's data set dist to continue.

x.range <- as.numeric(c(-169.00, 172.00))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(-53.00, 79.00))  # min/max latitude of the interpolation area
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1), y = seq(from = y.range[1],
to = y.range[2], by = 0.1))  # expand points to grid, from dataframe to spatialpixels

coordinates(grd) = ~x+y        # from dataframe to SpatialGrid

library(automap)
kriging_result = autoKrige(g~dist, p, new_data=grd)

Error in model.frame.default(terms(formula), as(data, "data.frame"), na.action = na.fail): invalid type (S4) for variable 'dist'

Have you any idea what the problem is?

I have tried kriding in other packages and I have had similar problems.


texts: How to avoid corrupt geometries when translating between TopoJSON and GeoJSON?

texts:
I'm trying to display some geojson vector features on a leaflet map.
The size of my GeoJSON is too big (25.4mb of data, 8.47mb gzipped), so I'm trying to use a TopoJSON format to reduce the size without loosing detail.

I converted my geojson data with topojson tool and just basic params:
topojson --id-property okato_code -p name -o adm.topojson adm.json

That significantly reduced the file size (down to 1.54mb of data, 201.3kb gzipped).
But when I translate it back ( using topojson.object() ) to GeoJSON on the client, it looks overgeneralized and like it has some incorrect topologies when zooming in.
I thought there was no generalization happening during translation to/from topoJSON, then why are the geometries different?
How to make them identical? My goal is just to reduce the size of the file without any simplification of geometries.

Here's the map - http://basemap.ru/json/
You have to wait for the big GeoJSON to download, it's included for comparison,
but after it loads, you can see a distinct difference on high zoom levels.
GeoJSON layer is orange and TopoJSON layer is white. I've included the layer switcher to turn topojson and geojson on and off.

The data itself is here:
http://basemap.ru/json/json/adm.geojson
http://basemap.ru/json/json/adm.topojson
The script with translation is http://basemap.ru/json/mapscript.js

If it makes sense, my map is in Asia Lambert Conformal Conic projection (http://spatialreference.org/ref/esri/102012/)


texts: Resolving non-keyword arg after keyword arg error in Python script?

texts:
I tried the following script to get drainage basin of specific points:

import sys, os
import fiona
import networkx as nx
import itertools
from shapely.geometry import Point, LineString
import numpy

sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
print 'import GRASS ok'
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'nl'
mapset = 'nl'
gsetup.init(gisbase, gisdb, location, mapset)

H = nx.Graph()
for line in fiona.open('tc_line.shp'):
    for seg_start, seg_end in itertools.izip(line['geometry']['coordinates'],line['geometry']['coordinates'][1:]):
        H.add_edge(seg_start, seg_end)

edge_node = []
i = 1
for node in H.nodes_iter():
    if H.degree(node) > 2:
        for edge in H.edges(node):
            print edge
            g.run_command('r.water.outlet', overwrite = True, input = 'dra', basin = 'netx' + str(i), coordinates = edge[1][0], edge[1][1])

And the script ran fine until I add g.run_command('r.water.outlet', overwrite = True, input = 'dra', basin = 'netx' + str(i), coordinates = edge[1][0], edge[1][1]), and the error showed up:

File "C:\Users\Heinz\Desktop\netx_test.py", line 29
g.run_command('r.water.outlet', overwrite = True, input = 'dra', basin = 'netx' + str(i), coordinates = edge[1][0], edge[1][1])
SyntaxError: non-keyword arg after keyword arg
[Finished in 0.1s]

So I am wondering that could Fiona work with GRASSGIS? And I want to know how to solve this problem.

I am working under Win10 64bits with python 2.7.12, Fiona 1.7.0, networkx 1.11 and GRASS GIS 7.0.5 installed via OSGeo4W pacakge.


texts: GEE Select images that have certain number of pixels

texts:
I have code that calculate NDVI for Landsat 5,7 and 8. I would like to know how  can I "filter" my image collection so I have only images that have a certain number of pixels. For that I thought maybe I can first count the number of pixels in the image and then use gt but I'm not sure how to put it inside the imageCollection.

To do that I think I need to take my imageCollection and change it into a list and then iterate through this list- for every image to count the numbers of pixels using something like this-

//count the number of total pixels
var countpixels = image1.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: table,
  crs: 'EPSG:4326',
  scale: 30,
  });

and then filter according to the number of pixels.

My problems are-

I have around  250 images in my imageCollection. I don't want to create list and to type by myself all the numbers like this-

var listOfNumbers = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19];

How can I filter image by number of pixels so I don't get very small rasters?

My end goal is to calculate statistic for my polygon for many years and for that I can't have the small rasters or the ones that look like cheese with many holes due to the cloud mask.

Here is one of my codes-

/**
 * Function to mask clouds based on the pixel_qa band of Landsat SR data.
 * @param {ee.Image} image Input Landsat SR image
 * @return {ee.Image} Cloudmasked Landsat image
 */
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2).divide(10000)
  .copyProperties(image, ['system:time_start']);
};

//Create LANDSAT7 dataset

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filterDate('1999-01-01', '2013-04-29')
                  .select('B1','B2','B3','B4','pixel_qa')
                  .filterBounds(geometry)
                  .map(cloudMaskL457);

//RGB visualization

var visParams = {
  bands: ['B3', 'B2', 'B1'],
  min: 0,
  max: 1,
  gamma: 1.4,
};

//clip the dataset according to the geometry
var clippedCol=dataset.map(function(im){
  return im.clip(geometry);
});

// Get the number of images.
var count = dataset.size();
print('Count: ',count);
//print(clippedCol);
print(dataset,'dataset');

//function to calculate NDWI in LANDSAT7
var addNDVI = function(image) {
  var NDVI = image.normalizedDifference(['B4', 'B3'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(NDVI);

};

//NDWI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var NDVIcolor = {
  min: -1,
  max:1,
  palette: ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'],
};
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.mean(), 30));
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.max(), 30));
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.min(), 30));
print(ui.Chart.image.series(withNDVI, geometry, ee.Reducer.stdDev(), 30));


texts: Remove duplicate labels for parallel line features in QGIS 3

texts:
I've been really struggling to get my highways to label neatly. The idea is to have one route number label every 1 km on the route. The only issue is that there are duplicate separated carriageways with the same route number. And then the label gets added twice, for both directions of travel:

What I'm actually looking for is this (one label, placed on either side of the separated carriageway):

I've tried tweaking the variables in the labeling options, but nothing gives me this result. I also tried dissolving based on Route Number (the label field) and then labeling without placing a label on every multipart feature. This works only if there are no lines like my blue line. If the same Route Number repeats somewhere else, like in my example image, then only the one or the other gets labeled.


texts: Compute mean attribute value of point features in same category?

texts:
I've got a point shapefile wih a column with categorical variable and a column with numerical values like the one below

Point | crop     | NDVI

1     | wheat    | 0.9
2     | barley   | 0.7
3     | wheat    | 0.2
4     | barley   | 0.5

How can I aggregate the points with the same crop showing also the mean NDVI (in a new column) like the result reproduced below in QGIS?

Point | crop     | NDVI | meanNDVI

1     | wheat    | 0.9  | 0.55
2     | barley   | 0.7  | 0.6
3     | wheat    | 0.2  | 0.55
4     | barley   | 0.5  | 0.6


texts: What channel initiation threshold value should I use when creating a stream network for drainage density calculation?

texts:
I need to calculate drainage density for two southern AL watersheds (D’Olive and Fly Creek).

What channel initiation threshold value should I use when creating the stream network?

I am following these instructions http://gis4geomorphology.com/watershed/ to create my stream network and will calculate drainage density using this http://gis4geomorphology.com/drainage-density/ calculation.


texts: How can I strip spaces from a moustache template field in CartoDB infowindow?

texts:
I have a custom infoWindow on a cartodb.js map which uses moustache templates to populate the fields in the infoWindow.  In it I have a link to tweet a correction to the map which I populate with the relevant fields. I want to make the name of the {{country}} a hash tag but in order to do this I need to strip the spaces from the field in the case of countries with more than one word e.g.  Democratic Republic of Congo which I would like to appear as #DemocraticRepublicofCongo.  Is there a way to perform javacript operations like .replace(/\s+/g, '') on a moustache template? The infoWindow code follows and the full web page can be found at https://afterfibre.nsrc.org/af_dev.html with full source at https://github.com/stevesong/nsrc-afterfibre/blob/master/af_dev.html

  x

{{operator}}

Phase

{{phase_name}}

Country

{{country}}

Owner

{{orgname}}

Type

{{type}}

In Service

{{go_live}}

Source

map

Contribute

Tweet a correction!


texts: No hillshade effect when no-data area is on screen

texts:
I'm facing this weird problem when hillshade effect partially brakes if a no-data area appears on screen.
If I move the area to hide it off screen, the hillshade effect works.
I tried changing symbology, colors,classes etc without change.


texts: Unicode problems while writing selected features to a new shapefile with Fiona

texts:
I have a shapefile with a large number of geographic features (marine protected areas).  I'm trying to read through all of the areas in the shapefile, select those that fit a certain criterion, and write a new shapefile with only the selected areas (areas with reefs).

Many of the areas have names with non-latin characters (i.e. in Swedish or Croatian).  I believe I have properly encoded everything in utf-8.

import shapefile
import fiona
from shapely.geometry import Point, shape
import pandas as pd

#read in habitat code list
hab_data = pd.read_csv(habitat_data.csv)

#read in shapefile
with fiona.open("all_protected_areas.shp") as source:
    polygons = [pol for pol in source]
    with fiona.open('reef_habitats.shp', 'w', encoding='utf-8', **source.meta) as outfile:

        #loop over polygons
        for i, poly in enumerate(polygons):

            #match the site code with the habitat code from csv
            sitecode = poly['properties']['SITECODE']
            habitatcodes = hab_data.loc[hab_data.SITECODE==sitecode]['HABITATCODE']

            #choose only reefs (1170)
            if '1170' in habitatcodes.values:

                outfile.write(poly)

The code works fine until I run into a polygon with a record containing non-latin characters.  The first one is: ('SITENAME', b'Obala izme\x1au rta \x8ailo i Vodoto\x1a').  The shapefile fails to write on the final line with:

ValueError: Invalid field type <class 'bytes'>

I'm using Python 3.  From my own investigations I've seen that there are some tricky differences in the way that Python 3 and 2.7 handle unicode.  There also seems to be some weirdness in writing unicode with fiona.

By using Python 2.7 and various combinations of .encode() and .decode(), I've managed to get errors like UnicodeEncodeError and UnicodeDecodeError but no solution.

I'm open to using a different package if there is an easier way to do this feature selection.


texts: Overwriting existing field names during MapInfo to ArcGIS conversion?

texts:
I have a large number of MapInfo .TAB files which all need converting to ESRI Shapefiles. As far as I know, the Universal Translator tool within MapInfo is the simplest method to carry out this task. However, I'm also aware that this conversion will truncate all field names to a maximum of 10 characters. Unfortunately, I cannot bypass this limitation by using an alternative file type; ESRI Shapefile is our only option.

If I were to obtain a list of shortened field names in an Excel file (for example, wouldn't necessarily have to be Excel), would it be possible to import these (into either MapInfo or ArcMap) and overwrite the existing field names?

Bonus points for if this can be implemented into the conversion process. I'm trying to limit the amount of time I spend converting these files.

When using Universal Translator, I would select multiple .TAB files to convert at a time, so that I can let it run in the background while I work on something else; minimising the time actually spent on the task. It would be great if I could still do this, as opposed to selecting each .TAB file individually, and then selecting the respective file containing the list of field names.


texts:
I wondering if there is a way I can take my field data collected on a Trimble R10 RTK GPS and import it into ArcMap for further processing, I have done some research and I know of Trimble GPS Pathfinder Software, is this the only way?

If there is a way to bring it into Quantum, that would also be alternatively helpful.

Thanks Always


texts:
I am trying to convert a python code from field calculator to arcpy (pycharm).
All i am wanting to do is convert the below python code into a way python in pycharm will read it.
!Owner_Ab!.title()

This is my attempt however it does not work.
feature_class = CAD_WA
FukMoO = 'Owner_Ab'

with arcpy.da.UpdateCursor(feature_class, FukMoO) as cursor:
    for row in cursor:
        row[0] = row[0].title
        cursor.updateRow(row)


texts: Label legend interval with number of elements in ArcMap?

texts:
I have made a series of choropleth maps in ArcMap 10.
When I classify the interval the bottom of the dialogue tells how many elements are in the class.

Is it possible to have this number automatically places beside my class labels?


texts: Pyscripter crashes with EOFError timed out

texts:
I am not sure whether someone ran into similar problem.
I've been using PyScripter with Arcpy for the past two years and I have been happy with it most of the time.

Yesterday I ran into a kind of weird problem. I wrote a pretty straight forward and simple script. Basically it reads through text file, gather information and then based on the information gathered running first a attribute search on a given layer (arcpy.SelectLayerByAttribute), then union all the geometry that have been selected, then insert to a newly created featureclass.

I have to do this a lot of times for the whole US nation.

The codes runs fine in Pyscripter then it crashes at about 18505 times of selecting layer, with an error message saying 'EOFError Timed Out', then pyscripter freezes. I cannot shut it down but clicking 'x', it kept telling me it's in debug mode. I have to shut it down using task manager.

I tried to re-run at the data from where it stopped and this time it runs fine all the way till it hits ~18505 again. I tested about 10 times on different data and it always crashes at ~18505 times doing SelectLayer.... I've been stuck there for a whole day.

Then finally I decided to run the script in command line and it runs all the way till all the data are done.

I wonder if there's a limit set somewhere in pyscripter? I like the ide and want to continue to use it. But this is very frustrated. I don't want to run it to this and waste a day again.


texts: Retrieve data from multiple selected cells in QTableWidget for a QGIS Plugin?

texts:
I am creating a plugin which shows the list of location of vector/raster layers. I am using QTableWidget to show the data. When I select one row/cell and click a QPushButton it opens the specific vector/ratser layers in the QGIS. I use the following code to retrieve the data from the selected cell:

row = self.tableWidget.currentRow()
path = (self.tableWidget.item(row,3).text())

How do I select multiple cells at the same time to open multiple files?


texts: Change LineString type of Osm2Po edge and vertex file

texts: Is it possible to change the linestring type from 2 to 3 dimensions in osm2po, and if so, how? Currently the generated table column looks like this: geom_way geometry(LineString,4326), but I would much rather like to have geom_way geometry(LineStringZ,4326),


texts: Can 2D polygon be used as hard breakline for terrain in ArcGIS Desktop?

texts:
I've been working with TINs for years and have only ever been satisfied when using 3D polylines for breaklines.

I am currently transitioning to using terrains and after reading the documentation in ArcGIS, it seemed that I ought to be able to apply a polygon as a breakline.

Intrigued by the possibility, I decided to try it with a LiDAR data set I received from a project partner along with a 2D polygon feature class of wetted stream channel planform.

I specified the wetted channel polygon as a hard breakline, built the terrain, and ended up with two very deep (0 elevation) troughs along the extent of each bank.

See below:

Based on my reading the software should create nodes where the polygon perimeter intersects triangle edges and interpolate z-values along those edges, but that does not seem to be the case.  Is this actually possible?

I am using ArcGIS 10.0 SP4 with an ArcInfo license and 3D Analyst and Spatial Analyst extensions.


texts:
I'm having trouble registering a table with a Geodatabase in an SDE, the option is greyed out. This is a weekly task for me and I have never had this happen before.

I receive a new table weekly from our accounts department and have to check if there is any new entities in it to add to my parcel database. When I right click on the table and go to mange the option to register with database is greyed out.

I checked for locks on the table and there are none. I previewed the table and it appears to be the same as every week

I am the only one working off the database at the moment, has anyone run into similar problems before? Any suggestions on things to try?


texts: Integrating shapefile in spatialite DB with spatialite gui

texts:
For a project, i need a to display datas of a shapefile in an Android app (Geopaparazzi). I'm doing it by integrating the shapefile into a spatialite DB with spatialite gui 1.5 (1.6 having issues with Geopap).

It usually works fine, but with my new shapefiles in wgs 84, it gives me an error saying it can't open the shapefile.shx file, while it's in the directory with the .shp .

It still works with other files, so i don't really understand and i can't reach the person who created the new shapefiles right now.

Any ideas ?


texts: PostGIS - Merging linestrings into multilinestrings in a particular order

texts:
I have a set of linestrings, which I want to group into a single multilinestring, and in a specific order. However, when I use ST_Collect, it merges them out of order, like this:

Geometries as linestrings:

Order at which the multilinestring is being generated:

I know this "order" is only in my head really, and the function will use some parameter to determine the order. A quick research has given me that ST_Collect does its magic ordering by the row id, but for reasons I won't bother going through here I can't reorder my ids in the "aggregating order".

I do, however, have a separate table with the correct order, that's referencing the geom ids, so a simple join would suffice to put them in order, but ST_Collect is ignoring that. Is there a way to aggregate linestrings into a multilinestring using the ordering from this second table?

(the drawings are merely illustrative, but depict a real case, that is, some of the lines are not connected, and must remaind so. Therefore, I can't use ST_LineMerge to "force" a correct order out of it)


texts: How to get feature for a given bbox from shapefile by ogr

texts:
I have some features point,line,polygon with shapefile format, now I want to get features for a given boundbox, is it possible to remove the features out of the bbox?

And for a polygon, I think it is necessary to close it.

I wonder if this is possible?


texts: How to maintain polygon boundary location during overlay in ArcGIS

texts:
I am working with ArcGIS Desktop 10.1 (Advanced license).

I am trying to overlay two polygon layers (buildings and landuse) via an Update overlay such that all building polygons are retained and the remaining areas have the landuse classification.

I want to collapse any areas less than 5 metres in width, and I do not want the building polygon boundaries to shift.

i.e., If there is a 3m gap between the edge of a building and a landuse boundary, I want the landuse boundary to collapse to the building edge.

Does anyone know how to do this?

If I use the Update tool with a 5m tolerance, all boundaries (landuse and buildings) are shifted.

Thanks for any suggestions.


texts: eval() statement in QGIS data defined override expression editor gives unexpected results

texts:
In QGIS I've set up some project variables called rv_farbe1, rv_farbe2 and so on. The contain color values in html notation like '#ff0080'.

Now in rule base labels I would like to setup data defined override by selecting the text color based on an attribute value sektor_lfd_nr. The idea is to concatenate a string @rv_farbe and the sektor_lfd_nr and evaluate the result in order to access the variable value like

eval('@rv_farbe' ||  "sektor_lfd_nr")

As far as I understand, in expression editor preview I should get the color value, but preview (german: 'Ausgabevoransicht') displays nothing:

Furthermore, the text with data defined override set for it's color does not display.

Both using the variable directly (@rv_farbe1) as well as using the html style color code directly ('#ff00bb') works, so do I miss something about eval() statement, or what might be wrong here?


texts: Is it possible to have transition between switching base layers

texts: I use OpenLayers to display several maps. While transitioning between Google Map layers, the new layer is superimposed onto the old one. But for OpenStreetMaps, there is this empty void in between transitions. Is there a setting for making individual tiles of the old layer stay until the tiles of the new layer are available?


texts: What Spatial Reference System do I store Google Map's Lat/Lng in

texts:
I am new to PostGIS, Spatial Reference Systems and projections and want to store lat/lng coordinates retrieved from Google Maps and its services like the Geocoder. I'm working with data in America and the main calculations are for nearest neighbors. Results will be plotted back on a web map (eg: via Google maps API or leaflet)

Question: What spatial reference system is Google map using? I figured out so far that its the WGS 84 datum with the Mercator projection in the geographic coordinate system. Should I store the location data as it is, or transform to US National Atlas Equal Area EPSG:2163?

Since I will be calculating distances and finding nearest neighbors, if I guess correctly that doing nearest neighbor search requires transforming the entire table to EPSG:2163, then will the solution be to store the data in both Mercator and EPGS:2163?


texts: How to disable overlay check box after a certain zoom level?

texts: How do I disable a particular overlay checkbox in the Openlayers Layerswitcher?


texts: Automating production of GPS-compatible map tiles from georeferenced raster maps

texts:
I have a set of TIFF maps that I want to load onto Garmin GPS 62s units. The maximum size for each tile on the GPS is 3Mb and less than 1 million pixels. I manually produced one tile, following this workflow:

Convert whole TIFF map to JPG using IrfanView
Manually crop the image to be less than 1 million pixels (i.e. 1024X1024) and compress to be less than 3Mb.
Add JPG to QGIS and georeference using the georeferencer tool.
Add same JPG image to Google Earth as an Image Overlay.
Use the Extent information from the georeferenced JPG in QGIS to fill out
the extent boxes in Google Earth.
Save as KMZ file.
Add KMZ file to Garmin Basecamp.
Add KMZ to GPS unit with Basecamp.

This is obviously a slow and fiddly process which I am reluctant to repeat 100 times to produce all the tiles needed.

Is there a way to automate sections of this?

Specifically I imagine the follow are possible, but don't know of any tools to do it:

Automatically divide one georeferenced raster into multiple tiles of given dimensions, retaining geocoded information, in QGIS.
Export multiple tiles into a KMZ file in QGIS

This would make the process a lot quicker. The closest tools I have found is the GarminCustomMaps plugin for QGIS, but this just takes a JPG image of your current map canvas - great for using with multiple vectors etc, but when zoomed out with large raster files, I just get a very pixelated image which is useless for reading on a GPS. Zooming in more to increase the resolution negates the point of it, as you'd have to repeat it multiple times.


texts: Longitude values in GeoJSON out of range -180 to 180

texts:
I've extracted a geojson data file from QGIS for a transportation layer.
Now in the file, I can see multiple entries like the following:

        {
            "geometry": {
                "coordinates": [
                    [
                        [
                            482745.218106,
                            6656907.959327
                        ],
                        [
                            482785.825278,
                            6656836.299613
                        ],
                        [
                            482820.460806,
                            6656815.996028
                        ],
                        [
                            482967.36322,
                            6656599.822557
                        ]
                    ]
                ],
                "type": "MultiLineString"
            },
            "properties": {
                "_col": 4194,
                "_row": 5456,
                "_zoom": 13,
                "brunnel": null,
                "class": "service",
                "oneway": 0,
                "ramp": 0,
                "service": null
            },
            "type": "Feature"
        }

I guess in the coordinates are the latitudes and longitudes. Normally, the range of values should be between -180 and 180, but here it isn't. Why?


texts: Arcgis online and Collector: Controling shared maps and data

texts:
I am using ArcGIS Online and Collector to conduct field surveys. I have a 'Web Map' (call this 'Web Map 1') that is 'served' to Collector which has a single point feature class in it. This point feature class collecting one item/feature that has a time date stamp to it. This all works fine.

In order to view the data I have also created a 'Web Mapping Application' to view the data via web browser. As the data I am collecting has time/date stamp I need to show records collected in the last 24 hours, 48 hours and all records. To do this I created 'View' layers of my original 'Hosted Layer' and applied a filter. I have added these 'View' layers to a new 'Web Map' (call this 'Web Map 2' as I want to keep my original 'Web Map' simple and just for collecting the data in Collector. The 'Web Mapping Application' use 'Web Map 2' for it's data. This works fine.

The issue I have is I only want staff to see 'Web Map 1' in Collector and not 'Web Map 2'

Is there anyway to only share (or setting) a 'Web Map' so that it is not visible in Collector but still visible to users as 'Web Mapping Application'?

Alternatively is there another way/workflow to do this?


texts: NULL data affecting weighted mean calculations in QGIS?

texts:
I am using QGIS. I've been calculating weighted means based off a database. The database was created by combining two different sets of Data together. The Attribute Table then shows Null values for two reasons: (1) mismatch while joining the two databases together and (2)missing values within the databases themselves.

I have tried to find a logic in the effects Nulls values have on the weighted means but have failed to come to a solid conclusion. They definitely do not act as a ZERO value.

I came across this post in reference to my search:
Null values in Geographically Weighted Regression Analysis

It helped explain a lot of things, however this sentence "A null means something entirely different than zero" is the answer to which I am looking for.

What does a Null value mean in weighted mean centers?


texts: Deleting shapefile and related extension files from folder when not in use

texts:
I have to delete two shapefile layers from my folder. I use these shapefiles to run processing Algorithms after which i want to delete them. The code I use use is as given below.

for filename in os.listdir(os.path.dirname(inter)):
    if (nameNewdel+'.shp' in filename) or (nameNewdel+'.shx' in filename) or (nameNewdel+'.dbf' in filename) or (nameNewdel+'.prj' in filename) or (nameNewdel+'.qpj' in filename) or (nameNewdel+'.cpg' in filename):
        os.remove(os.path.join(os.path.dirname(inter),filename))

where inter is path to the file and nameNewdel is filename without extension.

I am able to delete one shapefile with this code but the same code gives WIN 32 error when trying to delete another shapefile. The shapefile is not in use. I am able to delete the shapefile with this code in the python console of QGIS but not able to delete it using pyqgis. When I introduce try and except for os.remove in the above code it deletes all the files of the second shapefile except .shp and .dbf.

I have referred to Releasing PyQGIS file locks?, Deleting Shapefiles using QGIS? and WinError 32 error when trying to delete shapefiles with QGIS python
Accordingly I have added and removed the layer from QGIS and still not able to delete .shp and .dbf file of the shapefile. I also set the shapelayer to none and tried deleting it using QgsVectorFileWriter.deleteShapefile.

for filename in os.listdir(os.path.dirname(inter)):
    if (nameNewdel+'.shp' in filename):
        dellayer = QgsVectorLayer(inter, nameNewdel,"ogr")
        dellayer = None
        QgsVectorFileWriter.deleteShapeFile(filename)

.shp and .dbf files are still not deleted. Is there any alternative to work around this?


texts: What is the scope of DIVA-GIS geodata that are available for download?

texts:
In the 'Free Spatial Data' download section it says the geodata is "free" - is there a concrete license? Is it public domain, creative commons? Can it be used for commercial purposes?

When I download it, there is an attached license text file that says: "These data were extracted from the GADM database (www.gadm.org), version 2.5, July 2015. They can be used for non-commercial purposes only.  It is not allowed to redistribute these data, or use them for commercial purposes, without prior consent."


texts: ArcGIS Silverlight API 2.4 support for Microsoft Silverlight 5.0?

texts: Please let me know whether ArcGIS Silverlight API 2.4 version officially supports or CERTIFIES Microsoft Silverlight 5.0.?


texts: Importing Bing Maps into application using ArcGIS Engine?

texts:
I'd like to create an application using ArcGIS Engine 10.

I want to import Bing Maps into my application as base map. I surf the internet but couldn't find any tools or API to import these maps to my application.

Is there any way to import the maps to application using ArcGIS Engine?


texts: How do I implement PyQGIS multithreading to update a progress bar

texts:
I am developing a PyQGIS plugin using python. I facing a problem during the implementation of threading in my plugin. I have written a Python script which contains a class wrapping all geo-processing methods in use. The problem I am facing is processing the methods in threads. To give a better picture of my problem, I have a UI which calls the method calcNDVI() as:

self.lst.calcNDVI(self.RedBand, self.NIRBand, None, self.outputRaster, self.sensorType, self.dataType, self.addToQGIS)

This function call, calls a calcNDVI() method which does all the magic and adds the result to a QGIS project. The method call is also associated with passing some arguments.

I have gone through many sources but I can't get to understand how I can implement my algorithms in a threaded way. The code I have come through implements everything in the run method which makes it difficult for me as I have more than one method to be processed in my class, and I would also like the user to be able to choose which method to be run through the UI.

I am processing Landsat datasets which take too much time to process.


texts: Export data from "Identify Results" table to Attribute Table

texts:
I am working on a street network in QGIS. I needed to increase the resolution of the street segmentation in an area of my network, so I copied and paste this area from another layer.

I selected all features from this network:

And paste it to my original network:

The street segmentation is now fine (I got the resolution I wanted in the area I was interested to). The problem is that in the resulting layer, in the areas where the network I pasted overlaps with my original one, I completely lost all of the data in the attribute table.

As you can see, I selected a random street segment in the overlapping area and the corresponding columns in the attribute table are all NULL.

However, the good news is that the data I wanted to keep (street ID and traffic data TJM_2010) are still present in the Identify Results:

Do you know a way to export this data into the attribute table? Since the street segments that now lack information in the attribute table are quite a lot, it would be great if the method was automatic.


texts: Why do these coordinates not display correctly (or intersect correctly) with this boundary file?

texts:
This is a follow up to my previous question. I wrote code that imports the boundaries of a few US states, and, for some arbitrary point within those boundaries, buffers a circle around that point and then draws both on a map.\

Even though the latitude and longitude specified fall squarely within the boundaries, they don't appear this way on the map. Am I using the wrong projection? I chose the Albers projection per this question. I guessed that this may be the problem because the axes aren't in longitude and latitude either; they seem arbitrary.

This is the output, which is incorrect. The blue point/buffer should be in the northwest quadrant of the red map.

This is the code:

library(sf)
library(tidyverse)
library(USAboundaries)

albers_proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

# download county boundaries
counties <- us_counties(map_date = "2000-01-01", resolution = "high",
                        states = c("Utah", "New Mexico", "Colorado", "Arizona"))
counties <- counties %>% st_transform(albers_proj)

# example point
lon <- -111.101099
lat <- 40.477140
elev <- 150
tower <- st_point(x = c(lon, lat, elev), dim = "XYZ")
tower <- tower %>% st_sfc(crs = albers_proj)

# draw the radio horizon around the tower
scalar <- 3.57
radius <- scalar * sqrt(tower[[1]][3]) * 1000
tower_buffer <- st_buffer(tower, dist = radius)

xlim <- c(-2e06, 1e06)
ylim <- c(-6e05, 6e05)

png("map.png", width = 1100, height = 850, units = "px")
plot(st_geometry(counties), col = adjustcolor("maroon", alpha=0.5), xlim=xlim, ylim=ylim, axes=TRUE)
plot(st_geometry(tower_buffer), col = adjustcolor("blue", alpha=0.5), add=TRUE)
dev.off()


texts: ArcGIS Model Builder Copy Features tool fails with 10.2 fGDB as input

texts:
So far I cannot detect any difference between the two input feature classes other than that anything in pre 10.2 copy features will work. If I set the input to be a FC from a 10.2 file GDB, I get the 'ol 000210 error. Projections aren't an issue, both inputs reside in the same level of the file structure, no differences in OS permissions. I've created a fresh fGDB for testing, same behavior. I can even have both FC versions in the same MXD and just switching the input in the model decides success or failure.

To be fair, when i say non-10.2 GDB, my older ones are almost all 9.3. I'll try to test with an intermediate version as well. I have tested with different output locations and no problems there. Its odd because the error raised is related to output...not input. I am currently writing to %scratchFolder%

And to honor the Stack Exchange protocol, here is my question:

Has anyone seen similar if not identical behavior from this? Again, it appears to only affect 10.2 fGDB using the copy features tool.


texts: Pandas to GeoJSON (Multiples points + features) with Python

texts:
I have this working code, data2kml works ok, but data2geojson is incomplete.
How can I add the elevation, name and the description in the geojson? I can't understand how FeatureCollection works.

Obs: I can't use geopandas because of some dependency conflicts.

import json
import simplekml
import geojson
import pandas as pd

def data2kml(df):
    kml = simplekml.Kml()
    df.apply(lambda X: kml.newpoint(
            name=X["name"],
            description=unicode(X["description"].decode('utf8')),
            coords=[( X["long"], X["lat"], X["elev"])]
                                    )
            , axis=1)
    kml.save(path='map.kml')

def data2geojson(df):
    points = []
    df.apply(lambda X: points.append( (float(X["long"]),
                                       float(X["lat"]))
                                    )
            , axis=1)
    with open('map.geojson', 'w') as fp:
        geojson.dump(geojson.MultiPoint(points), fp, sort_keys=True)

col = ['lat','long','elev','name','description']
data = [[-29.9953,-70.5867,760,'A','Place a'],
        [-30.1217,-70.4933,1250,'B','Place b'],
        [-30.0953,-70.5008,1185,'C','Place c']]

df = pd.DataFrame(data, columns=col)
data2kml(df)
data2geojson(df)


texts: Disable map panning using the arrow keys in the keyboard

texts:
I have a Leaflet map that has map.dragging.disable(); which works for mouse dragging. But when the map container is focused using the Tab key, pressing the arrow keys on the keyboard will still pan the map.

What I want is to completely disable map panning whether using the mouse, touchscreen, or keyboard, so that the map appears like a static image. Is there any solution for this?


texts: Convert a data frame in R to a polygon for viewing in Google Earth

texts:
I have a dataframe in R which contains numerous columns of lat long coordinates. Each row of the data frame contains 9 x lat/long coor points.

I wish to plot each of these rows in Google Earth using the plotKML function or any other R package for that matter since I can only use R V2.14.0!!

Ideally I would like to continue using v3.x.x but plotKML doesn't work then so if there's an alternative R->KML package please shout.

I've read up about the differences between polygon and polygons here but I am having difficulty converting my dataframe to a spatialpolygon in order to plot the data in Google Earth. I suspect it's because I have multiple columns of lat long coordinates and plotKML doesn't know (or I haven't specified correctly) what to do with them.

Can anyone point me in the right direction about how to go about this?

I have a sample of my data here in dput() format if it's of any benefit.


texts: Google earth plug in loading old ground overlay

texts:
Im having a problem were i have updated Groundoverlay images and the links, but the plug in is either showing the old overlays or just greyed out, when viewed on my PC

I know the coding is ok as is i view it on my laptop everything is correct so is it a problem with the plug in cacheing old data or do i need to make use of the fetchkml code?


texts: Is there a way to check and repair a GeoPackage GPkg file?

texts:
I have downloaded imagery in GeoPackage format from My DigitalGlobe EV which is used offline on a non-internet connected system. However the GPkg seems to be non-standard.

I'm trying to load the GPkg with GeoServer with the GeoPackage extension. However, when I add a GPkg data store I get a Java exception, which I've filed a bug on.

I've tried using GDAL_translate to convert it to a GeoTIFF, and even convert it from GPkg to a new GPkg, but both cause gdal_translate to hang with 100% CPU usage. GDAL also complains about an invalid 'application_id' in the GPkg file.

Is there a script to repair or rebuild a raster GeoPackage file? Before I contact DigitalGlobe, is there anything else I can try?


texts: Getting area of multipart features using SearchCursor() of ArcPy?

texts:
I am tasked with programmatically getting the area of each individual part of a multi-part feature. I have a script already to read each individual part and print its XY coordinates and it works fine but for some reason can't seem to modify it to also print the areas too.

I have tried inserting the SHAPE@AREA token but receive errors that the "Float object is not iterable". I'm sure it's something simple I am just missing but I can't seem to find an example of this for multipart features. I can write a simple code that gives me total area but not the area of each part. Any suggestions? Here is my script before I include the SHAPE@AREA

import arcpy
from arcpy import env
env.workspace = "C:/Data/Exercise08"
fc = "Hawaii.shp"
cursor = arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"])
for row in cursor:
    print("Feature {0}: ".format(row[0]))
    partnum = 0
    for part in row[1]:
        print("Part {0}: ".format(partnum))
        for point in part:
            print("{0}, {1}".format(point.X, point.Y))
        partnum += 1


texts:
I'm fairly new at this, so maybe this is a simple issue.

I'm trying to make a map that shows a grid of 50 meter by 50 meter boxes.  I've done the Data Frame  Properties for a measured grid and made sure that the coordinate systems match, but no grid appears.


texts:
What is the meaning of "10" in this code below?
I've been searching for it for a long time and I have no idea what it's actually about.


texts: Leaflet GridLayer Vector Slicer with GeoDjango GeoJSON DOM

texts:
Overview

I am trying to use Leaflet's GridLayer Vector Slicer to split up a geojson with complex geometries being passed to my webmap via GeoDjango. I can't seem to get the tiled GeoJSON to display, however.

Documentation for VectorGrid is here: https://github.com/Leaflet/Leaflet.VectorGrid

I can tile the same vector layer from a static framework when I assign the geojson text to a JS variable in a separate geojson.js file, then load that script and tile the vector using VectorGrid.

Static framework that works:

For example, from geojson.js:

var soil_ph = {"type": "FeatureCollection",... }

Then, from index.html:

<script type="text/javascript" src="./geojson.js"></script>

<script type="text/javascript">

        var vectorGrid = L.vectorGrid.slicer( soil_ph, {
            rendererFactory: L.svg.tile,
                }).addTo(map);
...
</script>

This works fine and the display of the vector is snappy.

What I need, but doesn't work:

When I try to call a serialized GeoJSON view in GeoDjango, however, it doesn't seem to work.

From views.py:

def soil_ph_view(request):
    soil_ph_json = serialize('geojson', soil_ph.objects.all(), geometry_field="geom")
    return HttpResponse(soil_ph_json, content_type='json')

From urls.py:

urlpatterns = [
    url(r'^soil_ph/', views.soil_ph_view, name='soil_ph'),
    ...,
]

From index.html:

var soil_vector = $.ajax({
    url: "{% url 'soil_ph' %}"
});

var soil_ph = L.vectorGrid.slicer(soil_vector2, {
    rendererFactory: L.svg.tile,
});

I have also tried the following to pull the geojson data from the Django view:

var soil_vector = new L.GeoJSON.AJAX("{% url 'soil_ph' %}");

But no luck in getting it to tile and display.

How do I resolve this?


texts: Georeference plugin Qgis 3.x not working

texts:
I cannot open the Georeferencer in both Qgis 3.0.2 Girona and 3.2.3 Bonn.
The plugin worked fine in both 3.x versions until today. It still work in version 2.18.24.

It doesn't seem to open any windows at all, so I don't think it is an issue of windows being in another screen as in some other case.
I also had a look at solution suggested here: http://gis.stackexchange.com/questions/83334/georeferencer-plugin-will-not-load-in-qgis , but they do not apply to my case since I do not have a "Plugin-GeoReferencer" folder in my registry for versions 3.x.

Reinstalling did not solved the issue for qgis 3.2.3

any ideas of how to fix the problem?

+++ UPDATE +++

Same problem with the just installed 3.4 MAdeira version


texts:
How can I convert polyline to polygon when the command MAPPOLYLINETOPOLYGON doesn't work at all?

I have the map, which is based on the raster PDF. The most reasonable option is redrawing these shapes and make them as polygons for further georeferencing and export as the .kml files.

Unfortunately, I am stuck with converting the closed polyline as a polygon.

https://knowledge.autodesk.com/support/autocad-map-3d/learn-explore/caas/CloudHelp/cloudhelp/2019/ENU/MAP3D-Use/files/GUID-6129F222-C2C8-4340-AF13-51803B178EAB-htm.html

https://knowledge.autodesk.com/support/autocad-map-3d/learn-explore/caas/CloudHelp/cloudhelp/2015/ENU/MAP3D-Use/files/GUID-C782FDA0-F05F-42EA-B93D-338B2BF062F8-htm.html#GUID-C782FDA0-F05F-42EA-B93D-338B2BF062F8

The MAPPOLYLINETOPOLYGON option is simply not working, because it's unknown in AutoCAD 2020.
I found one of the reason here:

https://forums.autodesk.com/t5/autocad-forum/mappolylinetopolygon-not-working-in-autocad-2020/td-p/9105292

https://knowledge.autodesk.com/support/autocad/troubleshooting/caas/sfdcarticles/sfdcarticles/AutoCAD-Error-Unknown-command-when-entering-commands.html

but there is no solution for 2D option, which I am dealing with.

I marked my polyline as close, but it didn't help.
https://www.cadtutor.net/forum/topic/41524-how-to-convert-polylines-to-polygons/

Does anyone have any suggestions?


texts: Can Elevation Profile data in Arc's Flex Viewer 3.0 be exported?

texts:
a related thread here: google earth profile used in mashup or flex

Would like to export the elevation profile data (csv, txt, etc.) from Arc's Flex Viewer 3.0 as it creates one:
Also, if there is a way to do this export using Google Earth's Elevation Profile tool that would be helpful too.
I know GPS Visualizer can do this with kml/gpx files but would like to have the ability to extract the data directly from Flex.


texts: Exporting a WMS as .TIF fails to grab full extent

texts:
Similar question to this post, however, I am using ArcGIS Pro 2.0.0 and have connected to a WMS.

I've zoomed to the region I wish to export, and when I go to Share > Export > Map, it successfully exports a GeoTIFF, however it looks like a 8.5x11 style rather than the extent of my current map. This is a problem because the region of data that I want  is wider than this extent.

I cannot seem to find a method to change the export extent. I've got to Map > Properties > Extent and attempt to change the settings but the never apply:

Custom Extent: Current Visible Area
Check: Clip layers to extent

Map reloads, then I inspect the properties and they have reverted.

How can I export a WMS as a .TIF with the extent of my current zoom?


texts:
I have a python script tool that I have published to a geoprocessing service. Every time I run the service it fails. Here's what I'm getting for errors:

The python code in the tool is below:

import arcpy, arcgisscripting, sys, os
MapIndexFC = r'C:\Users\p04729\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Testing.sde\MapSets_D04_ALL_HFs_MapIndex_WGS84'
MapIndexLyr = r'\\sv04esri\Records\GIS\MapSheetImages\MapIndexSelected.lyr'
ImageOutput = r'\\sv04esri\Records\GIS\MapSheetImages\Testing.gdb\Image.lyr'
MI_FC_field = 'filename'
out_buffer_FC = r'C:\Users\p04729\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Testing.sde\Buffer'
MapIndexSelect = r'C:\Users\p04729\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\Testing.sde\Selected_Features'

arcpy.env.overwriteOutput = True

Input = arcpy.GetParameterAsText(0)
Fset = arcpy.FeatureSet()
Fset.load(Input)

arcpy.Buffer_analysis(Fset,out_buffer_FC,'5 FEET','FULL','ROUND','ALL')

arcpy.MakeFeatureLayer_management(MapIndexFC,MapIndexLyr)

arcpy.AddMessage('selecting features in Map Index layer...')
arcpy.SelectLayerByLocation_management(MapIndexLyr,'INTERSECT',out_buffer_FC)

matchcount = int(arcpy.GetCount_management(MapIndexLyr).getOutput(0))
if matchcount == 0:
    print('no Map Index features were selected...')
    arcpy.AddError('You must select a feature from the Map Index layer, exiting script')
    sys.exit()
else:
    arcpy.CopyFeatures_management(MapIndexLyr, MapIndexSelect)

    for row in arcpy.da.SearchCursor(MapIndexSelect,MI_FC_field):
        selected_val = row[0]
        img_name = selected_val.split('.')[0]
        arcpy.AddMessage('iterating over selected Map Index features. ..')
        arcpy.AddMessage('looking for image...')
        for dirpath, dirnames, filenames in os.walk(r'\\sv04esri\Records\GIS\MapSheetImages'):
            for file in filenames:
                if file.endswith('.tif') or file.endswith('.sid'):
                    if img_name in file:
                        Out = arcpy.MakeRasterLayer_management(dirpath+'\\'+file)
                        arcpy.SetParameter(1,Out)
         del row, selected_val, img_name

The error is telling me that my Buffer output does not exist or is not supported. I know it does exist, and I'm unsure why it would not be supported?


texts: Does USGS 3DEP data use a hash or checksum?

texts:
I'm downloading some USGS 3DEP elevation data, e.g., the file at the following URL:

https://www.sciencebase.gov/catalog/item/563afea8e4b0d6133fe7357b

The "ID" of this file is given in both the URL and the metadata as the 24-length string 563afea8e4b0d6133fe7357b.

My question is, does this ID correspond to a hash function or checksum applied to the file? I'd like to verify the downloaded file's integrity by comparing some checksum to this ID. Both MD5 and SHA1 produce much longer hashes than 24 characters, and I'm not sure what else to try.


texts:
I am struggling to understand how to render vector tiles in Openlayers 3. this is probably due to me not understanding Z/X/Y tiles.

I have created an mbtile source and am serving it using tileserver-gl-light from a geojson data source using the EPSG:4326 projection.

tileserver-gl-light tiles.mbtiles

Works fine when I view them at (tileserver-gl-client):-

https://server:8080

And I get a tile downloaded with:-

http://server:8080/data/tiles/0/0/0.pbf

How can I now display this layer in openlayers? I have looked at the example here but there is no mention of projection. I tried the following code

var map = new ol.Map({
    layers: [
      new ol.layer.VectorTile({
        source: new ol.source.VectorTile({
            format: new ol.format.MVT(),
            tileGrid: ol.tilegrid.createXYZ({maxZoom: 22}),
            url: 'http://SERVER:8080/data/tiles/{z}/{x}/{y}.pbf'
        })
      })
    ],
    target: 'map',
    view: new ol.View({
      center: [125,83],
      zoom: 8
    })
  });

But have no idea how to centre it to the area I want  which is -4.2376,52.0528 Zoom level 8 in EPSG:4326

When I zoom or pan I get errors like:-

GET http://SERVER:8080/data/sn_all/8/128/123.pbf 404 (Not Found)
GET http://80.229.34.247:8080/data/sn_all/8/129/123.pbf 404 (Not Found)

Which makes me think my projection is not right?


texts: How to correct bad georeference data in a raster?

texts:
I'm doing a bit of vegetation mapping and
 I'm relatively new to GIS.

My issue is that the aerial photographs I have are georeferenced with world files (.jpgw) but they are ever-so-slightly off. The images are referenced in WGS 84 Pseudo Mercator (EPSG: 3857), and come from qimagery.information.qld.gov.au/. The project is in same CRS as rasters (better to project vectors, right?), and I'm using QGIS.

(Green = NP boundary, Blue = Creeks (Aus. Albers, GDA94 respectively)

I could use the GDAL Georeferencer, but this can subject the image to warping which I'd like to avoid. As far as I can tell they don't need any transformation, just a small shift on the lon/lat. It's probably frowned upon, but I just need to click and drag into position -- if my project/target CRS matches the raster I assume it's safe (especially on a small scale like this).

Could I edit the lon/lat in the world file/georef data?
(This ideal because I could perform in batch)

Or is it something simple like the ellipsoid?
(Currently set to WGS 84 to match raster)


texts:
I'm working with a point layer (many many features) that has some attributes, like the elevation of each point.

I created a polygon grid (5 km x 5 km each cell) with an unique ID and I made an intersection of the point with this polygon. This way each point that belong to a certain cell has now the unique ID.

Now, how can I display only the point with the maximum elevation for each different IDs? What is the expression I have to put in the field calculator of the rule renderer?

So in other word, how can I put into formula something like:

for each "ID" show max("ele_feet")


texts: Implementing if-then-else logic for number of points or area greater than x

texts:
I wish to alter this script (see below) that I use in model builder to choose which analysis to run next. I want to make "is state plane" process to run if "number of points (shapefile) is equal to or grater than 1". I also wish to have a second edition of the script that run the "is state plane" process if the area of a polygon (shapefile) has an area which is greater to or equal to 500m2 for example, how would I go about doing that?

  # Import modules
import arcpy
import sys
import traceback

# Set local variables
prj = ""
indata = "C:/ToolData/well.shp"
dsc = arcpy.Describe(indata)
sr = dsc.spatialReference
prj = sr.name.lower()

try:

   # check if indata is in StatePlane, has no PRJ, or one other than StatePlane
   if prj.find("_stateplane_") > -1:
       # Set the Is Unknown parameter to FALSE, and the Is StatePlane parameter to TRUE
       arcpy.SetParameterAsText(1,"false") #The first parameter refers to the "Is Unknown" variable
       arcpy.SetParameterAsText(2,"true") #The second parameter refers to the "Is StatePlane" variable
       arcpy.AddMessage("Coordinate system is StatePlane")

   elif prj == "unknown":
       # Set the Is Unknown parameter to TRUE, and the Is StatePlane parameter to FALSE
       arcpy.SetParameterAsText(1,"true")
       arcpy.SetParameterAsText(2,"false")
       arcpy.AddMessage("To continue, first define a coordinate system!")

   else:
       # Set the Is Unknown parameter to FALSE, and the Is StatePlane parameter to FALSE
       arcpy.SetParameterAsText(1,"false")
       arcpy.SetParameterAsText(2,"false")
       arcpy.AddMessage("Coordinate system is not StatePlane or Unknown")

except Exception as e:
    AddPrintMessage(e[0], 2)


texts: Installing ArcGIS Runtime SDK for .NET and what else to be able to display offline maps?

texts:
We have WPF application that is intended to run on Windows laptops carried in the field, that includes some mapping functionality. We need to be able to provide some mapping functionality, even when the laptops are offline.

With respect to mapping, our customers fall into three groups:

Customers who are heavy ESRI customers, who have defined their own maps using ESRI's software. For these, we expect we'll be able to display the maps in the application using ESRI's ArcGIS Runtime SDK for WPF.
Customers who aren't ESRI customers, but who have their own facility layers defined that they want to be able to load on the laptop and display in the application.
Customers who don't have their own facility layers defined, but need to be able to display streetmaps.

We know we can use ESRI's ArcGIS Runtime SDK for WPF for the first group of customers.

What we don't know is whether we can use the ArcGIS Runtime SDK for the second and third groups.

What we need for the second group is the ability for the customer to copy their map data onto their laptops, in some standard format, so that we can display their map data within the application.

What we need for the third group is to load a standard landbase onto their laptops so that we can display a street map within the application.

Looking at the ESRI documentation, it looks as if the Basic level would suffice for the second and third groups: [ArcGIS Runtime SDK for .NET - License your app][1].

But it's not clear to me what is involved in deploying an app.

With the ArcGIS Runtime SDK for .NET, what needs to be installed on the laptop, besides the application itself, to be able to display offline maps?  [1]: https://developers.arcgis.com/net/desktop/guide/license-your-app.htm


texts: Generate stream network for every pixel of flow direction as a polyline shapefile based on the raster cell value

texts: I have a flow direction raster (coded with ESRI flow directional encoding scheme)  and I want to generate stream network for every grid cell of flow direction raster as a polyline feature based on the pixel value in the flow direction without using the flow accumulation raster and threshold value for stream network generation. Is there any tools or package available to do that? I have tried using SAGA GIS for the purpose (Terrain analysis-channels-D8 flow analysis which is found only in earlier version of SAGA GIS) but it does not generate stream network for every grid cell. Sometimes, the stream network for multiple cells are combined to give only one polyline feature in SAGA.


texts: Does anyone know the equivalent of FloodRisk plugin for QGIS 2 in QGIS 3? Or alternatives to FloodRisk that work in QGIS 3?


texts:
BACKGROUND:

At 3.8.3 I have several line shapefiles, each with CRS EPSG:2927.  This is a CRS that I am familiar with, whose units are in feet.

When I use the field calculator with $length to create a new field containing each record's length the resulting values are in meters (or something close to that), when I was expecting feet.

This looks suspiciously like a CRS issue, but I confirmed the following:

All layers really are EPSG:2927
I restart QGIS and load just one of the layers.  Project Properties reports EPSG:2927 (plus EPSG:2927 is displayed in the lower-right corner).

Furthermore, if I run the Add geometry attributes tool, the results are correctly in feet.  When I use Arc Desktop's (10.7.1) field geometry calculator, the results are identical to the QGIS Add geometry attributes tool.

Finally, if I use the QGIS field calculator with $length on a meter-based CRS (EPSG:3338) layer, the results are close (within one percent), but not correct.  Again, the Add geometry attributes tool creates the meter correct values.

SUMMARY:

$length is not generating the correct line lengths.  It appears to ONLY generate output in meters, or something close to that.

In comparison, the QGIS Add geometry attributes tool and the Arc Desktop field geometry calculator generate identical correct results, confirmed with the QGIS measuring tool.

Either $length has a bug, or I'm using it incorrectly.  Thoughts?


texts: Fetching node(nd) reference data of way in Overpass query result

texts:
I wanted to check the feasibility of a doubt I am stuck with in Overpass query. I have created a simple overpass query to fetch data of one road,

[timeout:300];
(way[highway]({{bbox}}););
out geom;

and I got result as which gives 2 nodes with a ref id but those reference id are not present in details anywhere in the output.

<way id="256582511">
    <bounds minlat="53.6052931" minlon="10.2005595" maxlat="53.6053272" maxlon="10.2044249"/>
    <nd ref="331340954" lat="53.6052931" lon="10.2044249"/>
    <nd ref="331278759" lat="53.6053272" lon="10.2005595"/>
    <tag k="cycleway:left" v="track"/>
    <tag k="highway" v="tertiary"/>
    <tag k="maxspeed" v="70"/>
    <tag k="name" v="Hauptstraße"/>
    <tag k="ref" v="K 107"/>
    <tag k="surface" v="asphalt"/>
  </way>

Now for same bounding box-bbox, I exported data through osm.org I get data as below. This data contains details for the reference id, which is not the same in above query result.

<node version="2" lon="10.2040022" lat="53.6049916" uid="149203" user="uebi64" timestamp="2010-05-12T21:42:16Z" changeset="4681402" visible="true" id="731150522"/>
<node version="2" lon="10.2042253" lat="53.6049897" uid="149203" user="uebi64" timestamp="2010-05-12T21:42:16Z" changeset="4681402" visible="true" id="731150523"/>
<node version="2" lon="10.2040041" lat="53.6050694" uid="149203" user="uebi64" timestamp="2010-05-12T21:42:16Z" changeset="4681402" visible="true" id="731150532"/>
<node version="2" lon="10.2042272" lat="53.6050675" uid="149203" user="uebi64" timestamp="2010-05-12T21:42:16Z" changeset="4681402" visible="true" id="731150534"/>
<way version="2" uid="149203" user="uebi64" timestamp="2010-05-12T21:55:20Z" changeset="4681531" visible="true" id="58993281">
<nd ref="731150532"/>
<nd ref="731150522"/>
<nd ref="731150523"/>
<nd ref="731150534"/>
<nd ref="731150532"/>
<tag v="yes" k="building"/>
</way>

I am trying to use the data from Overpass query for further processing but because it has no data related to the reference id. Is it possible in the Overpass query to fetch the details of the node references?


texts: How can I convert polygon features stored in Shapefile, GeoJSON or GML (over WFS) into obj file

texts: I am currently working on an augmented reality application using ARKit. The required file formats are DAE or OBJ. My clients would like to display information stored within Shapefiles, GeoJSON or GML file formats. After much research I can't seem to find a simple solution on how to convert the polygons into the desired file formats. The best solution would be access to a free web-service that can complete the conversion for me. Mapbox provides some functionality with these files, however I am unaware if it can export to OBJ. Can anyone direct me to an existing solution (web-service) or a set of steps in order to complete the conversion?


texts: How to use an expression to count intersections in QGIS

texts:
I have a line layer and a point layer. The points intersect the lines. How do I write an expression for the Field Calculator in the line layer which will tell me how many points intersect each line?

I tried aggregate('point', 'count', intersects(geometry(@parent), $geometry)) but the result of this operation is '11' (total count of all points) for each line.


texts: Incorrect display of national characters from Spatialite in QGIS, Windows

texts:
I used this tutorial by underdark (thanks much :-) ) to convert my area of interest from OSM file (whole Czech Republic downloaded from Geofabrik.de server). Used OSGEO4W shell (saying "GDAL 1.11.0, released 2014/04/16") and following command:

C:\GIS\CZ_OSM_data_2014-07-22T20-21-02>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat 12.50004 50.46600 12.93001 50.15486 nejdek_OSM_2014_07.db czech-republic-l
atest.osm.pbf

I loaded the exported DB in QGIS and found out that the places labels are not displayed correctly in QGIS 2.4.0-Chugiak (GDAL/OGR 1.11.0, SpatiaLite 4.1.1) running on Windows 7 64bit while the same project is ok when opened in 2.4.0-Chugiak (GDAL/OGR 1.10.0, SpatiaLite 4.1.1) running on 32bit Kubuntu linux on my laptop.

Tried to find and encoding settings or something similar but no success...

Any advice or idea what am I doing wrong?


texts: How to merge/combine files in Mapshaper using combine-files command?

texts:
I'm trying to simplify several geojson and topojson files at once to output them into one file with a consistent level of detail. I'm using Mapshaper because I like it's map simplification features.

I see in the docs that it has a  combine-files feature:

  combine-files Import multiple files to separate layers with shared topology. Useful for generating a single TopoJSON file containing multiple geometry objects.

That's exactly what I want to do, but there's no example and everything I've tried doesn't work.

I've tried applying it in the input commands box before upload, at the same time as upload, after upload, with and without quotes, and from the console and output box after upload. They all just show the separate files as separate files.

How do I use combine-files to upload / input multiple files into one?

I noticed that the browser interface seems to automatically treat multiple files as layers, so I tried -merge-layers name="another layer" target="selected layer", or -merge-layers name="another layer", but it replaced the layer instead of merging.

Basically I just want to be able to see and simplify the contents of multiple files at once, to ensure consistency and no gaps between layers that border each other, then output them as (ideally) one file with one layer for each original file.


texts:
I downloaded a layer from https://www.arcgis.com/home/item.html?id=5064154431bf4df5aa444f3b863dfd89

I can add the layer to ArcGIS but I am not sure the file type.

How can I convert this to use a buffer tool on this (like 5-mile buffer from the wetlands)?


texts:
I would like to move a marker cursor at the same position and the same time as the default cursor over a Mapbox map when is moved active map cursor.
It's an application with two divs with a map in every div, but by default only appears one cursor over the active map, see JSFiddle: https://jsfiddle.net/6q2mdn3f/
I need to show two cursors, attached image with left map showing an arrow but, how can I do it?

  <title>2 Mapes - Mapbox</title>
  <meta charset="UTF-8">
  <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

  <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v1.13.0/mapbox-gl.css' rel='stylesheet' />


texts: Generate near tables from one layer to all in mxd

texts:
I am looking to create a script where user can choose a layer, then find all the layers within a set distance of it. I have set up a script to iterate through the mxd but I keep getting an error for the Generate Near Tables line saying my parameters are invalid and that my value is not a feature layer. I even added in a line to convert all to feature layers but this error won't go away. is it something wrong with my loop or input classes? I am new to coding so this is very overwhelming. This is the code I have been working with:

import arcpy
import sys

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, '')[0]
ctr1 = 1
ctr2 = 1

destPath = r"E:/File_Organization/Data/Working_Folder"
arcpy.env.overwriteOutput = True
for lyr in arcpy.mapping.ListLayers(mxd, '', df): # Loop through layers
# Output layer to shapefile
    outLyr = sys.argv[1]
    arcpy.MakeFeatureLayer_management(lyr, outLyr)
    for outLyr in arcpy.mapping.ListLayers(mxd, '', df):
        inFc = sys.argv[2]
        outT = sys.argv[3]
        distance = sys.argv[4]
        outL = sys.argv[5]
        outTable = outT + str(ctr1)
        FinalLyr = outL + str(ctr2)

        arcpy.GenerateNearTable_analysis(inFc, outLyr, outTable, distance, "LOCATION", "NO_ANGLE", "ALL", "", "GEODESIC")
        ctr1 = ctr1+1
        arcpy.MakeXYEventLayer_management(outTable, "NEAR_X", "NEAR_Y", FinalLyr)
        ctr2 = ctr2+1


texts: Set the feature class for a Python script to act upon in ArcMap

texts:
I want to use Python in ArcMap 10.5 Python window to give me the longest string in a field. I need to trim down some shapefile (they have to be shapefiles) file sizes. I want to reduce field lengths where possible. I am also going to run simplify tool. I am asking for Python help in this post, not ways to decrease shapefile size.

I don't understand how to specify which shapefile in my workspace for the script to run upon. Here is the script, most of it copy and pasted from another StackExchange post.

import arcpy

arcpy.env.workspace = "W:\WildCAD6\GIS"

#HOW DO I SET THE SHAPEFILE THAT I WILL NEXT CHOOSE FIELD FROM?

def longestWord(!FIELD!): #IS THIS WHERE I SET THE FIELD?
    count = 0    #You set the count to 0
    for i in !FIELD!: # Go through the whole list
        if len(i) > count: #Checking for the longest word(string)
            count = len(i)
            word = i
    return ("the longest string is " + word)


texts: How to use the OGR/GDAL configure options in QGIS Python environment

texts:
I'm looking for a solution how to use the OGR/GDAL configure options (--config) in the QGIS Python environment. http://trac.osgeo.org/gdal/wiki/ConfigOptions

Setting special configure options is necessary to load some data correctly (OGR Driver).

QgsVectorLayer(filename, title, "ogr")


texts: Saving to KML in QGIS producing bizarre results

texts:
I have a shapefile in WGS84  and I'm trying to save to a KML. Something is going wrong in this save and my KML file contains all of the polygons stacked in a sort of column shape.

Why is this happening and how can I solve it?


texts: ogr2ogr mssqlspatial string imports as a single special character

texts:
When importing a shapefile (in this case a public census tract shapefile) into SQL Server 2016 using ogr2ogr, the string fields are set to a single special character. The geography and integer fields import without issue.

The first import will be ʳ, the second will be ǿ, a new character each time I overwrite. The mssql field creates with the nvarchar(x) datatype and the appropriate length that matches the shape fields. I've also tried a .vrt file with specifics, but same thing happens.

ogr2ogr -overwrite -f "MSSQLSpatial" "MSSQL:Driver={SQL Server Native Client 11.0};server=(local);database=Spatial;trusted_connection=yes" tl_2010_06_zcta510.shp

Same unanswered issue here - Encoding issue with ogr2ogr and SQL Server

This happens on all shapefiles I try to import, no matter the source. String fields turn into single special characters when loaded. I am using the most recent GDAL 2.1.3, but have also tried the dev 2.2.0.

The data appears when querying "Results to Text", but has the special single char like in the snippet below.


texts: Formatting dynamic text element in Data Driven Pages?

texts:
Is there a way to format the Name element generated as dynamic text from Data Driven Pages?

This element serves as my title and it is too long and needs to be wrapped onto another line.


texts:
I wrote the following code to start using QGIS APIs. in a very short script seen below, the "Error: Algorithm not found" is thrown and "print processing.alglist()" is showing none.
Any idea why this is happening?

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
import nltk
import sys, os
import shapely
import math
from collections import defaultdict
from shapely.geometry import LineString
import json
import datetime
import time
import itertools
import glob
import sys

app = QApplication(sys.argv)
qgis_prefix = "C:/Program Files/QGIS Valmiera/apps/qgis"
QgsApplication.setPrefixPath(qgis_prefix, True)
QgsApplication.initQgis()
layer1 = QgsVectorLayer("D:/Ehsan/Test/InputData/SHP/OSM_roads_no_psv_latin1_id.shp",     "test" , "ogr")

if not layer1.isValid():
    print "Layer failed to load! :("
print layer1.isValid()

import processing
processing.alglist()
processing.runalg("qgis:clip","D:/Ehsan/Test/InputData/SHP/OSM_roads_no_psv_latin1_id.shp","D:/Ehsan/Test/InputData/temporary/tmpotile.shp","D:/Ehsan/Test/LM_Feat_test.shp")

Notice: the imports list is to show that I have almost everything imported and added to the path and pythonpath...

Please help me to figure out how to resolve it.
So far, the QGIS has seemed to be so buggy and problematic outside its python console.
@Lukas Graf


texts: How to correct uneven raster cell sizes?

texts:
I am conducting a spatial analysis in ArcGIS.   I have noticed my raster's are not 'aligning'.  I have data from several sources which I converted to raster, all the same projection (NAD 1983 UTM Z10).  My DEM is cell size is 24.21477324 m, and some other data sets are 24.215m.  I rounded but I do not think this 1000th of a difference should make a huge difference, however when I zoom in to view the cells, they are still off.  I am including pictures below.

I have tried resample tool, and making both 24.215 but this did not make any change.

How can I align these cells?? This is important for my work as I am calculating the cell values together and I am concerned that my output map has incorrect values.


texts: using "abs()" function in QGIS field calculator?

texts:
I'm pretty new to GIS in general and I'm trying to take an absolute value of the difference between two fields in field calculator when creating a new field, but using ABS () returns an invalid expression in Preview. ABS () is not listed on the supported expressions--is there a way around this?
Thanks!
John


texts: Is the ST_ prefix appropriate for functions not included in SQL/MM Part 3?

texts:
I was reading a thread on the geospatial extension for Presto in this Github issue, where a function, line_locate_point, was introduced. It was based off PostGIS's ST_LineLocatePoint function, which returns a float representing the fraction along a line of the closest point on that line to a given location.

The question was brought up why it was named line_locate_point and not ST_LineLocatePoint like the PostGIS version. The response was that this function does not exist in the SQL/MM Part 3 standard, and so it should not start with ST_.

Reading quickly through the standard, I don't see any comments on how to handle cases where you introduce a spatial function to your database that is not in the standard. Is the spirit of the ST_ prefix to differentiate spatial functions from non-spatial functions (as seems to be the case with PostGIS), or is it to indicate that the function complies with an equivalent function in SQL/MM Part 3?

Looking at the current state of Presto's API, I have to say that the latter approach looks less clean and introduces some confusion as to why the names aren't consistent, but perhaps this could be addressed by a simple note at the top.

My question, then, is whether there is some aspect of the standard that I'm overlooking that allows for extensions to it beyond the defined set of spatial objects, or alternatively, if this is expressly forbidden by some written or unwritten rule of following standards.


texts: Finding tables related to layer using ArcPy?

texts:
We have an ArcGIS map service that includes about 50 layers, and many of them have related tables. I need to build a schematic of this service that includes layer name, source feature class (or table), attributes that are included -- but I also need to find which table(s) are related to a given layer, which is where I am stuck.

I have access to the source MXD. So in theory, I can open each layer in ArcMap and look at it. But, with 50 layers, I'm looking for a Python method -- especially one that I can reuse if this MXD/service changes in upcoming months (layers or tables added/removed, e.g.) to re-create my list of referenced data.

I can use the arcpy.mapping module to list the layers/tables, but I do not know how to use it to ask about relationships to other layers/tables.


texts:
At QGIS startup I get the message: "QGIS couldn't find your GRASS installation.
Would you like to specify path (GISBASE) to your GRASS installation?" (Windows 7) No amount of rooting around the directory structure in the GRASS directories solves it.  Tried uninstalling/reinstalling, no joy, but I also noticed that my plugins were still there, so the uninstall does not appear to have been fully clean.

I have a perfectly working installation on another Windows 7 machine that I installed just last week, so no clue why there's a problem on this machine.

Ideas?


texts: Reprojecting Google Maps (latlong) to Irish Grid TM65 using spTransform in R?

texts:
I am trying to transform points taken from Google maps (as latitude and longitude) to Irish Grid TM65 using the spTransform in R, and I'm confused on something.

Everything I find on the net tells me google using a mercator projection as so:

"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

For the Irish Grid I am using a proj4 string extracted from Irish Census shapefiles

"+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +datum=ire65 +units=m +no_defs +ellps=mod_airy +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15"

So that is fine - I plugged all that into R, but was getting wacky answers. So then I fiddled about with different CRS in QGIS to let me see what was going on. Anyhow I found that instead of using the above projection for google, instead by using this projection for the CRS for the points taken from googlemaps then gave me correct projections on my Irish Grid CRS:
""+proj=longlat +datum=WGS84 +no_defs"

How come the projection the entire internet told me to use gave me wacky results and the one I picked almost randomly worked?

This worked in both QGIS and R by the way. I'm reasonably happy with the results - but clearly I'm missing something and I need to understand this properly for the work I am doing


texts: Running a function on subset of a raster image instead of whole image in R

texts:
I have a wide swath raster image. I want to detect cloudy pixels based on criteria such that a pixel in raster layer only consider its 25*25 neighborhood not the entire image.

My problem is here that I don't want to consider whole image when checking a pixel for cloudy.

I want a R code to do the following:

library(raster)
raster_layer = raster("path_to_raster_layer")
for (i in 1: nrow(raster_layer))
{
     for (j in 1:ncol(raster_layer))
     {
        if(raster_layer(i,j)  < ( ( Q1 - 1.5 * ( Q3- Q1 )) of a neighborhood
         of 25*25 )
         {  remove this pixel   }
        else { keep this pixel }
     }
}

Q1  = first quartile of neighborhood dataset
Q3 = third quartile of neighborhood dataset


texts: Efficient way to store multiple values of an x,y location in Python?

texts:
I am trying to read values from a stack of images and store them in a structure that will allow for easy readability. My current idea is to create a loop

for x in width:
    for y in length:
        for z in height:

and try appending the z value to a list at [x, y]. The issue with this technique is that I do not know how to prevent losing the x and y position through the iteration and it seems rather inefficient.

The input data will consist of z values extracted from tiff images. I am reading them by first creating a stack through a loop statement

for i, fname in enumerate(image_files):
img = cv2.imread(fname, cv2.CV_LOAD_IMAGE_UNCHANGED)
img_stack[:,:,i] = img

which I found here: https://stackoverflow.com/a/20056062/1342580 (CV2 is only returning nan for me but that's a different question)

After the stack is created I used to loop at the beginning of my question to go through all of the pixels (which at the moment is 3858 x 3742) and append the z value of an x,y location from each image to the x,y value (for example if I had 4 images in the stack the result at position 1, 1 would be [1, 1, Z1, Z2, Z3, Z4]. Then I append that to another list that will contain a list for every x,y position

Is there a better technique for doing this?


texts: keep getting error while trying to insert row using insert-cursor

texts:
I'm trying to add a text file of points to a feature class using insert-cursor and I can't seem to figure it out. I'm new to python and arcpy so this is a little confusing for me. the error message I get is below. if anything else is needed I will try and provide it.

the error message reads:

  Traceback (most recent call last):   File
  "L:\400\lab9data\InsertCursor_GPSPts.py", line 61, in
      theInsertCur.insertRow([newRow])   File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\arcobjects\arcobjects.py", line
  106, in insertRow
      return

convertArcObjectToPythonObject(self._arc_object.InsertRow(*gp_fixargs(args)))

RuntimeError: Error in InsertRow

#import arcpy
import arcpy

#set up an insert cursor to write to the feature class you created in part 1. Variable name should be theInsertCur.

#1 ----------Add one line of code below-----------------------------
theInsertCur = arcpy.InsertCursor("l://400//lab9data1.1//lab9data//cursor.gdb//GPS_Pts//myyGPSpoints")

#a counter used to skip first heading line
theCnt = 1

#open the file and read each line using a loop
#2 ---------------Change file path below------------------
for theLine in open("l://400//lab9data1.1//lab9data//GPS_collection.txt", 'r'):

    #We are skipping the heading line and then reading all other lines in the text file
    if theCnt > 1:

        #split the line into a list of three values
        threeValues = theLine.split(',')

        #print threeValues list to get an idea of what is happening
        print(threeValues)

        #get the x coordinate, y coordinate, and speed values from the threeValues list
        theX = float(threeValues[0])
        theY = float(threeValues[1])
        theSpeed = float(threeValues[2])

        #create a new point object with theX and theY as its coordinates. Variable name should be newPoint.
        #see Professional Library - Geoprocessing - The ArcPy site package - Classes - Point
        #3 ----------Add one line of code below-----------------------------
        newPoint = arcpy.Point(theX, theY)

        #add a row to InsertCursor. Variable name should be newRow.
        ##see Professional Library - Geoprocessing - The ArcPy site package - Functions - Cursors - InsertCursor
        #4 ----------Add one line of code below-----------------------------
        newRow = arcpy.InsertCursor("l://400//lab9data1.1//lab9data//cursor.gdb//GPS_Pts//myyGPSpoints")

        #set the shape field value to the newPoint. The shape field name is 'Shape'
        ##see Professional Library - Geoprocessing - The ArcPy site package - Classes - Row
        #5 ----------Add one line of code below-----------------------------
        newPoint = theInsertCur.newRow()

        #set the speed field value to theSpeed. The speed field name should be 'Speed'
        ##see Professional Library - Geoprocessing - The ArcPy site package - Classes - Row
        #6 ----------Add one line of code below-----------------------------
        theSpeed = theInsertCur.newRow()

        #insert newRow using theInsertCur
        ##see Professional Library - Geoprocessing - The ArcPy site package - Functions - Cursors - InsertCursor
        #7 ----------Add one line of code below-----------------------------
        theInsertCur.insertRow(newRow)

    #advance counter
    theCnt = theCnt + 1

#delete references to cursor and row to clean up
#uncomment the code below if you are confident you have set the variables above
#to match what is in the code below.
##see Professional Library - Geoprocessing - The ArcPy site package - Functions - Cursors - InsertCursor
#8 ----------Uncomment below-----------------------
del theInsertCur
del newRow


texts:
Has anyone determined process times in a question in GRASS?
I would like to know how much, on average, it could mean to cross a spatial join in GRASS, with two SQLite bases, of approximately 10 million records and another of 4 million.
I only cross three fields.
It is best to do it from GRASS, or perhaps from another method.
I get the impression that the time is surely excessive.


texts:
I'm trying to export positions of a set of points using pyqgis

I saw the numerous posts concerning exporting shp dataset to csv but the command used in all of them is not working.

How to write vector layer to csv file with geometry by PYTHON?
Create a csv with geometry as WKT in QGIS (and choosing the field delimiter)
Save as .csv with coordinates in both QGIS 2.2 and PyqGIS

My code line is:

QgsVectorFileWriter.writeAsVectorFormat(final_palmtrees_layer, result_path+'_trees_pos.csv', "utf-8", None, "CSV", layerOptions ='GEOMETRY=AS_XY')

There are two issues with this line

I wrote it in the console and i'm getting a "7" error code in the console, but I cannot find what it means.
I put it in a script and i'm getting the csv file but only the attributes, not the points position + excel tells that the file is corrupted


texts: Adding ArcGIS Online vector tile service in Leaflet map?

texts:
I am trying to add an arcgis online served vector tiles service to a Leaflet map. I am using ( https://github.com/Leaflet/Leaflet.VectorGrid ) to add the vector tiles to a Leaflet v1.3.1 map application.

I am able to add ESRI layers into the Leaflet map using the following simplified code:

var esriTilesUrl = "https://basemaps.arcgis.com/v1/arcgis/rest/services/World_Basemap/VectorTileServer/tile/{z}/{y}/{x}.pbf";

var esriTilesPbfLayer = L.vectorGrid.protobuf(esriTilesUrl);

esriTilesPbfLayer.addTo(map);

However, when I try to add my arcgis online vector tiles to the leaflet map using the same code

var esriTilesUrl =  "https://tiles.arcgis.com/tiles/l4TwMwwoiuEVRPw9/arcgis/rest/services/CaryHouseHolds/VectorTileServer/tile/{z}/{y}/{x}.pbf";

var esriTilesPbfLayer = L.vectorGrid.protobuf(esriTilesUrl);

esriTilesPbfLayer.addTo(map);

I get the following error:

GET https://tiles.arcgis.com/tiles/l4TwMwwoiuEVRPw9/arcgis/rest/services/CaryHouseHolds/VectorTileServer/tile/13/3222/2303.pbf 404 ()

What is the correct url to use to consume arcgis online vector tiles?


texts: GeoTools GML parser sorted features collection

texts:
I am using the following code to parse GML file in java:

File gmlFile = new File("path/to/gml/file");
InputStream gmlInStream = new FileInputStream(gmlFile);
GML gml = new GML(Version.WFS1_0);
SimpleFeatureCollection featureCollection = gml.decodeFeatureCollection(gmlInStream);
SimpleFeatureIterator iterator = featureCollection.features();
while (iterator.hasNext()) {
       SimpleFeature feature = iterator.next();
       //Do something with the feature
}
iterator.close();

I need to have the returned features collection be sorted by a specific feature's attribute. How to do that?


texts: ASCII file with latitude, longitude and data to GeoTiff using Python

texts: I have a .asc file with three columns of data: longitude, latitude and a data value.  Each latitude and longitude point represents the location of a centriod of a 1/8 degree unprojected grid.  I want to convert this file to a raster (.tif) using gdal in Python.  I then want to use projected datasets to subset this large grid.  Any thoughts on how to begin here?


texts: GeoServer SLD - size by scale and color by attribute

texts:
I would like to write an SLD that defines point size by zoom level using zoom denominators and color by a comparison operator. I am having trouble combining the two into rules because I cannot find evidence in the SLD cookbook that logical (And) operators are supported for anything other than combinations of spatial and comparison operators. Here is what I have for my first zoom level (denominator based off Bing maps scales, this part works fine), can anyone suggest a work-around or a missing step? [Note: the only thing that is changing in the rule filter is the color based off of whether the point belongs to Company A or Company B.] Thanks for any insight!

<NamedLayer>
<Name>default_point</Name>
<UserStyle>
  <Title>Default Point</Title>
  <Abstract>A sample style that draws a point</Abstract>
  <FeatureTypeStyle>
    <Rule>
      <Name>Local company A</Name>
      <Title>Local Company A</Title>
      <ogc:Filter>
        <ogc:And>
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>company</ogc:PropertyName>
            <ogc:Literal>A</ogc:Literal>
          </ogc:PropertyIsEqualTo>
          <MaxScaleDenominator>144447.93</MaxScaleDenominator>
        </ogc:And>
      </ogc:Filter>
      <PointSymbolizer>
        <Graphic>
          <Mark>
            <WellKnownName>circle</WellKnownName>
            <Fill>
              <CssParameter name="fill">#000080</CssParameter>
            </Fill>
          </Mark>
          <Size>8</Size>
        </Graphic>
      </PointSymbolizer>
    </Rule>
    <Rule>
      <Name>Local Company B</Name>
      <Title>Local Company B</Title>
      <ogc:Filter>
        <ogc:And>
          <MaxScaleDenominator>144447.93</MaxScaleDenominator>
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>company</ogc:PropertyName>
            <ogc:Literal>B</ogc:Literal>
          </ogc:PropertyIsEqualTo>
        </ogc:And>
      </ogc:Filter>
      <PointSymbolizer>
        <Graphic>
          <Mark>
            <WellKnownName>circle</WellKnownName>
            <Fill>
              <CssParameter name="fill">#696969</CssParameter>
            </Fill>
          </Mark>
          <Size>8</Size>
        </Graphic>
      </PointSymbolizer>
    </Rule>
  </FeatureTypeStyle>
</UserStyle>


texts: Merging two rasters into one image mosaic using ArcGIS Desktop?

texts:
I would like to merge two raster into one image. Both rasters have value range different from each other. On one hand, the values range from  (Image 1) 21.5 to 15.7. One the other hand are 18.15 to 13.3 (Image 2).

Image 1 is found outside of Image 2 and they do not have overlapping areas. But I want both images to be displayed as one with the same value range and color ramp.

I am using ArcGIS 10.


texts: GeoTools exception: failed to connect to EPSG database

texts:
I downloaded the GeoTools 11.2  binary and am attempting to get it to run in my Netbeans java application.
I tried to run this part of code:

...
SimpleFeatureSource source = dataStore.getFeatureSource("morteza:final");

    SimpleFeatureCollection fc = source.getFeatures();
    try {
        SimpleFeatureIterator iterator;
        iterator = fc.features();
        while (fc.features().hasNext()) {
            SimpleFeature sf = fc.features().next();
            Point geom = (Point) sf.getAttribute("the_geom");
            System.out.println(geom.getCoordinate().x);
            System.out.println(sf.getAttribute("the_geom"));
        }
    } catch (Exception ex) {
        System.out.println("hiiii");
        ex.printStackTrace();
    }
...

and I want to get features of wfs of geoserver.
at the line iterator = fc.features(); I faced this exception:

org.opengis.referencing.FactoryException: Failed to connect to the EPSG database.
at org.geotools.referencing.factory.epsg.ThreadedEpsgFactory.createBackingStore(ThreadedEpsgFactory.java:441)
at org.geotools.referencing.factory.DeferredAuthorityFactory.getBackingStore(DeferredAuthorityFactory.java:133)
at org.geotools.referencing.factory.BufferedAuthorityFactory.isAvailable(BufferedAuthorityFactory.java:235)
at org.geotools.referencing.factory.DeferredAuthorityFactory.isAvailable(DeferredAuthorityFactory.java:119)

......
...
Caused by: org.postgresql.util.PSQLException: Connection refused. Check that the hostname and port are correct and that the postmaster is accepting TCP/IP connections.

how can I fix this exception?
I don't want to fetch features from Postgresql but it tries to fetch from it.
please help!


texts: Apply "natural breaks" classification to a raster symbology in QGIS?

texts:
Does anyone know how can I apply "natural breaks" to a raster symbology in QGIS?

Because I only see the "Equal interval" option..

I'm using QGIS 1.9