texts:
I am building a model in ArcMap ModelBuilder that will take an input dataset folder, iterate through each shapefile in that folder and use that shapefile as an input for another tool. I then want to base the output file name of the new tool to be based on the shapefile name used in the input (from the iteration). One of the restrictions for the output filename is that it must be less than 8 characters.

Can I use "%FeatureClass%" or something similar when putting the filename path in the output filename textbox?


texts:
I'm very new to remote sensing in general, and aren't sure if I've made an error.

I've calculated EVI2 for a Landsat 8 and Sentinel 2 image. Both Top of Atmosphere. I've used Google Earth Engine to calculate it.
The range of EVI2 has come out as -0.72 to 2.06 for the Sentinel image, and -0.20-0.94 for the Landsat image.

Every paper I've seen which shows EVI2 maps has a scale of -1 to 1. I've tried using ArcMap's stretch tool, with minimum-maximum set to -1 and 1, but the resulting image is mostly a block of white because many of the values are bigger than 1.

I'm confused! Should I present my EVI2 as -0.72 to 2.06, which is a nice image visually?

I have a similar problem with NDVI, although all my NDVI values fall between -1 and 1, the maximums are beneath them.


texts: How do you determine when user is clicking the end point of an edit sketch using ArcObjects?

texts:
I am using vb.net working on ArcObjects. One requirement is to pop up a toolbox once the user click the end point (double click the point under editing status). So I want to know how can I know when user clicking the end point. Is it possible?

By the way, will this possible by using python?


texts: How to obtain the maximum two values of a field, without having to do it manually in the attribute table?

texts:
I have a layer of points corresponding to the vertices extracted from a polygon, and from these I need to obtain the two maximum values and save them in a new layer of points.

Because I have more than 150 layers of points I would like to know if it is possible to do it without having to open the table of attributes of each of these and select the values that interest me.

I tried using the field calculator using the expression "maximum(expression[,group_by][,filter])" but that would not have the second maximum.


texts:
I have 2 layers under EPSG:3857 and one map under EPSG:4326 to represent the world.
When I preview my aggregation of layers, here's what is shown :

In fact, this is what it should be displayed, on North America coast and Europe:

Check the scale which is insane, I had to zoom in a lot.

Why did my layers don't fit the world? My two layers are SQL View Layers and the config is the same :


texts: issues with projection Rgdal versus Manifold

texts:
I am pretty new to the topic and OK with R, but spatial data is still terra incognita

I got a Mapinfo file and i had to convert it to ESRI shp and then within r into spatial lines data frame and use with leaflet+OS maps epsg:3857

whatever i tried within R with spTransform (rgdal) it did not work - either was not visible or was visible with offset (approx 20 meters). Proj4 did not want to work with SLDF.
I took the file into Manifold, converted it to 3857, returned to R, transformed it to epsg:4269 and it works perfectly OK with leaflet + OS map 3857. But i have to automate this process to work within R, without any Manifold intervention.

I checked how the files differed:
i selected one item (the same of course) from the file and to make it easy to present created a centroid for it. Then I checked projection, xy, and latlong.

the item created from the file without using Manifold:
proj4string

"+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"

centroid:

528685 181836.2

centroid after spTransform to 4269

-0.1450149 51.52028

the item from the file pre-processed by Manifold:
proj4string

"+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"

centroid:

-16321.64 6713938

centroid after spTransform to 4269:

-0.1466198 51.52079

the latter is correct of course. Just to emphasize - it is the same file and the same item

Few questions:

any idea what happens in Manifold that fails in rgdal?
is there any mechanism in R that I can apply to fix the error (to match the 'with Manifold' output)? Maybe a kind of offset or work on prj files?

I tried to use the latter CRS for transformation

P5 <-spTransform(ptp.points1P, CRS( "+proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"))

and it returned still incorrect values
the xy is t least inverted but shifted

 -16142.98 6713847

the latlong is of course incorrect, but the same as without any attempt to transform

-0.1450149 51.52028

content of prj file not transformed - failing to correctly show

PROJCS["Transverse_Mercator",GEOGCS["GCS_Airy 1830",DATUM["D_unknown",SPHEROID["airy",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

content of prj manifold transformed file - good to go

PROJCS["Mercator_2SP",GEOGCS["GCS_unnamed ellipse",DATUM["D_unknown",SPHEROID["Unknown",6378137,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["standard_parallel_1",0],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]


texts: Removing tidal flat and coastline after applying Landsat Cloud Masking

texts:
I am a new Google Earth Engine user.
I have tried to remove clouds and cloud shadows using pixel_qa band.
However, sometimes the cloud masking makes tidal flat or coastline to remove.
I don't know why it happens.

How can I remove only clouds, not coastline and tidal flat?
Besides using pixel_qa band, is there another way to remove clouds?

Code:

function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2016-05-01', '2016-05-31')
                  .filterMetadata('WRS_PATH', 'equals', 116)
                  .filterMetadata('WRS_ROW', 'equals', 34)
                  .map(maskL8sr);

Map.addLayer(dataset)


texts: Capture GDAL projection and create matching output (OGR or GDAL) shapefiles?

texts:
Working in C++, using GDAL and OGR: Depending on what the user has supplied for image input, I need to output point (or other) shapefiles in the same projection as the input raster images.

The output shapefiles must be readable by ESRI as defined to be in the same projection as the input image.  We are working with many input images and output shapefiles, so the projection process must be automated.

How do I best transfer the input GDAL projection information to create a corresponding output shapefile of the same projection?

Do I have to go to OGR and use well-known text protocols to output the shapefile?

Can I instead use the GDAL projection information to write out the shapefile using GDAL somehow?


texts: ArcObjects IDataset.Copy causes program to crash

texts:
My implementation of IDataset.Copy continues to crash and I cannot figure out why. The program attempts to copy a shapefile from a temporary folder in one directory to a final folder in another directory. Here is the relevant code, where strOutVoronoiPath is a filepath containing the destination folder and g_temporary is a filepath containing the temporary folder:

        //create a shapefile workspace factory
        ESRI.ArcGIS.Geodatabase.IWorkspaceFactory pSWSF = new ShapefileWorkspaceFactoryClass();

        //create a temporary IWorkspace from Shapefile IWorkspaceFactory object
        ESRI.ArcGIS.Geodatabase.IWorkspace pTempWSP;
        pTempWSP = pSWSF.OpenFromFile(g_tempFolder, 0);

        //create a final IWorkspace from Shapefile IWorkspaceFactory object
        ESRI.ArcGIS.Geodatabase.IWorkspace pFinalWSP;
        pFinalWSP = pSWSF.OpenFromFile(strOutVoronoiPath, 0);

        //create a temporary IFeatureWorkspace from Shapefile IWorkspaceFactory object
        ESRI.ArcGIS.Geodatabase.IFeatureWorkspace pTempFWS;
        pTempFWS = (IFeatureWorkspace)pSWSF.OpenFromFile(g_tempFolder, 0);

        //create a final IFeatureWorkspace from Shapefile IWorkspaceFactory object
        ESRI.ArcGIS.Geodatabase.IFeatureWorkspace pFinalFWS;
        pFinalFWS = (IFeatureWorkspace)pSWSF.OpenFromFile(strOutVoronoiPath, 0);

        //Some code goes here where I do some casting and operations
        //to produce a shapefile inside the temporary folder. I now have
        //a shapefile I want in temporary folder (at g_tempfolder, which is a string
        //filepath) and I want to copy it over to a final folder at strOutVoronoiPath
        //I have taken this shapefile and stored the data from it inside an IFeatureClass.
        //called pFeatureClass. The following code attempts to take the feature class data and
        //cast it into IDataset, from which I call IDataset.Copy. The program crashes
        //on the final line, where I invoke the pDs.Copy() method.

         pFeatureClass.FeatureClass = pTempFWS.OpenFeatureClass(strOutVoronoiFile);
         pDs = (IDataset)pFeatureClass; //cast the feature class into a dataset
         pDs.Copy(strOutVoronoiFile, (IWorkspace)pFinalFWS);

I'm thinking I might have cast something wrong in the final three lines of code?


texts: The features don't appear at the first time in my exported PDF

texts:
my problem consists of exporting to PDF format. The first time I run my code (Text Editor of QGIS) the PDF only has the legend. But the second time, it has the features.
Changing the line of QgsMapLayerRegistry doesn't affect this issue. What do you think the problem is?

Here is my code:

from qgis.core import *
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from PyQt4 import Qt

municipios = QgsVectorLayer("ROUTE\\MUNICIPIOS.shp","Municipios","ogr")
QgsMapLayerRegistry.instance().addMapLayer(municipios)

el_renderer = iface.mapCanvas().mapRenderer()

#AJUSTE DE ZOOM EN UN CUADRO
rect = QgsRectangle(municipios.extent())
rect.scale(1.9)
#el_renderer = iface.mapCanvas().mapRenderer()
el_renderer.setExtent(rect)

composition = QgsComposition(el_renderer)
composition.setPlotStyle(QgsComposition.Print)
composition.PaperOrientation(QgsComposition.Portrait)

lst = [municipios.id()]  # add ID of every layer
el_renderer.setLayerSet(lst)
#el_renderer.setOutputSize(QSize(composition.paperWidth(), composition.paperHeight()), composition.printResolution())
el_renderer.setOutputSize(QSize(210,297), composition.printResolution())
x, y = 0,0
w, h = composition.paperWidth(), composition.paperHeight()
composer_map = QgsComposerMap(composition, x ,y, w, h)
composer_map.setOffset(-50, -60)
composer_map.setNewExtent(el_renderer.extent())
composition.addComposerMap(composer_map)
composition.addItem(composer_map)

print "\n QgsComposerMap funciona"

# LEGEND

legend = QgsComposerLegend(composition)
legend.model().setLayerSet(lst)
legend.setItemPosition(180, 230, True)
composition.addItem(legend)
print "\n leyenda actualizada"

#FINISH LEGEND

#la ultima parte

printer = QPrinter()
printer.setOutputFormat(QPrinter.PdfFormat)
printer.setPaperSize(QSizeF(210,297), QPrinter.Millimeter)
#printer.setOutputFileName("ROUTE\\my_pdf.pdf")
printer.setOutputFileName("ROUTE\\DATOS\\my_pdf.pdf")
#printer.setPaperSize(QSizeF(composition.paperWidth(), composition.paperHeight()), QPrinter.Millimeter)
printer.setFullPage(True)
printer.setColorMode(QPrinter.Color)
printer.setResolution(composition.printResolution())

pdfPainter = QPainter(printer)
paperRectMM = printer.pageRect(QPrinter.Millimeter)
paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
composition.render(pdfPainter, paperRectPixel, paperRectMM)
pdfPainter.end()

print "\n se ha exportado a PDF"


texts: Inputing postcode when geocoding with MMQGIS Web Service

texts:
I frequently use MMQGIS Plugin (Geocode >> Geocode CSV with Google / OpenStreetMap) in QGIS 2.18.10 (64bit) to geocode adresses from a CSV file. In particular for Germany, but I would imagine this issue cuts across country borders.

My CSV file contains a seperate column with postcode information, but there is no field to input it in the dialog window. I need to use the postcode information for geocoding. Should I concatenate it with another field like City (for example - 10117 Berlin)? Does anyone know where to input the postcode?


texts: How to style a vector layer in QGIS using hexadecimal color code stored in attribute?

texts:
I have a polygon shapefile storing hexadecimal color codes in a specific attribute.

  gid    |   label   |   color
---------+-----------+------------
1        | Crops     | #FF5F2C
2        | Forest    | #00FF00
3        | Pond      | #0000FF
4        | Forest    | #00FF00

I was looking for a way to say QGIS to read these values in order to style the polygons with them. I can't find anything about this.

I hoped it would be possible without deeping into the Python Console (which still frightens me...).


texts:
I am trying to write a plugin with a input control such as textbox on it. I looked around quite a bit for building  plugins. But I can only find examples of how to build a plugin with buttons/icons in its toolbar.

How can I add other widgets such as a textbox on my plugin's toolbar?

For the button, there isn't much to handle except for a callback when it's clicked. I imagine it would be quite different for a textbox as one needs to associate the input control to e.g. a global variable??

For now, I need interaction with the textbox (QLineEdit). I want to implement a plugin for going to the next feature where I can display the row number of the current record, and jump to any given row by entering its row number.

I am using QGIS2.99 and Qt5, the standard Plugin Builder and pb_tool to create and deploy the plugin template. I suppose I need to add the QLineEdit in initGui() of the main py file(as shown below). I don't know much about Qt5, and only figured I might need addUserInputWidget() and/or connect it with the add_action() function? But I don't know how to put the pieces together, and register the events etc.

Here is what I have so far:

def initGui(self):
    """Create the menu entries and toolbar icons inside the QGIS GUI."""

    # Create & add a textbox
    self.textbox = QLineEdit(self.iface.mainWindow())
    self.textbox.resize(80,20)

    self.txtAction = self.toolbar.addWidget(self.textbox)
    #should I use addUserInputWidget(self.textbox) #??

    self.txtAction.setToolTip(self.tr(u'Current Row Number'))
    #callback=self.runTextChange   how to register an event handler?

The UI shows but there is no interaction.

Can someone help explain how to link the QLineEdit to the toolbar so that I can get/change the text value?

Also, a minor issue I can't figure out is that the textbox takes up the full width of the toolbar no matter what I use in its .resize() method. How should I limit its width?


texts: How to subset temporal range of CRU TS 3.0 rasterstack in R?

texts:
I have downloaded open data of CRU TS 3.0 from the website https://data.ceda.ac.uk/badc/cru/data/cru_ts/cru_ts_3.00/data in which cru_ts_3_00.1901.2006.tmp.nc.gz was downloaded. I open this data in R through the following code
library(ncdf4)
getwd()
nc.tmp <- nc_open("cru_ts_3_00.1901.2006.tmp.nc")
print(nc.tmp)
tmp1 <- brick("cru_ts_3_00.1901.2006.tmp.nc", varname="tmp")
tmp1
plot(tmp1)

> tmp1
class      : RasterBrick
dimensions : 360, 720, 259200, 1272  (nrow, ncol, ncell, nlayers)
resolution : 0.5, 0.5  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : C:/Users/XX/cru_ts_3_00.1901.2006.tmp.nc
names      : X372, X373, X374, X375, X376, X377, X378, X379, X380, X381, X382, X383, X384, X385, X386, ...
z-value    : 372, 1643 (min, max)
varname    : tmp

Now, I would like to take a mean of rasters from year 1960 to 1990 (1960-1990). I can't seem to make sense out of the names in tmp1 as it starts with X372, X373 .... and ending at X1643. So I am confused that which year(containing months) is 1960,1961,1962 etc?


texts: Calculating the time taken to complete the overlay display in map?

texts: Currently, I am able to calculate the server response time by logging the difference of current date to initialize the map with the time when layers is added. But this gives me server response time. My POI is to calculate time taken when all WMS layers are loaded successfully(let's say render time) for certain bounding box.Is there any way I can find the  render time using javascript after the addlayer event is successfully executed for all layers?


texts: What's the syntax to skip optional saga:catchmentareaparallel parameters?

texts:
In PyQGis, I want to use saga:catchmentareaparallel algorithm, in the parameters there is an optional arguments, and I don't want to use them, what is the value should i put on the argument? i try to put null and "" but unfortunately it's not working.

processing.runalg('saga:catchmentareaparallel', elevation, sinkroute, weight, material, target, step, method, dolinear, linearthrs, linearthrs_grid, chdir_grid, convergence, carea, cheight, cslope, accu_tot, accu_left, accu_right, caspect, flwpath)


texts: Calculating geometry area in ArcGIS Desktop using ArcPy?

texts:
I have a polygon shapefile and I want to write a Python script to add a field named "area" and calculate its geometry areas.

The script is like this:

arcpy.AddField_management("wrfd01_layer_intersect","area","LONG","","","","","NON_NULLABLE","NON_REQUIRED")
arcpy.CalculateField_management("wrfd01_layer_intersect","area","!shape.area@squaremeters!")

When I run it the field adds, and you can see the process of calculating the geometry area.

However, in the attribute table, there are no values in the field "area".

How to solve this problem?


texts: Confusion matrices for two lists of raster objects

texts:
I have two lists of raster objects. One contains validation plots and the other contains results from a random forest binary classification of those plots.  I would like to loop through both lists with caret::confusionMatrix to assess accuracy of the classifier. However, some rasters only have one class present, and caret::confusionMatrix returns the

Error: there must be at least two levels in the data.


texts: Reclassifying string variable in ArcGIS Field Calculator with Python Parser?

texts:
I am trying to calculate a roadway pedestrian restriction based on its functional class, using if/else logic in Python. Every time I enter this, it only turns up 'No'. The categories are strings, and I am hoping to return strings. I have double-checked that the arguments for the values in FUNC_CLASS are accurate, so I am not sure why it will only turn up 'No's. I'm new to Python.

Codeblock:

def Reclass(function):
  if (function =="U / FWY OR EXP") or (function =="R / INTERSTATE") or (function =="U / INTERSTATE"):
    return 'Yes'
  else:
    return 'No'

Expression:

Reclass ( !FUNC_CLASS! )


texts: Merge Contour Boundary for Polygons QGIS

texts:
Hoping this is fairly simple but haven't been able to find much online. I am trying to create some elevation polygons from a DEM model.
Right now I am getting this after getting the contours, converting to multipart, then polygons. But you can see its messing up some of the unbounded contours

Is there any way I can merge a boundary box and then get the polygons to show up correctly? Based on elevation


texts: Convert GeoJSON with multiple Polygons and MultiPolygons into single MultiPolygon?

texts:
Is this even possible?

I have a GeoJSON file with multiple polygons and multipolygons, and I need to set styling on them (in Leaflet.js). I'd like to create a single, large MultiPolygon from the individual MultiPolygons and Polygons.

How can I do this (preferably in Node.js or at the Linux command-line)?

Things get really tricky when there's a MultiPolygon consisting of Polygons AND other Multipolygons.

Can this be "cleaned up" without any loss of information?


texts: ArcMAP path distance function handling of polygon and raster inputs

texts:
I am using a path distance function to calculate distances from a set of source cells on a land use raster. Unfortunately, my land use raster and my other inputs have different cell sizes. The land use raster has 300 m pixel sizes and my other inputs have ~90m pixels. I would prefer to avoid resampling. I would also like the entire area of the source pixel to be represented as 0 on the resulting distance raster. I will be pulling out individual pixels from the land use raster to use as source cells, and the easiest way to do this is to convert to points . However, given the size differences in the pixels, I'm not sure that it would be accurate to represent the source cells as points. I was thinking of using polygons instead (which can be expanded to the 300 m size of the pixel), but I'm unsure of how the tool will handle the polygons when calculating distance. If the tool will use the centroid of the polygon rather than the edges, then using a polygon is no improvement over using a point to represent the  area. I had also considered using a raster as an input, but, again, I'm not sure how the tool will handle the size differences.

Essentially: Does anyone know how exactly ESRI's distance tools handle different inputs during execution when the source inputs may not match the analysis cell size or may span across multiple pixels in the analysis cell size? How does the program handle resampling and/or source cell definitions in these cases?


texts: How to obtain id of the entity on the map?

texts:
Let's say I am looking the map where there is a lake or the forest. How to find out find an id of it? If the lake/forest is divided into smaller chunk, how I can select given chunk and get id of it?
Selecting chunks (or entire objects) is possible in OSM editor, I can see plenty of description there, but I don't see id of it (just for the record, I don't intend editing data, I used editor just to try to read the id).


texts: Extract population density data from a raster by shapefile

texts:
I'm trying to extract population density data from a raster by shapefile. I get NaN values, but I think my problem is further upstream at the raster level. I don't do much GIS work, so I'm likely making a rookie mistake here.

The population data source is interesting. It comes from an effort by Facebook to create population density maps from satellite imagery.

What I do below is grab the tif, use the raster package to load the file, and extract and average the density values by constituency (an administrative unit) with a shapefile I load.

The first problem is that the raster plot is blank (but the legend appears).

The second problem probably flows from the first: I am not able to extract any data.

library(raster)
library(velox)
library(sf)

# download the raster data
# https://data.humdata.org/dataset/highresolutionpopulationdensitymaps-ken
  download.file("https://data.humdata.org/dataset/2964b369-c10c-4b55-94a8-495de3fc9858/resource/1271b53b-30be-4c4b-aa27-1bd104506636/download/population_ken_2018-10-01.zip",
                destfile = 'popden.zip')

# unzip the file
  unzip(zipfile = "popden.zip",
        exdir = 'popden')

# load raster
  r <- raster("popden/population_ken_2018-10-01.tif")
  plot(r)

# download the constituency shapefile
# http://www.opendata.go.ke/datasets/constituencies-2009/data
  ken_con_sf <-
    read_sf("https://raw.githubusercontent.com/mikelmaron/kenya-election-data/master/data/constituencies.geojson") %>%
    st_transform(4326)
  plot(ken_con_sf)

# extract
  r_v <- velox(r)
  ken_con_den_mean <- r_v$extract(ken_con_sf, fun=mean)

@Spacedman correctly pointed out that my initial raster plot is not blank. Here's my new attempt to set breaks and visualize.

# get breaks as quantiles
# very slow so commenting out and giving values
  #rdf <- rasterToPoints(r) # very slow
  #rdf <- data.frame(rdf)
  #colnames(rdf) <- c("X","Y","den")
  #rb <- as.numeric(quantile(rdf$den, prob = seq(0, 1, length = 11)))

  rb <- c(0.002665759, 0.003351455, 0.005186487, 0.017159642, 0.048804275, 2.909804751, 3.859476157, 4.869144348, 5.943431054     7.815839071, 13091.488689846)

# plot
  plot(r, breaks = rb,  col = rainbow(100))

It's a little better, but still seems off to me.

To @Spacedman's second suggestion, I updated r_v$extract to run r_v$extract(ken_con_sf, fun=function(x){mean(x,na.rm=TRUE)}). This appears to have done the trick.


texts: Find points within a radius in ArcGis 10.2

texts:
I have a feature service which is created from a view in a geospatial geodatabase.
I want to retrieve all the point/geometries than fall within a specific distance.
I am using arcgis 10.2
I know that i could do it using the JS api like in this  Example.
However I would prefer to have geo-processing task that would retrieve the a point and calculate which point fall within a buffer of a specific radius around this point.
The problem is that I cannot find a tool to give me the points within the buffer.
Is there any simple way to do that ?

In other words my inputs are

the radius
A feature class which in my case in a view in a geodatabase
a point around which I want to find the nearby points.

And my output should be all the nearby points around the input point


texts:
This is a really specific question, but anyone with experience trying to pull in outside information into a map (live or static) can probably add something to the question.

A long winded post is here  http://www.thadwester.com/1/post/2012/04/mapping-attainable-bitrates-for-rural-telcommunication.html

Does anyone have experience mapping attainable bitrates?

OR

how were you successful in pulling data from a non GIS source and integrating it into your map.


texts: Merging two isolines into one in PostGIS

texts:
How can I sum two isolines A and B (closed linestrings with values) into one isoline C with summed values on common parts?
For each line string there is a value and I would like to sum on common part.

If it is easier I can try to sum isopolygons made from isolines.


texts: Where is the log for advanced Python calculator errors?

texts: There is a specific log file for Advanced Python calculator in QGIS 2.4.0? I got the tab message below point to a log file there I wasn't able to find by myself.


texts: Creating point shapefile from GeoTIFF using GDAL/OGR?

texts:
I have a 0.05 deg global maps (~8000 files) in GeoTiff format. It's a kind of a mask containig "1" and "0". I need to create a point shapefile for each raster with points located over "1" values of raster.

I tried to do this in QGIS via this tutorial http://www.qgistutorials.com/en/docs/sampling_raster_data.html, but creating a shapefile with ~26M of points and sampling it over each raster is quite slow.

Is there a better way to do this, preferably with gdal/ogr?


texts: How can i add a north arrow in QgsComposition?

texts:
I'm trying to add an svg image of north arrow in QgsComposition using this code:

mapRenderer = iface.mapCanvas().mapRenderer()
theComposition = QgsComposition(mapRenderer)
theComposition.setPlotStyle(QgsComposition.Print)
northarrowIcon = QgsComposerPicture(theComposition)
northarrowIcon.setPos(QPointF(50,110))
northarrowIcon.setPictureFile("./svg/north_arrows/north-arrow_3_simple_symmetric_triangular.svg")
theComposition.addItem(northarrowIcon)

i don't have an execution problem but the arrow is not shown.

Can someone check and give me an idea what can be the cause of that.


texts: Pattern recognition - machine recognition of plausible shape of a polygon

texts:
I have a schoolproject where the task is to locate all possible shapes (of a certain complex polygonic structure) in a layer of several heterogenous polygons.

And now the task is to find all the ones that assimilate (more or less)

Starting from above the exact tools, I think we are dealing with a machine-learning concept here. I've looked at Grass-GIS GeoMorphon and it seems that it has the potential to start using as a tool. However I have never tried such an advanced task before. I have also considered OpenCV to maybe provide me with a solution, however the libraries for pattern recognition is based on raster photos and perhaps not suited for a vector-scaled dataset.

So far I've been filtering out all polygons <= 6 vertices after having done a (dangerous) simplification of the polygons. That reduced the number of polygons left to about 60 %. But there is still a long way to go and I have no idea of how to progress.


texts: QGIS 3.4 Processing Python Scripts folder default configuration is lost

texts:
I have lost the default configuration of the 'Processing Python Scripts folder' on QGIS 3.4 and in consecuence some Python processing scripts and plugins cannot work properly because they cannot find the path to their resources.

Does any one knows wich are the default configuration paths of the processing scripts folder?

This is the error msg:

2018-12-19T10:51:17     WARNING    warning:C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\PyQt5\uic\objcreator.py:152: DeprecationWarning: 'U' mode is deprecated
              plugin = open(filename, 'rU')

         traceback: File "", line 1, in
          File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\utils.py", line 335, in startPlugin
          plugins[packageName] = package.classFactory(iface)
          File "D:/Users/gasato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\DigitizingTools\__init__.py", line 28, in classFactory
          from .digitizingtools import DigitizingTools
          File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\utils.py", line 672, in _import
          mod = _builtin_import(name, globals, locals, fromlist, level)
          File "", line 983, in _find_and_load
          File "", line 967, in _find_and_load_unlocked
          File "", line 677, in _load_unlocked
          File "", line 728, in exec_module
          File "", line 219, in _call_with_frames_removed
          File "D:/Users/gasato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\DigitizingTools\digitizingtools.py", line 33, in
          from .dtDialog import DigitizingToolsAbout
          File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\utils.py", line 672, in _import
          mod = _builtin_import(name, globals, locals, fromlist, level)
          File "", line 983, in _find_and_load
          File "", line 967, in _find_and_load_unlocked
          File "", line 677, in _load_unlocked
          File "", line 728, in exec_module
          File "", line 219, in _call_with_frames_removed
          File "D:/Users/gasato/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\DigitizingTools\dtDialog.py", line 27, in
          os.path.dirname(__file__), 'ui_about.ui'))
          File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\PyQt5\uic\__init__.py", line 198, in loadUiType
          winfo = compiler.UICompiler().compileUi(uifile, code_string, from_imports, resource_suffix, import_from)
          File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\PyQt5\uic\Compiler\compiler.py", line 53, in __init__
          qtproxies.QtWidgets, CompilerCreatorPolicy())
          File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\PyQt5\uic\uiparser.py", line 177, in __init__
          self.factory = QObjectCreator(creatorPolicy)
          File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\PyQt5\uic\objcreator.py", line 85, in __init__
          if self.load_plugin(filename, plugin_globals, plugin_locals):
          File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\PyQt5\uic\objcreator.py", line 152, in load_plugin
          plugin = open(filename, 'rU')


texts: Adding a layer on top of eSearch widget using web appbuilder

texts:
I added one layer on top of eSearch widget and its working in fine when we launch in WAB. I saved and  downloaded the code of updated eSearch widget. I found the same code as before adding layer to it and also I have seen in server/apps/2/widgets/eSearch/config.json, but 'config.json' file showing the same code with some other services layers.

Where I can find the updated code after adding layers ?


texts:
What is the best way, if any, to downgrade LAS files to an older version?

LAS1.4 works well with ArcGIS 10.3.1, however our clients that are using ArcGIS 10.2.2 or earlier are unable to load LAS1.4 files into an LAS dataset; the file format is undetectable.

Prefer a method in ArcGIS or a free open-source software.

Also, if downgrading the LAS1.4 to LAS1.3 or earlier is even possible, what are the implications of such a process?
Will there be an expected degradation of the data itself or certain statistics that are lost?


texts: Why is the first element of my WFS layer None instead of a feature?

texts:
I have made a WFS connection and got a layer from it. I want to use the first feature from that layer.

driver = ogr.GetDriverByName("WFS")
wfs_url = "WFS:http://example.com/wfs"
wfs = driver.Open(wfs_url)
layer = wfs.GetLayerByName("mylayer")

Now I want to use the first feature from it.

>>> print(sessionlayer[0])
None
>>> type(sessionlayer[0])
<class 'NoneType'>

Why does layer[0] not give me the first feature?
What is that layer[0] element?

If I loop over the features of the layer, I do not get that weird None element:

for i, feature in enumerate(sessionlayer):
  print(i, feature)

gives me

0 <osgeo.ogr.Feature; proxy of <Swig Object of type 'OGRFeatureShadow *' at 0x7f33cca12690> >
1 <osgeo.ogr.Feature; proxy of <Swig Object of type 'OGRFeatureShadow *' at 0x7f33cca126c0> >
...
...

Shouldn't the first printed element be that None again?


texts: Line appears broken into two geometries in QGIS

texts:
I have a shapefile with line features. One of my lines was causing problems in a python script and in the QGIS interface it appeared as two segments (one long one with a brief gap between two nodes).

I used the calculator with num_geometries( $geometries) and it showed this entry has two geometries (all the others had one). I imagine this is the problem. Is there any way to fix this and merge them back to a single geometry? Deleting and re-digitizing would be a lot of effort. I didn't see anything in the Validity checker and it didn't complain, but that might be because I messed up a parameter.

The merge tools in QGIS digitizing bar is not relevant. It is for merging multiple features, not multiple gometries in one feature -- it complains I have only one feature selected.


texts: I want to convert openlayers interface code to dojox code for pie chart use.So can we use same controls of openlayers map in dojox.If yes, how?

texts: Is their any way to add same openlayers map control to dojox map like navigation,panel,panzoom


texts: WFS GetFeature request blank XML for some FeatureTypes

texts:
I have following URL request to GeoServer. The first request results blank page but the second page results in number of feature count. Why it there blank XML for some FeatureTypes while other gives number of features.

http://192.168.70.65:7070/geoserver/wfs?request=GetFeature&typename=GisData:sites&version=1.1.0&resultType=results&service=WFS&

http://192.168.70.65:7070/geoserver/wfs?request=GetFeature&typename=NCellGisData:sites&version=1.1.0&resultType=hits&service=WFS&


texts: Geopandas ImportError: DLL load failed (in pyviz tutorial geographic data examples) win7

texts:
I'm trying to run the geographic data examples from pyviz tutorial. I'm getting the error list in modules below. My Geopandas module work fine (and the others dependent modules) with 3.5 python version in another project made with pycharm interpreter. I think it could be a python version problem...not sure.

ImportError                               Traceback (most recent call last)
 in
    ----> 1 import geopandas as gpd
    2 world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
    3 world.head()

C:\ProgramData\Miniconda3\lib\site-packages\geopandas__init__.py in
    3 from geopandas.geodataframe import points_from_xy
    4
    ----> 5 from geopandas.io.file import read_file
    6 from geopandas.io.sql import read_postgis
    7 from geopandas.tools import sjoin

C:\ProgramData\Miniconda3\lib\site-packages\geopandas\io\file.py in
    2 from distutils.version import LooseVersion
    3
    ----> 4 import fiona
    5 import numpy as np
    6

C:\ProgramData\Miniconda3\lib\site-packages\fiona__init__.py in
   67 from six import string_types
   68
    ---> 69 from fiona.collection import Collection, BytesCollection, vsi_path
    70 from fiona._drivers import driver_count, GDALEnv
    71 from fiona.drvsupport import supported_drivers

C:\ProgramData\Miniconda3\lib\site-packages\fiona\collection.py in
    7
    8 from fiona import compat
    ----> 9 from fiona.ogrext import Iterator, ItemsIterator, KeysIterator
   10 from fiona.ogrext import Session, WritingSession
   11 from fiona.ogrext import (

ImportError: DLL load failed: Impossível localizar o módulo especifica


texts: Efficient way to renumber source and target using pgRouting

texts:
pgRouting is a powerful tool for users to handle routing tasks, and I am now use pgRouting 2.1 to route a river (line shapefile).

Before starting to route, the first thing to do is always to create a routable table, and I use pgr_createTopology to create source and target of any edge.

Here is the problem, if I want every direction of 'from-source-to-target' in all edges to be the same, I need to check and renumber them manually.

Is there an efficient way to achieve this?

P.S. I have seen the following posts, and still have no idea to solve this problem.
PGR_CREATETOPOLOGY

How does assign_vertex_id() in pgRouting work?

pgRouting assign number of source and target


texts: ArcSDE Feature View in a GeoDatabase other than the Feature's Database

texts:
here is the problem I have been struggling with:

We have our spatial features in an SDE enabled DB (call it DB1) and have a client DB which is SDE enabled as well (meaning geo-DB call it DB2) and have some non spatial tables too.

Without copying features & using spatial features from DB1, we would like to generate Spatial Feature Views in DB2.

Is that possible? What options do we have to come over this?

Thanks a lot in advance


texts:
I am quite new to R and am having some trouble setting up my data for some point pattern analysis.
What I want to do: conduct a point pattern analysis on NYC arrest data and see if there exists a spatial dependence between arrests and Covid-19 cases.
What I've done so far: downloaded data in the form of shapefiles
https://data.cityofnewyork.us/City-Government/Borough-Boundaries/tqmj-j8zm (the ZIP code boundaries)
https://www1.nyc.gov/site/nypd/stats/crime-statistics/citywide-crime-stats.page (year to date data for arrests in NYC by zip code)
Code:
library(readxl)
library(rgdal) #Brings Spatial Data in R
library(spatstat) # Spatial Statistics
library(lattice) #Graphing
library(maptools)
library(raster)
library(ggplot2)
library(RColorBrewer)
library(broom)

# Load nyc zip code boundary polygon shapefile
s <- readOGR("/Users/my_name/Documents/project/zip","zip")
nyc <- as(s,"owin")

### OGR data source with driver: ESRI Shapefile
Source: "/Users/my_name/Documents/project/zip", layer: "zip"
with 263 features

# Load nyc arrests point feature shapefile
> s <- readOGR("/Users/my_name/Documents/project/nycarrests/","geo1")

### OGR data source with driver: ESRI Shapefile
Source: "/Users/my_name/Documents/project/nycarrests", layer: "geo1"
with 103376 features
It has 19 fields

#Converting the dataset into a point pattern
arrests <- as(s,"ppp”)

### Error in as.ppp.SpatialPointsDataFrame(from) :
  Only projected coordinates may be converted to spatstat class objects

This gave me the error above.
I know the error has to do with the coordinates not being in the Cartesian coordinates.
How can I convert my sp object to have (projected) Cartesian coordinates in order to convert it to a point pattern (Poisson point process) ?


texts: How to Project or Reproject Map Using csv file data in R?

texts:
I have a csv file loaded in R as the follow example below. I trying to assign projection. Then, plot map and draw pie for species variable for each grid cell.

Month    long      lat  variable   value
3       48.0000 -12.8800 SKJTotal  10.990
3       47.8800 -13.0800 SKJTotal   5.127
3       48.0000 -13.0000 YFTTotal   4.576
3       48.0000 -13.0200 SKJTotal   4.240
5       42.6800 -12.5700 SKJTotal  10.838
5       41.8700 -11.4800 SKJTotal  13.615
5       45.3200 -14.6300 BETTotal  13.615
5       45.6800 -14.1300 SKJTotal  21.675
5       42.5500 -11.7200 SKJTotal   6.130
5       45.5000 -13.7000 BETTotal  43.350
5       41.7800 -11.6300 YFTTotal  74.878

I used the codes below to assin the projection.

coordinates(mpma) = ~long+lat # trying to assign the projection
proj4string(mpma) <- CRS("+proj=longlat +datum=WGS84")
is.projected(mpma)
[1] FALSE

The summary from overall csv file this one below.

summary(mpma)
Object of class SpatialPointsDataFrame
Coordinates:
        min   max
long  34.00 49.35
lat  -27.33 -9.93
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0]
Number of points: 27054
Data attributes:
     variable        value
 SKJTotal:9018   Min.   :  0.0000
 YFTTotal:9018   1st Qu.:  0.0000
 BETTotal:9018   Median :  0.4575
                 Mean   :  6.1566
                 3rd Qu.:  6.0000
                 Max.   :412.4290

Then, I try to plot the map using normal plot as follow:

plot(mpma$long,mpma$lat)
map(add = T,fill=T,col="grey")

The returning map is not projected.

How do I fix the projection?

I have been tried many codes without success.  I also tried to use the world basemap setting the long=c(30,50) and lat=c(-30,-9). But here the output map as very small. And I can't find the code to zoom or increase the resolutions. I intend to draw grids and draw pies for each grid cells.


texts: GRASS GIS, r.li modules: How to fix "ERROR: Illegal areas disposition"?

texts:
While trying to run commands in examples of GRASS GIS r.li module, I get the error.
This is the output I get from commands in the examples
grass_cmd_output(Language:Japanese)
(Fri Sep  4 13:12:18 2020)
g.gui.rlisetup
(Fri Sep  4 13:12:21 2020) コマンド終了 (2 秒)
(Fri Sep  4 13:12:24 2020)
g.region raster=geology_30m -p
projection: 99 (Lambert Conformal Conic)
zone:       0
datum:      nad83
ellipsoid:  a=6378137 es=0.006694380022900787
north:      228500
south:      215000
west:       630000
east:       645000
nsres:      30
ewres:      30
rows:       450
cols:       500
cells:      225000
(Fri Sep  4 13:12:24 2020) コマンド終了 (0 秒)
(Fri Sep  4 13:12:27 2020)
r.li.patchdensity input=geology_30m conf=my_conf out=patchdens
ERROR: Illegal areas disposition
(Fri Sep  4 13:12:28 2020) コマンド終了 (0 秒)

Following this article, there was the similar(=same?) problem only on Windows OS. I'm using Windows 10 64bit and GRASS 7.8.3 (Language:Japanese).
Could you suggest any solutions to fix this error?


texts: Using carto mobile iOS sdk download offline map having satellite map style instead of normal map style

texts: Is there any way to use satellite map style for maps downloaded for offline usage using Carto Mobile SDK?


texts: Browser display nothing when debugging leaflet api file to connect to geoserver

texts:
actually here is url of geoserver on which my map is published:
http://localhost:8081/geoserver/test/wms?service=WMS&version=1.1.0&request=GetMap&layers=test:city_limit&styles=&bbox=-114.315796136219,50.842822033593,-113.859905274451,51.2124253170889&width=768&height=622&srs=EPSG:4326&format=application/openlayers

And i have written a html code using leaflet api to connect to geoserver:

<!DOCTYPE html>

    <title>WMS example - Leaflet</title>

    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">

    <link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />

    <link rel="stylesheet" href="https://npmcdn.com/leaflet@1.0.0-rc.3/dist/leaflet.css" />


texts: Downloading GitHub QML style file for use in QGIS

texts:
I am trying to download a QML-style from GitHub and add it to QGIS styles for use, this is how it looks like:

How can I do this? I'd really like to use this style.


texts: Where does "Update measurements" put its updated measurements? (MeemsTools)

texts:
I have several .shp files loaded into MapWindows GIS.  I go to plug-ins > Meems Tools and click on "UpdateMeasurement Plugin", and it loads.  I click on the Update measurements, and in the the Update measurements window that appears, I click the "Update measurements" button, without changing any of the options in the window.  The "select attribute" section grays out, and the bottom of the window lists "Finished updating measurement values", next to the "MapWindow Wiki help" link.

  Where does it put these updated measurement values?

possibly relevant information:

I'm using Windows 7 version 6.1 .

The files are loaded as multiple layers in the same window of MapWindow GIS.

The files are simplifications using MapShaper.org of the "district" files from "Boundary" on http://hawaii.wr.usgs.gov/maui/data.html .


texts:
Using geoserver I have some layers in CRS EPSG: 20935 (Arc 1950 UTM 35S). In order for them to work with base map layers (in leaflet) I set the native projection to 20395, the declared projection to 4326 (WGS 84) and ask geoserver to reproject to declared. This works and I can see my layers on leaflet map. The base maps are in web mercator but as I understand it 4326 will work.

However, using ArcGIS doing the same reprojection there are ~9 options as to which transformation to use. Given the location of my data number 7 would be most appropriate. Testing ArcGIS's "layer to kml" tool I discovered that ArcGIS defaults to reproject using the first transformation (ARc1950 to WGS84(1)) which leads to an East West shift of 30 - 70 metres compared to transformation(7). So in Arc it is necessary to reproject to WGS84 first (selecting the desired transformation). How can I find out what transformation geoserver is using?


texts: Code to get featureclass name from existing topology

texts:
I am creating and validating topology through vb.net and arcobject code.Things are running fine as I am able to add rule and validate it.

Problem is when once the program is run and topology is created. Second time when I am running my code, I have written that if topology exist then delete it and create new topology but I observed that deleting topology every time will corrupt memory.

So Second approach I used is to first delete rule in the existing topology and then featureclass from it and not the topology.

Then new rule and featureclasses are added to same existing topology.
But the problem is I am not able to get the featureclass name in the existing topology. How can I get the name of feature class participated in topology?

The code is as follows-:

If pdataset.type=esridatasettype.featureclass then
    Dim fcnme as string
    Dim fcarray as Arraylist=New Arraylist
    fcnme=pdataset.name
    fcarray.add(fcnme)
    for i as integer=0 to fcarray.count-1
        pname=fcarray.item(i).tostring
        Dim gp as Geoprocessor=New geoprocessor
        Dim removefc as DataManagementTool.RemoveFeatureclassFromtopology=New DataManagementTool.RemoveFeatureclassFromtopology
        removefc.in_topology=topologyname
        removefc.in_featureclassnme=pname
        Dim sv as object=nothing
        try
           GP.execute(removefc,nothing)
        Catch ex as exception
        end try


texts: How to update only a single vector layer in QGIS map canvas?

texts: I want to update a single vector layer as I am working on 8 vector layers and they are getting fast data from serial port so for updating I use mapCanvas().refresh() which refresh all the layer and I have almost 20 static S57 layer so the refresh takes a long time.


texts: Grouping and sorting a list of equipment using an aggregate function in the print composer

texts:
I'm using QGIS as a spatial database for a construction site and we are now checking the installed equipment. The inspection requires us to have a list of the equipment ID to inspect for each type of line (for example: all valves and pipes for process water in a certain building).

In order to generate the inspection sheets I use the atlas function of QGIS (which then allows me to generate reports as well). The atlas uses an inspection layer but the ID list is generated from an equipment layer which is why I use an aggregate function to generate the list and not a concatenate function.

I managed to generate the correct lists using the following code:

aggregate('Equipment nomenclature final','concatenate',expression:="COMPLETE TAG",filter:= ("Zone" is  attribute(  @atlas_feature , 'Zone' ) and  "Fluid  code" is attribute(  @atlas_feature , 'Fluid code' ) and  "PID" is attribute(  @atlas_feature , 'PID' )) ,concatenator:='; ')

And the result looks like this:

However the IDs aren't in a convenient order for readability on site and the aggregate function doesn't accept group_by parameters. I would like to have them grouped by 'Element Code', the first 2 or 3 letters of the ID (ex: EID, AV, AS, etc.) and in incremental numbers for the 3rd part of the ID. Ideally, the list would even return to next line when the Element code changes:

23 AS 011, 23 AS 012
23 AV 010, 23 AV 011
etc.


texts: This example http://openlayers.org/en/v3.5.0/examples/export-map.html on OpenLayers 3 does not work if you add more than two layers. I'm reading my layers from MapServer and the error is Uncaught DOMException: Failed to execute 'toDataURL' on 'HTMLCanvasElement': Tainted canvases may not be exported. Can you assist


texts:
I am currently writing a app in R / Shiny where I use Leaflet.  Background layers are WMS and it works fine for most services without specifying any projection anywhere.  When using fiddler I see that the bounding box is specified in meters which is logical since Leaflet uses ESPG 3857.

Then to the problem: When I try to access NASA's WMS services with addWMSTiles I dont get it to work since NASA's wmses require a lat/lon bounding box.

Does not work:
https://neo.sci.gsfc.nasa.gov/wms/wms?service=WMS&request=GetMap&layers=MOD_LSTD_CLIM_M&styles=&format=image%2Fpng&transparent=true&version=1.3.0&CRS=CRS%3A84&width=256&height=256&bbox=-20037508.342789244,-20037508.34278071,20037508.342789244,20037508.34278071

Works:
https://neo.sci.gsfc.nasa.gov/wms/wms?service=WMS&request=GetMap&layers=MOD_LSTD_CLIM_M&styles=&format=image%2Fpng&transparent=true&version=1.3.0&CRS=CRS%3A84&width=256&height=256&bbox=-180,-90,180,90

I have tried to define the projection and force Leaflet to use this projection without success.

Defining WGS84

epsg4326 <- leafletCRS(
 crsClass = "L.Proj.CRS",
 code = "EPSG:4326",
 proj4def = "+proj=longlat +datum=WGS84 +no_defs ",
 resolutions = c(8192, 4096, 2048, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2, 1, 0.5 ),
 origin =c(0, 0))

Then initiating Leaflet with:

leaflet(options=leafletOptions(crs=epsg4326, worldCopyJump = FALSE)) %>%

So the question is, how do I get Leaflet to send a bounding box in lat/lon??


texts: Creating polygons for a choropleth map using PostGIS

texts:
I have a grid of points, each with a value to it as shown in the image below. The data was retrieved from PostGIS and rendered in Quantum GIS. The final result should be a polygons of a choropleth map that will be used to render on the browser.

Question: How can PostGIS use the grid of weighted points to create the choropleth polygons?

Initial

Desired result

Polygons to be created in PostGIS, rendered in Leaflet. Classes group according to some rule like:

d > 1000 ? '#800026' :
d > 500  ? '#BD0026' :
d > 200  ? '#E31A1C' :
d > 100  ? '#FC4E2A' :
d > 50   ? '#FD8D3C' :
d > 20   ? '#FEB24C' :
d > 10   ? '#FED976' :


texts: What kind of QT widget is the QGIS drawing area?

texts:
I mean when open QGIS , and add a vector layer or other type of data. The content will show on the "Drawing Area", i want to ask, the drwaing area is some type of widget or just in the QMainWindow?

I download the qgis source code(1.8 or 1.9). and searche the src/ui/ folder, do not find any named similar to  "drawing "   *.ui file. And i use qt-creator and open the "qgisapp.ui" and judge by the size of the window, the drawing area is in this window. But i am not sure, because i am not familiar with qgis code.

What i want to do is want to mime QGIS to code little GIS drawing toy using Qt. So anyone could answer my question?


texts: How to calculate distance points sequentially on the same shapefile?

texts:
I have a points layer that came from a flight path and I want to calculate the distance between each point sequentially. Ideally I'd need a solution that is do-able inside Graphical Modeler, and doesn't require any extra plugins etc.
I've looked at some of the other answers but nothing has solved this for me.
The flight path loops back on itself multiple times (i.e. it's not a straight line) so I can't use nearest neighbour.


texts: Where does qgis store calculated field from field calculator?

texts:
When adding a (computed) field for a PostGIS table, where is the new field stored?

I thought it would be with the original table, but it does not seem to be the case for QGIS 3.0.1/3.0.2 (Ubuntu 16.0.4). When I create a calculated column label by copying from a column called name, the label column is displayed in QGIS but does not show in PostgreSQL (10.3) /PostGIS (2.4).

I was wondering where is the new column stored?

The content of table originally was:

=> select * from  tt;
 ogc_fid | name |                        geom
---------+------+----------------------------------------------------
       1 | i    | 0101000020416900000000000000000040000000000000F03F
       2 | i    | 01010000204169000000000000000000400000000000000040
(2 rows)

After the field calculation, it stays the same (without the new column).

The definition of the table as reported by pg_dump is:

--
-- PostgreSQL database dump
--

-- Dumped from database version 10.3 (Ubuntu 10.3-1.pgdg16.04+1)
-- Dumped by pg_dump version 10.3 (Ubuntu 10.3-1.pgdg16.04+1)

SET statement_timeout = 0;
SET lock_timeout = 0;
SET idle_in_transaction_session_timeout = 0;
SET client_encoding = 'UTF8';
SET standard_conforming_strings = on;
SELECT pg_catalog.set_config('search_path', '', false);
SET check_function_bodies = false;
SET client_min_messages = warning;
SET row_security = off;

SET default_tablespace = '';

SET default_with_oids = false;

--
-- Name: tt; Type: TABLE; Schema: public; Owner: db
--

CREATE TABLE public.tt (
    ogc_fid integer,
    name character varying(10),
    geom public.geometry(Point,26945)
);

ALTER TABLE public.tt OWNER TO db;

--
-- PostgreSQL database dump complete
--


texts: Scene boundary issue when reducing sentinel 1 in natural values vs dB

texts:
I have written a function that smooths sentinel 1 images by reducing across clusters created from a sentinel 2 median image, and then reduces the sentinel 1 collection to the p5 value. When I run the function on dB values it works as expected, but when I use natural values strange things happen at the scene boundaries.
An example is here:
// The location of the scene boundary I am looking at
var geometry = ee.Geometry.Polygon([[25.356019971276755,-15.045220885450785],[25.450433728600974,-15.045220885450785], [25.450433728600974,-15.014549928161024], [25.356019971276755,-15.014549928161024], [25.356019971276755,-15.045220885450785]]);

Map.centerObject(geometry,14);

// Cloud mask the S2 collection and reduce to median
var s2image = ee.ImageCollection('COPERNICUS/S2_SR').filterBounds(geometry)
.filterDate('2020-03-01','2020-05-01').filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 25)
.map(function (img){
  var cloudProb = img.select('MSK_CLDPRB');
  var scl = img.select('SCL');
  var cloud = cloudProb.gt(0).or(scl.eq(3));
  var mask = cloud.lt(1);
  return img.select(['B4','B3','B2']).updateMask(mask);
}).median().divide(1000);

Map.addLayer(s2image,{max:2},'S2 median',0);

// Compute clusters from the s2 median image
var clusters = ee.Algorithms.Image.Segmentation.SNIC(s2image,25,3).select('clusters');

Map.addLayer(clusters.randomVisualizer(),{},'Clusters',0);

// Filter the s1 collection and convert to natural values
var s1coll_dB = ee.ImageCollection('COPERNICUS/S1_GRD').filterBounds(geometry).filterDate('2019-03-01','2020-02-28').select('VH');
var s1coll_nat = s1coll_dB.map(function(image){return ee.Image(10.0).pow(image.divide(10.0)).rename('VH')});

// Function to compute the cluster VH means and reduce the collection to p5 value
function  reduceS1(coll){
  var means = coll.map(function getMeans(image){
    var clusterMeans = clusters.addBands(image).reduceConnectedComponents(ee.Reducer.mean());
    return image.addBands(clusterMeans,['VH'],true);

  });
  return means.reduce(ee.Reducer.percentile([5]));

}

Map.addLayer(reduceS1(s1coll_dB),{min:-25,max:-10},'dB');
Map.addLayer(reduceS1(s1coll_nat),{max:0.02},'natural values');

https://code.earthengine.google.com/858b09f6a49e36c00de39d676f7e0d61


texts:
I'm building a runtime layer and til a couple of days ago my code was working. Since today it started giving me an exception that is causing subsequent code not to work.

Uncaught ReferenceError: layerView is not defined

The problem is caused by this call:

cdb.vis.Vis.addInfowindow(map, sublayer1, ['cartodb_id','ppe_ma']);

I guess something has changed in the api but i dunno how to track this changes as these functions are not yet documented.


texts: Retrieving attributes from multiple polygons intersecting a given zone (QGIS)

texts:
I am trying to find a way to display relevant information for a given zone based on all the layers that it intersects.
For example, I'm studying five zones and I'd like to know for each of them which type of soil can be found within their boundaries.
I have tried to "join features by location" but I end up with multiple features when different type of soils are present in the area. This is not ideal because I'll have to repeat the process for many other layers (soil, geology, hydrography, ...)
In the image below for example, I would need to know which colors can be found within the red dotted area (yellow, blue, dark green and light green). It doesn't necessarily need to be stored as attributes, if I can have a .txt file with the info that's also fine.


texts: How can I batch export all bookmarks to Illustrator?

texts:
My map has about 40 bookmarks. I want to export each to an individual Illustrator file. I've tried to adapt code I saw elsewhere, but I'm a total python novice. Instead of getting 40 exported files, it's re-exporting the same file name over and over.

import arcpy.mapping
Folder = r"H:\My Documents\Maps"
mxd = arcpy.mapping.MapDocument("CURRENT")
df= arcpy.mapping.ListDataFrames(mxd, "Layers") [0]

for bkmk in arcpy.mapping.ListBookmarks(mxd, "", data_frame=df):
    df.extent = bkmk.extent
    dfname = bkmk.name
    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToAI(mxd, Folder + "\map" + dfname + ".ai", resolution=300)

Edit November 15, 2017: this code now seems to work, at least once. I'm not sure what problem I had going before; maybe it was related to the output folder. I have found, however, that if something breaks in the script (e.g., one time I had a backslash in a bookmark name) I can't rerun it even after fixing the issue. I have to restart my ArcMap session.


texts: Using Workspace To Raster Dataset tool in ArcGIS Desktop?

texts:
I am going to use the "Workspace to raster dataset" tool to merge several tif files. I have 4 test images (same area, same geographic, just different in time).

The input workspace has been set up, it is a folder with four images.
The target Leicester dataset was used after copying and renaming one of the four images.

I set "Mosaic" to "MEAN" and "Mosaic colormap Mode" to "Match" among the options. In this setting, I think that output Leicester's one arbitrary grid should be the average value on the same grid as the four images I entered. However, output was not the average of the four images entered. I've tried everything, changing all the options.

I knew this would create Leicester with the average value of the Leicester I entered. What is the problem?


texts: Concatenating part of string and sorting it using QGIS

texts:
I'd like to know how to concatenate part of a string and sort it.
For example, i have a FIELD1 (string) and i should have the numbers of the string separated with ';' in the FIELD2. Only the numbers are repeated in the records.
Moreover, i should sort the numbers ascending about their values.

How do I create an array string, in the field calculator of QGIS, to do this?


texts: Using Expressions to underline part of a label

texts:
I am using QGIS (sadly only 2.8) and I have a label on my dataset. It currently looks like this:

    "MP_Lookup_Constituency" || '\n' ||  "MP_Lookup_Name" +' ('+"MP_Lookup_Party"+')'

The Constituency appears above the rest of the text, but I'd like to (not sure if possible) just underline this part of the label. So the Constituency is underlined, but the rest isn't. I've tinkered with different things, but no joy, so any help would be amazing!


texts: Fixing address errors using Python Parser of ArcGIS Field Calculator?

texts:
I am working on a task that will calculate a field to fix errors stemming from the Standardize Addresses tool and human error with the input data.

The script is below. I'm getting Error 000989 "Python syntax error: Parsing error SyntaxError: invalid syntax (Line 6)"

Note: ADDR_SN and ADDR_SD seems to be a truncated output of Street Name and Street Direction (There is some confusion as the Standardize Address tool is new to me)

Any pointers to help fix this code?

Expression:

AddSN2( !ADDR_SN!, !ADDR_SD!)

ADDR_SN is the Street Name and ADDR_SD is the street direction.

Code Block:

def AddSN2(StNa, STDir):
    if (StNa == "BEND"):
        return (STDir+ +StNa)
    elif (StNa == "CEDAR GARDEN"):
        return (StNa.replace(" ", "")
    else:
        return StNa


texts: How to improve cartographic visualization of overlapping polygons?

texts:
I need some advice about how to improve the visualization of cartographic information.

User can select different species and the webmapping app shows its geographical distribution (polygonal degree cells), each specie with a range of color (e.g darker orange color where we find more info, lighter orange where less info).

The problem is when more than one specie overlaps. What I am currently doing is just to calculate the additive color mix of two colors using http://www.xarg.org/project/jquery-color-plugin-xcolor/

As you can see in this image, the resulting color where two species overlap (mixed blue and yellow) is not intuitive at all.

Someone has any idea or knows similar tools where to get inspiration? for creating the polygons I use d3.js, so if more complex SVG features have to be created I can give a try.

Some ideas I had are...

The more data on a polygon, the thicker the border (or each part of the border with its corresponding color)
add a label at the center of polygon saying how many species overlap.
Divide polygon in different parts, each one with corresponding species color.


texts: Newly added layers to QGIS 2.18 will not display in the print composer

texts:
I have an existing QGIS project that until today was just fine.  I added a new shapefile to the project but these new data fail to display in the print composer.  All other existing data in the project display just fine as do the other map elements in the print composer.  I have made no other alterations to the project or the print composer layout other than adding these data.

I have created a new blank project, added these data in question to the new project, created a new print composer, and confirmed that the data draw in this new project.  This leads me to believe this particular dataset in not the issue.

I have duplicated the print composer layout and this new layout has the same issue.

I have added other data for this area to the project and cannot view these data as well leading me to believe that there is some corruption in the print composer.

Please note that all data share the same CRS, and there were no alteration to any settings or elements in either the main display or the print composer.  This is happening on QGIS 2.18 with Ubuntu 16.04

What steps can I take to try to get my data to display in the print composer?


texts: What PostgreSQL version to use with QGIS 2.8 (on Linux/Fedora 21)

texts:
I'm looking to install QGIS on my Linux/Fedora 21 laptop.

QGIS ver. 2.8 seems to be the new/stable version (Feb. 2016)

What PostgreSQL version should I use with that?


texts: How to extract the maximum value from an array in the QGIS Field Calculator?

texts:
I'm trying to get the maximum value from an array, using the max() function. The syntax for the max() function is this:

  max(value1,value2…)

I have an array: [ 1, 2, 3, 4, 5, 6, 7, 8 ] which I got by using this expression:

array_remove_all( array_agg( "fieldname"), maximum("fieldname"))

I need to obtain the maximum value from the array, but the max() function doesn't  work. When I calculate the max() of the array, it gives this error:

  Eval Error: Cannot convert '' to double

I can convert the array to a string using the array_to_string() function. When I calculate the max of the string version of the array, it gives this error:

  Eval Error: Cannot convert '1,2,3,4,5,6,7,8' to double

So it seems the max() function can't accept an array or a list in string format as input. I also tried the maximum() function, with similar lack of results. I don't see any other function in the Field Calculator for finding the maximum value

How to extract the maximum value from an array in the Field Calculator?


texts: Including .prj file for new CRS in QGIS?

texts:
I have a coordinate system in coal grid which is exclusively for coalfields of India in QGIS.

I want to transfer it from the coal grid to UTM.

I have .prj file of coal grid, But I don't know how to include it. How to set custom project system?


texts: Defining Projection Lambert Conformal Conical w custom GCS

texts:
I have a shapefile from Mexico´s INEGI site that doesn´t have a .prj file.  I am having trouble defining the projection in ArcGIS.  INEGI tells me it is projected as Lambert Conformal Conical with a ITRF92 GCS.  From what I understand, this GCS should be the same as WGS84 (I have used this with other INEGI files and it matched).  This combination of LCC & WGS84 doesn't exist in ESRI´s system, so I tried to create a custom projection coordinate system with the following information given by INEGI:

longitude of origin: -102.0,
latitude of origin: 12.0,
standard parallel 1: 17.5,
standard parallel 2: 29.5,
false easting: 2500000.0 m,
false northing: 0.0 m

But this doesn´t work - ArcGIS puts the file in Africa, not Mexico (though it seems to be the right scale).

I´m wondering if there is something going wrong with the GCS, for which INEGI gives this information:

Datum: ITRF92
Sphereoid: GRS80
Semimajor axis: 6378206.4
Inverse Flattening 294.9786982

These numbers are not the same as ESRI´s WGS84 or ITRF92, and I can´t figure out how to change them in the "Define Projection" tool.

This is the link to the site that contains the data: http://www.inegi.org.mx/geo/contenidos/recnat/edafologia/vectorial_serieii.aspx

At the bottom are the links to the date, I downloaded "Perfiles de Suelo".

The extent of the data is:

Top -2347872.613848
Bottom -340538.708859
Left -1075400.666415
Right  - 4073175.438209

The data is supposed to cover the nation of Mexico.


texts: Set CRS for Raster in QgsRasterLayer (pyqgis)

texts:
in a QGIS python plugin I get some code that writes a asc file, this is the header of the final file:

NCOLS  20
NROWS  15
XLLCORNER  664400.000000
YLLCORNER  4749500.000000
CELLSIZE   100.000000
NODATA_VALUE   -9999.0
NODATA_VALUE   -9999.0

as you can see I get the X and Y corners of the file.

I wanted also to load the file in QGIS and obviously it asks me for the CRS.

What I am looking for, is to avoid this message and the related pop-up window while using QgsRasterLayer.

I've found this 2 old questions:

How can I specify the CRS of a raster layer in pyqgis?

and

How to get CRS of a raster layer in pyqgis?

but I cannot use them in order to solve my problem.

Someone has an idea?


texts: QGIS asks for PostgreSQL credentials multiple times

texts: I have several QGIS documents that pull layers and tables from PostgreSQL / PostGIS databases.  When I open the documents QGIS asks for the the credentials several (5-8+) times.  Each document pulls its data from a single database, the credentials given have the same level of access to all layers and tables in the database.  Is there a way to get the number of requests down to one per session?


texts: Distance based random point selection in r?

texts:
To avoid clustering in subsequent analysis I try to randomly select points of a SpatialPointsDataFram with a given distance between the points. The remaining part should be devided for test analysis. I found Randomly sampling points in R with minimum distance constraint? from which I think it leads the right way, but I am not completely sure about the code.

Maybe there is function, like sample, that already includes a distance parameter?


texts: Save Rule Based Symbology as Default Layer Style

texts:
I am constantly moving layers in and out of my projects, I was able to set default style for layers and they all the layers I drag and drop to my project get the same style by default.

Lately I have discovered that I can create rule based symbology as a layer style which allows me to hide unselected features while I'm working on a selected feature.

I am wondering if it's possible to save rule based symbology as default for imported vector layers? When I use save symbol as default it saves a QGIS configuration file next to me shapefile but that's good for only one layer, as soon as I change the project I need to set everything up from scratch again by going to rule based symbology, create two rules, fill them in with the expression etc.


texts: ArcGIS Server running in different clusters?

texts:
I´m looking at different configurations of ArcGIS Server 10.2.1 and I can´t understand the installation of ArcGIS Server running in diferent clusters.

Is it only a configuration setting that use the resource of the machine or is it necessary for different machines servers with one license for each machine?

With the answer of Vince, is this correct?


texts: Why is Geoserver downsampling my geotiff?

texts: I created a KML file from a geotiff using Geoserver, but the image looks really pixelated in Google Earth. The image is a geotiff and is only around 7MB. Any reason why it might look so pixelated?


texts: Setting up files for Spline using ArcGIS Pro Spatial Analyst

texts:
I have converted Polygons to Raster so I can perform a spline, but I cannot get the raster as in input on Spline without Error 000840. Tracking down other solutions, I tried to "Make Feature Layer" from a .lyr file I made for the Raster, and that too has a Error 000840.
Any ideas on what I'm missing to get it all to fit in the Spline tool?


texts:
I'm a GIS professional. One of my clients (an ArcMap user) has asked our firm to teach some of his staff about ArcGIS online. In so doing, we are to provide an overview of GIS-- how to think spatially, how to visualize data, etc.

My first instinct when teaching "GIS 101" is joining data. You (the client/user) has a dataset (i.e. median income by zipcode), and you want to make a map. Is this possible using only ArcGIS Online?

I think this is one of the most exciting aspects of learning GIS, and I can think of a million ways to achieve this using ArcMap, python, geojson, etc., but I can't seem to figure out an easy to do this using ArcGIS Online. What are my options?


texts: Appending from file geodatabase to ArcSDE fails using ArcPy?

texts:
I'm trying to append a feature class from a file geodatabase to an Arcsde featureclass using this line of code:

arcpy.Append_management("parcel_copy", sde_path+"parcel" ,"NO_TEST")

This line works under Python console in ArcMap , but when I execute it under Pyscripter IDE, all the code works fine except the append statement, and it does not return an error.

I have tried to isolate the append statements in a little script, and it worked, but i still don't understand why it has not worked.

i wonder if i have not deleted a featurelayer refering to a featureclass in SDE, that locks this featureclass, so i cannot append to it.

any explanation?


texts: Retaining timestamp when saving gpx as shp in QGIS?

texts:
I'm somewhat new to QGIS.

I used my Garmin GPS unit to record my tracks while collecting plant specimens for a plant inventory. When I collected each plant I noted the time of collection. I want to get lat/long info for each collection time.

But when I import the GPX file and then save as an .shp file so that I can add lat/long columns, the "time" column is reduced to only the "date" - the hhmm has been removed.

I used XY tools to export the gpx and shp files to get a look at the time column.

The gpx timestamp format looks like this: >>PyQt4.QtCore.QDateTime(2016, 3, 19, 17, 33, 55)

The shp file, looks like this: >>PyQt4.QtCore.QDate(2016, 3, 19)

I'd like to not have to export files to .xls and combine from there, is there a way to retain the original timestamp when converting from GPX to SHP?


texts: Getting Google satellite images (via plugin) to behave well in QGIS with map created with different CRS

texts:
How to get Google satellite images (via the plugin) to behave well in QGIS with a map created with a difference CRS?

I have a problem using the Google Image layer plug in - at some scales, it won’t match the rest of the project (but does at others).

The project started with Tiff files imported from GB Ordinance Survey and set the CRS to EPSG:27700 - OSGB 1936 / British National Grid (which I think these tiles default to).

I then added a couple of shapefile layers (set to the same CRS). One has point data, the other is a polygon which defines a site boundary.

I then added (via the plug in) the google satellite image layer.

Whenever I add this layer to any project, this always defaults to a view where the whole world is repeated several times across the screen.  If I zoom in by one click from there, it seems to find itself ok and the site boundary matches the OS tile and the satellite image.

However on this project (and a previous one), when I zoom in too far, suddenly the site boundary doesn't match the same area on the google image (about 1:1000) - it is quite a jump.  I can zoom back out any amount and the polygon and google image match fine.

When the polygon and google image don’t match up, if I bring the OS tile above the Google image, I can see that the polygon and OS tile have remained lined-up.

I know Google use a different CRS - WGS 84 / Pseudo Mercator, which I assume is the root of the problem, but there must be a fix.  I’ve tried changing the google image layer’s CRS to match the rest of the project with the obvious result that it goes even further off (“sorry we don’t have image for here”) and right clicking on the OS layer and selecting “set project CRS from layer” doesn't help.


texts: Creating a new layer for free space between existing layers

texts:
How can I create a new layer that will only show the 'free space' between the raster and vector layers that are currently on my map?

Will I have to merge my existing raster and vector layers and use some kind of inverse selection tool? If so could someone explain how?

Any ideas?


texts:
My Google Maps is in English, and when I want to a share a location it gives me the address in English.

But how can I get this address string in another language?
I can't find a way to translate it.


texts: Exporting multiple tables to Drive Earth Engine

texts:
I have a list of FeatureCollections that I want to export all at once (instead of writing out the export module multiple times) using Javascript API. I have seen the posts by @Rodrigo E Principe such as here and others like this one showing the use of normal JS loops (admittedly the second link is for exporting to asset).

I am trying to do the exact same thing as the first link, just exporting a FeatureCollection instead of an ImageCollection. However, when I run the code I get no tasks to run nor am I getting any errors. I don't understand what's happening.

Here is my code:

// Make a list of Features.
var features = [
  ee.Feature(ee.Geometry.Point(97.58, 24.36), {name: 'sample1'}),
  ee.Feature(ee.Geometry.Point(99.45, 20.26), {name: 'sample2'}),
  ee.Feature(ee.Geometry.Point(93.24, 21.42), {name: 'sample3'}),
  ee.Feature(ee.Geometry.Point(96.13, 24.12), {name: 'sample4'}),
  ee.Feature(ee.Geometry.Point(96.88, 22.48), {name: 'sample5'})
];

// Create a FeatureCollection from the list and print it.
var points = ee.FeatureCollection(features);
print(points);

var l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate('2015-06-23', '2020-01-31')
                  .filterBounds(points);
print("Landsat", l8sr)

var index = ee.List.sequence(0,4)
var coords = points.geometry().coordinates()
print('coords', coords)

// do this for each index
var test1 = index.map(function(q){
    var sub = l8sr.filterBounds(ee.Geometry.Point(coords.get(q)))
    return sub
});

//run test for first one
var firstcoll = ee.ImageCollection(test1.get(0))
var img = firstcoll.first();
var bands = img.bandNames()

//get the value for each band for every image
var loop = index.map(function(ind){
  var getvals = function(img1){
    var date = img1.date().format('yyyy-MM-dd')
    var value = img1.select(bands)
   .reduceRegion(ee.Reducer.first(), ee.Geometry.Point(coords.get(ind)))
   return ee.Feature(ee.Geometry.Point(coords.get(ind)), {value: value, date: date})
  };
var newft = firstcoll.map(getvals);
return(newft);
});

print('loopresult', loop)
print('loop ex', ee.FeatureCollection(loop.get(0)))

// loop on client side
for (var i = 0; i<loop.length;i++) {
  var toexport = ee.FeatureCollection(loop.get(i));
  var nam = "test_point_" + i

      // Export
Export.table.toDrive(
  {collection: toexport,
  description: nam,
  folder: "dissertation",
  fileNamePrefix: nam,
  selectors: ['B1','B2','B3','B4','B5','B6','B7','B10','B11',
  'pixel_qa', 'radsat_qa', 'sr_aerosol']
  });
}


texts: Distances between points British National Grid (BNG)

texts: I have a point dataset expressed in British National Grid (about 10,000 points distributed all over the UK). I want to do some cluster analysis and I am aware that distance distortion might have an effect on the results. The precision I need is quite coarse (about 200m-500m). Do I need to take steps in the analysis or can I just calculate Euclidean distances?


texts: Efficiently dividing KML file by area, while using ogr2ogr?

texts:
I have a work process which uses an ogr2ogr text command to convert a MapInfo TAB file to a KML file. I'm very pleased with this because it maintains the original symbology (a simple 2 colour system used for lines in MapInfo - implemented as a line style when input, not derived directly from the data).

This command also selects for certain categories of line according to some data fields in the original file.

ogr2ogr -f "KML" OutputFileName.kml InputFileName.TAB -sql "select RNumber from InputFileName where ((RType = 1 or RType = 2) and OStatus >= 3)"

Unfortunately the original file is large, and the KML file therefore is large - too large for the application I'm using it in.

I'd like to do something simple and efficient to allow me to divide the KML file up into geographical regions. If possible this should be something simple enough to describe to others too - although that's not a necessity.

I could, of course, work with a process in MapInfo to extract the lines I want first - or I could do the same with QGIS. I've tried both and hit problems - which I'm sure could be overcome - but I'd like to know if I'm missing something more obvious/simple. For instance does ogr2ogr have the capacity to do this?

As background, I've faced an issue with QGIS which is that bringing the TAB file into the application seems to lose the line colour. I could re-create this based on a data field if necessary, but that adds extra steps (I know that I could use various workspaces/projects to simplify things). I'm also not sure how to ensure that QGIS codes line colour into the KML.

I also find the processes of selecting by geographical area in MapInfo frustrating - I'll be able to solve this if necessary I think, but I end up in a tangle of SQL query issues which so far hasn't given me what I need.


texts: When does a Python script in QGIS execute drawing of QgsRubberBand?

texts:
I have a script which has an unrelated error after the draw_polygon function. I don't care about the error at the moment because I want to see the current results. For some reason my draw function doesn't work. My thought was that the drawing of QgsRubberBand happens after the succesful completion of a function?

def draw_polygon(Polypoints):
    point = QgsVertexMarker(canvas)
    poly = QgsRubberBand(canvas, True)
    poly_points = []
    for feature in Polypoints:
        point.setCenter(QgsPoint(feature[0], feature[1]))
        point.setColor(QColor(255, 150, 0))
        point.setIconSize(20)
        poly_points.append(QgsPoint(feature[0], feature[1]))
    poly.setToGeometry(QgsGeometry.fromPolygon([poly_points]), None)
    poly.setColor(QColor(0, 150, 150))
    poly.setWidth(3)

I am using QGIS 2.18.


texts:
I have a geotiff file served with geoserver. The geotiff is styled with the following sld, which includes a spatial filter. I am sending a polygon geometry as an env parameter, which results in showing the part of the geotiff inside the polygon area.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld"
xmlns:ogc="http://www.opengis.net/ogc"
xmlns:xlink="http://www.w3.org/1999/xlink"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://www.opengis.net/sld
http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
version="1.0.0">
  <NamedLayer>
    <Name>raster_layer</Name>
    <UserStyle>
      <FeatureTypeStyle>
        <Transformation>
          <ogc:Function name="gs:CropCoverage">
            <ogc:Function name="parameter">
              <ogc:Literal>coverage</ogc:Literal>
            </ogc:Function>
            <ogc:Function name="parameter">
              <ogc:Literal>cropShape</ogc:Literal>
              <ogc:Function name="env">
                <ogc:Literal>shape</ogc:Literal>
              </ogc:Function>
            </ogc:Function>
          </ogc:Function>
        </Transformation>
        <Rule>
          <RasterSymbolizer>
            <Opacity>1.0</Opacity>
            <ColorMap type="values">
                <ColorMapEntry color="#e6004d" quantity="1" label="1" opacity="1"/>
                <ColorMapEntry color="#ff0000" quantity="2" label="2" opacity="1"/>
             </ColorMap>
          </RasterSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

I request the wms layer using the following code:

L.tileLayer.wms(WMS_ROOT_URL, {
    layers: layername,
    format: 'image/png',
    transparent: true,
    version: '1.3.0',
    styles: myStyle,
    env: 'shape:' + clipPolygon
});

While this works fine with smaller polygons, the polygons can get very large and complex, resulting in thousands of characters, which gives a 414 error: Request-URI Too Long.

I am trying to find a way to fix this, but I had no success so far. What I was thinking as possible solutions were:
1. Use post in wms request
2. Send the id of the polygon instead of its geometry and somehow use this to get the polygon geometry in the sld
3. Other solution?

I am using:
Geoserver v.2.11.1 and Leaflet v.1.2.0


texts: How do I get basic statistics on a field in a joined table?

texts: I have a shapefile to which I've joined a table, matching with an ID field. All the fields of the table show up correctly when I open the shapefile's attribute table. But when I go into vector/analysis tools/basic statistics I can only choose one of the fields in the shapefile itself; those I have joined are missing.


texts: QGIS Legend multiple name for one shapefile using styles

texts:
I am using styles and presets in QGIS to display the same shapefile on various composers.
The shapefile is representing townships with several fields such as population density, length of water supply pipelines per township, etc.

Each composer is displaying the shapefile with categories based on the fields mentioned above.

I would like to automatically change the name of the shapefile in the legend of the composer by using presets. Is it possible or do I have to manually change the name in each composer?

To do it manually, I use the following feature in Legend options

In short: I would like to have my layer name "displayed as" XXX for style 1 and YYYY for style 2. This would affect presets, styles and legend.


texts: Attribute table empty with custom driver in newer versions

texts:
I've written a custom ogr driver for one of our own formats.
When loading the data I get to see vector data nicely displayed in Brighton/Essen/Pisa. However, when opening the 'Attribute Table' only
Brighton is capable of showing the data. Using Identify on a feature I do get it's correct attributes in all 3 versions. All 3 versions also display the correct amount of features.

The attribute table shows 0 Features in Essen and Pisa.
On the other hand if i load a shp, and check the contents of the Attribute table
I do get to see all the attributes in all 3 versions.

This leads me to believe that somehow the attribute table does not get updated upon new feature request to the driver (or something of the sort). The Driver is an RTree implementation and will not return features unless they are inside the provided boundingbox / Envelope.

Any hints?

p.s using ogrinfo i'm also seeing al the correct attributes just like with identify.


texts: Two layers with same SRs and same coordinates don't overlap

texts:
i'm trying to do some analysis on road accidents.

I have :

a .csv with lat,long data of road accidents
a .shp with the road
network

When I create a new layer from .csv, it is automatically assigned to WGS 84 epsg 4326.
The road networks instead is  WGS 84 32N epsg 32632

The layers overlap when using the different SRs. The problem is when I change them to a uniform SR (that is, when creating a new shp of road accidents with epsg 32632), they don't overlap anymore and can only be visualized one at time.

Any ideas?

I'm using Qgis 2.18 (but the same with Qgis 3.0)


texts: I added layers from Geoserver, I have line objects. But when I open style editor in uDig, I have only WMS Named Style. I thought I could change for example colors by attributes, but I don't how to do it when it offers me only WMS Named Style. I'm working with uDig and Geoserver for the first time so please help me.


texts: is there a open source way to georeference a non geospatial vector (like eps, svg etc)?

texts:
Is it possible to reference a non geospatial vector file such as eps, pdf, svg, dxf etc.?
The vectors I'm talking about are points and lines without any georeferences, but the size, the position and the CRS of their bounding box is known.

Since there are some thousands of vectors, the rasterize -> georeference -> vectorize way is no solution.


texts: Replace the data source for several layers in several mxds - from shapefiles to the ARCSDE data source

texts:
Possible Duplicate:
  Permanently replace/update SDE server + service in each SDE Layer in an MXD

First let me say I have probably read every link regarding changing DATA Sources on the internet from the ARCGIS 10 tool to Python, however none have worked for me.  What i'm trying to do is very simple in concept.

Please don't recommend the ESRI Replacing data source help.

I'm updating to 9.3 to 10.1 - I need to replace my shapefile data source to my ArcSDE files.
I have over 4000 mxds with several shapefiles with the data source pointing to: R:\Shapefiles[layer folders*.shp
I need to change the shapefile datasource to point to C:\Users\lmuhammad\AppData\Roaming\ESRI\Desktop10.1\ArcCatalog\Connection to PGCGIS02.sde[gisdev.layername][gisdev.filename]
I need to repeat for over 4000 mxds automatically.

easy right.. HELP!!!!

------ UPDATE

This works ok xcept if the shapefile name differs greatly from SDE featureclass, =no connection any fixes?

import arcpy, glob

newSdePath = r"Database Connections\Connection to sdename.sde"

for file in glob.glob(arcpy.GetParameterAsText(0) + r"C:\*.mxd"):

mxd = arcpy.mapping.MapDocument(file)

for df in arcpy.mapping.ListDataFrames(mxd):

    for lyr in arcpy.mapping.ListLayers(mxd, ""):

        if lyr.isFeatureLayer== True:

            lyr.replaceDataSource(newSdePath, "SDE_WORKSPACE","", False)

mxd.save()del mxd


texts: Different scale for x and y in a QGIS map layout

texts:
I'm making some geological cross sections with the geoscience plugin. I have been able to do this and I can make a print layout where I have my cross section (which the plugin makes in x-y space where y is really z) and the plan view on top and it looks great. Now what I can't figure out if there is anyway to make a vertical exaggeration on the cross section (I don't generally like doing that because all the angles end up wrong, but still need to). I cannot find a reference to different scaling on the x and y for qgis for a layout, and the only thing I have come up with is to define a custom coordinate system where the y (which here is really z) is scaled by a certain amount relative to x.
So my question is, is there an easier way to do this or is there a simple way to scale a custom coordinate system (I'm trying to get the affine to work)
Thanks


texts: IApplication reference with ArcObjects Class Extension

texts:
I have an ArcObjects (9.3) Class Extension which has an event handler for OnCreate so that feature properties can be manipulated / assigned when a feature is first created.

Some of the functionality needs enclosing in an edit operation, which requires an IWorkspaceEdit object. As far as I know the only way I can get the IWorkspaceEdit is by getting the IEditor (extension) from IApplication, then grabbing the IEditor's EditWorkspace.

In an ArcObjects command (extending BaseCommand) I get the IApplication 'hook' passed on create (that is, the command's creation), so I can retain a reference to it. However, in the class extension's OnCreate method I am only passed the IObject which has been created. I can't see a way to get the IObject's workspace to work from that direction.

Can anybody please tell me

How to access the calling IApplication from a class extension or,
How to get an IObject's EditWorkspace or,
Any other way to create an edit operation within a class extension's OnCreate method?

Any help very much appreciated.


texts: ArcObjects Style of a symbol chosen with the ISymbolSelector interface

texts:
I am trying to store the properties of a SimpleMarkerSymbol symbol selected by the user using the ISymbolSelector interface.

I plan to store the size, colour and style number in a file.

I can get the size and colour of the selected symbol but cannot locate the style.

If the user selects one of the Simple Marker Symbols without changing the size or colour then my code works. If they change anything, the code cannot find the style in the style gallery.

Even though I modify the Style gallery enumeration so that they have the same size and colour as the chosen symbol the clone test fails.

Is there any other way of doing this?

Here is my code:

Dim sError As String
Dim pSymbolSelector As ISymbolSelector
Dim pMarker As ISimpleMarkerSymbol
Dim pColour As IColor
Dim pCol As System.Drawing.Color
Dim iColour As Integer
Dim sName As String = ""
Dim dSize As Double
Dim sStartEnd As String = comboStEnd.Text
Dim iStyle As Integer = -1
Dim pSymbol As ISymbol
Dim pMSym As IMarkerSymbol
Dim pStyleStorage As IStyleGalleryStorage
Dim pStyleGallery As IStyleGallery
Dim pMxDoc As IMxDocument
Dim pEnumStyleGall As IEnumStyleGalleryItem
Dim pStyleItem As IStyleGalleryItem
Dim pMarkerSym As IMarkerSymbol
Dim pMSymClone As IClone
Dim pMarkerClone As IClone

Try
    sError = ""
    StatusLabel1.Text = ""

    If String.IsNullOrEmpty(sStartEnd) Then
        StatusLabel1.Text = "Please choose Start or End"
        Exit Sub
    End If

    pSymbolSelector = New SymbolSelector

    pMarker = New SimpleMarkerSymbol
    If Not pSymbolSelector.AddSymbol(pMarker) Then
        Exit Sub
    ElseIf pSymbolSelector.SelectSymbol(0) Then
        pSymbol = pSymbolSelector.GetSymbolAt(0)
        pMSym = TryCast(pSymbol, IMarkerSymbol)
        If pMSym Is Nothing Then Exit Sub

        dSize = pMSym.Size
        pColour = pMSym.Color
        pCol = New System.Drawing.Color
        pCol = modMain.colourESRItoMS(pColour, sError)
        If sError.Length > 0 Then
            sError = "cmdChoose;" & vbCrLf & sError
            MsgBox(sError)
            Exit Sub
        End If
        iColour = pCol.ToArgb

        pMSymClone = CType(pMSym, IClone)

        ' Get the symbol style
        pMxDoc = m_pApp.Document
        pStyleGallery = pMxDoc.StyleGallery
        pStyleStorage = pStyleGallery
        pEnumStyleGall = pStyleGallery.Items("Marker Symbols", "ESRI.Style", "Default")
        pEnumStyleGall.Reset()
        pStyleItem = pEnumStyleGall.Next
        While pStyleItem IsNot Nothing
            pMarkerSym = pStyleItem.Item
            pMarkerSym.Size = dSize
            pMarkerSym.Color = pColour
            pMarkerClone = CType(pMarkerSym, IClone)
            If pMSymClone.IsEqual(pMarkerClone) Then ' True if pMSym has not been modifed for size or colour
                sName = pStyleItem.Name
                iStyle = pStyleItem.ID
                Exit While
            End If
            pStyleItem = pEnumStyleGall.Next
        End While

    End If
Catch ex As Exception
    sError = "ERROR: cmdChoose; Line: " & Erl() & vbCrLf & ex.ToString
    MsgBox(sError)
End Try


texts: Difference between Layer Definition File (qlr) and QGIS Layer Style File (qml)?

texts:
What is the difference between Layer Definition File (qlr) and QGIS Layer Style File (qml)? Which one is recommended?

I have only used qlr so far. It seems store all the data displayed in the "properties" of the layer.


texts: Can I use a startup.py in QGIS which is stored outside the user folder?

texts:
I would like to use a startup.py that is not stored in the default QGIS-user-folder like ~/.qgis2 or a folder used with the configpath option.

The startup.py I want to use should be located in a folder that is write protected for the users and should only be edited by sysadmins.

Boundless has implemented a interesting functionality called "QGIS Initialization scripts" ( https://connect.boundlessgeo.com/docs/desktop/latest/system_admins/init_scripts.html ).

As a workaround I check in the QGIS startup-file (e.g. qgis.bat) if "my" startup.py is located in the users .qgis2\python-folder. If there is no startup.py in the users folder I copy it into this folder. Otherwise I do a md5 check if the startup.py is identical with my startup.py and overwrite the users startup.py if it is not identical:

:start_qgis

set "file=%QGIS_UDIR%\python\startup.py"
set "referenzfile=%OSGEO4W_ROOT%\bwtacheck\startup.py"

if exist "%file%" goto start_bwta_check
copy /Y %referenzfile% %QGIS_UDIR%\python\startup.py
goto start_qgis2

:start_bwta_check
call %OSGEO4W_ROOT%\bwtacheck\md5.bat %file% file_md5
call %OSGEO4W_ROOT%\bwtacheck\md5.bat %referenzfile% referenz_md5

if "%file_md5%" equ "%referenz_md5%" (
echo MD5 identisch
) else (
echo MD5 nicht identisch
copy /Y %QGIS_UDIR%\python\startup.py %QGIS_UDIR%\python\startup_overwritten.py
copy /Y %referenzfile% %QGIS_UDIR%\python\startup.py
)

:start_qgis2
REM Run the central update commmandfile if the variable is set to a location and the file exist.
if *%QGIS_NETCMD%==* exit /b
if exist %QGIS_NETCMD% call %QGIS_NETCMD%
exit /b

Does the plain QGIS also support the use of a startup-script that is stored outside the typical userprofile-folder (the one used by qgis to store the plugins etc.)?

I tried to use a startup.py in the program-folder (like C:\OSGeo4W64\apps\qgis-ltr\python) but this did not work. Is this a bug or meant to be?

PS: It is important that the startup-script gets executed after QGIS has initialised.


texts: Using SearchCursor to access & export values in a raster table

texts:
I am desperately trying to export the values in a raster table to a .txt file. I've seen some other folk post about this, but for some reason can't get the code to work in my case.

I have a raster, with table already built, called "shed". I simply want to export two fields, "Value" and "Count", to a .txt or .csv file (called "try" in this case). So far, I haven't even been able to get a script that even looks at the raster's table. So far I have:

outfile=open("E:/Weathering_GIS/try.txt",'w')
rows=arcpy.SearchCursor('shed',"","","Value;Count","")
for row in rows:
    val=row.getValue('Value')
    count=row.getValue('Count')
    outfile.write(val)

But "val" and "count" always fail to become anything- they just remain empty values! Am I missing something? Maybe some sort of import I have overlooked?


texts: Why is pyproj shifted by 33 feet from QGIS?

texts:
I have written a script to take WGS84 coordinates and convert them into a projected system for importing into a CAD program. I thought I had everything working when I decided to use QGIS to confirm everything was working. But the results are about 33 feet off (mostly shifted in the x).

I'm converting everything to EPSG 26741

...
project_to = pyproj.Proj(init=epsg_code, preserve_units=True)
wgs84 = pyproj.Proj(init='epsg:4326')
...
easting, northing = pyproj.transform(wgs84, project_to, lon, lat)

Green is QGIS projection and yellow is Pyproj projection.