texts: How to create folders when exporting OL to KML

texts:
I have this function that export all plotted features on my OpenLayers (2.6) map, by using the OpenLayers.Format.KML as follows:

var featuresToExport = [];
for (var i = 0; i < MAP.__layers.length; i++) {
    if (MAP.__layers[i].getVisibility() && MAP.__layers[i].features != null) {
        for(var j = 0; j < MAP.__layers[i].features.length; j++){
            featuresToExport.push(MAP.__layers[i].features[j]);
        }
    }
}

var format = new OpenLayers.Format.KML({
    'maxDepth':10,
    'extractStyles': true,
    'extractAttributes': true,
    'extractTracks': false,
    'internalProjection': MAP.__map.baseLayer.projection,
    'externalProjection': new OpenLayers.Projection("EPSG:4326")
});

var kml = format.write(featuresToExport);

That way, all features are put at the same level on the KML. I would like to be able to arrange them in folders, for example for each layer.

How can I do that? The write function always write the complete KML structure, otherwise I could do something like this:

for (var i = 0; i < MAP.__layers.length; i++) {
    if (MAP.__layers[i].getVisibility() && MAP.__layers[i].features != null) {
        featuresToExport = [];

        kml += "<Folder><name>Layer " + i + "</name>";

        for(var j = 0; j < MAP.__layers[i].features.length; j++){
            featuresToExport.push(MAP.__layers[i].features[j]);
        }
        // write this features
        kml += format.write(featuresToExport);
        kml += "</Folder>";
    }
}


texts: Categorise vector layer based multiple attribute variables

texts:
I would like to make a classification of a dataset (WFS point layer) on the values of two different variables. I would like to select from the first column only a specific type of flower (Rose). Then, based on the second column, I need a condition which returns eighter yes- or no based on the values in this column.

If the value contains the '<' symbol it should classify as 'no', if it contains a value without the '<' it should return yes.


texts:
I am mid level GIS Analyst\GPS Coordinator who would like to become a GIS mobile developer at some point. At my workplace we work in the ESRI Platform,ArcGIS Server etc
Being a GIS Analyst without any programming experience under my belt.I have started taking java programming classes.At some point in time i would like to create native apps.

However,right now I am not sure how to get there. Should I start off by creating apps using ArcGIS server and then branch off in to native apps etc?

So I guess my question is where do I start? Formal training, languages etc.There is so much information but i am having a hard time figuring out where to go.

Hopefully some of you can shed some light from the prospects of someone who is experienced  professional in the field.
Thank you for your time and my apologies for the vagueness of my question.


texts: Is ANZLIC Metadata supported using the new approach of ArcGIS 10.x?

texts:
I like the new approach to metadata at ArcGIS 10 because I fall into the camp that prefers minimal metadata for many purposes, but like to have the option to work with more comprehensive standards when they are appropriate.
My understanding is that the ANZLIC Metadata Profile 1.1 (based on ISO 19115) is very similar to the ISO 19139 Metadata Implementation Specification so can the ArcGIS 10 metadata style called ISO 19139 Metadata Implementation Specification simply be used or is the ArcGIS/ANZLIC community awaiting the release of an ANZLIC specific metadata style for ArcGIS 10 by Esri or a third party?
A reply to an ArcGIS Discussion Forum thread suggests that it might be.


texts: Randomly replace a set number of cells with a value from a unique list

texts: I have an Excel/ASCII table (300x300) where I want to randomly select a certain number/percentage of cells and replace that value with a number from a list (0.131, 0.149, 0.052). Is there a way to do this? I am going to convert the ASCII file to a Raster. Or is there a way to do this in ArcMap?


texts: Installerv otb 6.6.1 under QGIS 2.18.24 in Windows 7

texts:
I had installed QGIS 2.18.24 on OSGeo4W, then the OTB stand-alone binary file. In QGIS, I manually entered the directory of OTB Applications folder

C:\Program Files\QGIS 2.18\apps\OTB-6.6.1-win32\lib\otb\applications)

and command line folder

C:\Program Files\QGIS 2.18.24\apps\OTB-6.6.1-win32\bin

and I checked the Enable button.

But nothing happens when I close the Processing window and no OTB menu appears in the Processing toolbox.

How can I solve this problem?


texts: Correcting a table's projection in MapInfo

texts:
Occasionally I have a MapInfo TAB file that doesn't match up with others. Presumably this is because MapInfo is interpreting it as a different projection to the one it was saved as. I am also assuming that it is possible to have two TAB files of different projections in the same workspace and have them match up.

As far as I can see the only way to change a tables projection is to save a copy which doesn't seem to help me.

I there a way to switch projections on the fly? Or at least to find out which projection the TAB file thinks it is and which projection MapInfo is interpreting it as? Am I just being thick?


texts: Earth Engine, Identify overlapping features

texts:
I have a feature collection of badly drawn geometries (a subset is shared with the script).  I need to identify the features that overlap in some way (for later processing).

The idea was to map each feature of the collection to every other feature in the collection and identify if they are overlapping. And set the number of overlaps in the property of the feature.  If the number is greater than 1, then it overlaps with more than itself.  However, I always end up with 12. (the number of features in my collection).

var fields = ee.FeatureCollection("users/JASPR/Results/fields");

var intersect = function(feature){

  var fc = fields.map(function(f){
    var intersects = feature.geometry().intersects(f.geometry())
    return(f.set({intersects: intersects}))
  })
  var status = fc.aggregate_count("intersects")
  return(feature.set({overlaps: status}))

}

var newFields = fields.map(intersect)

print(newFields)

Map.addLayer(newFields)

What do I need to change to get the correct number of overlaps?


texts: Split lines in segments based on distances in Oracle

texts:
I have a data set consisting of 2D linestrings which are stored in an Oracle database in table my_lines with an id and a geometry column (units are in meters).

In a separate table my_intervals with a line_id, a start and an end column, I have a set of intervals with start and end distances describing parts of the linestrings in my_lines - referenced by id and line_id. The intervals do not overlap.

I would like to split the lines in my_lines using the intervals in my_intervals so I get the line segments representing the intervals and the remaining line segments not represented by intervals.

Example:
In the following example I have a linestring which is 600m long, and two intervals, [0m; 200m] and [400m; 600m]:

So far I have looked at the SDO_LRS.CLIP_GEOM_SEGMENT but I must admit that I am lost on how the spatial parts of Oracle works compared to postgis, FME, OGR and similar, and furthermore how it can be combined into returning the different line segments.

Question:
How can this be done using some more or less magical Oracle SELECT statement?


texts:
I have some code in c++ that used gdal and process raster and create an output raster.
I want to create a button in ArcGIS Desktop (using ArcObjects SDK 10.2 c#) and pass the raster layer to my c++ code (library), break the gdal dependency, and send back the our raster to ArcGIS for display. Is there an example project which I could use for template.


texts: Reprojecting from MGRm to WGS84 using open-source tools

texts:
I would like to work with some of the Michigan state data (for example, this shapefile of cities ), but it is not in a common projection. However, when I try to reproject the file using ogr2ogr and any one of the possible EPSG codes, I get a map that is anywhere from rotated & skewed to flat:

ogr2ogr -s_srs EPSG:102121 -t_srs EPSG:4326 city_mi_84.shp city_mi.shp

Any suggestions on how I can make this work? I don't have access to Arc at the moment.


texts: Pull pixmap image source from ui edit text in QT

texts:
I'm putting together a custom form for QGIS identify using Qt Designer.
One attribute has a file name to the accompanying photo of the feature. I'd like to have the photo label in the form grab the file source data from the edit text box that shows the file name on load. Then if the file name is changed or added, the picture can also update.
More specifically, the edit text may show "001", though the file name may be "C:\Photos\IMG_001.jpg".
Can pixmap be set up in Qt Designer to do this?


texts: Labelling column plot with unit on top in QGIS?

texts:
I made a column representation in QGIS, which shows the effects on different points around point 24.

Is it possible to label the columns on top with the unit as you are used to from excel?

QGIS-Image:

Excel-Image:


texts:
I  tend to work with personal GDB and with complex labels, when done I convert them to annotations, creating the annotation feature class and the "Anno" relationship, both into the personal GDB.

How can I do if I want to create a copy of that set of annotations to change a few of them?

In the end I need to keep both, the original one and the amended one.

If I try to copy the annotation directly in ArcCatalog not only the Annotation feature class is copied, the "Anno" relationship and the original shapefile (where the labels originally come from) is copied as "Shapefile_1"


texts: Creating svg item on Map Canvas in QGIS API?

texts:
I'm working on a project with QGIS-2.4.0 API using C++. I need to load some svg files to the application for displaying as symbols. QGIS API has a class named QgsMapCanvasItem, which is an abstract class for items that can be placed on the map canvas. However, there are no subclasses related to svg symbol. So I created an object using QGraphicsSvgItem from Qt, here is the code:

int main(int argc, char ** argv)
{
QgsApplication app(argc, argv, true);

QgsMapCanvas * mypMapCanvas = new QgsMapCanvas(0, 0);
mypMapCanvas->enableAntiAliasing(true);
mypMapCanvas->setCanvasColor(QColor(255, 255, 255));
mypMapCanvas->freeze(false);
mypMapCanvas->setVisible(true);
mypMapCanvas->refresh();

QGraphicsScene *s = mypMapCanvas->scene();
QString fileName = "D:\\qgis\\svg\\symbol\\poi_tower_communications.svg";
QGraphicsSvgItem *m_svgItem = new QGraphicsSvgItem(fileName);

m_svgItem->setScale(0.1);
s->addItem(m_svgItem);

return app.exec();
}

But something is wrong, when I click the Minimize button or close button, a debug error occurs, the picture is the debug error.

I know the concepts of Qt Graphics view, scene, item. QgsMapCanvas is derived from QGraphicsView, QgsMapCanvasItem is derived from QGraphicsItem, but I couldn't find out which class is derived from QGraphicsScene?


texts: Stack - Mapnik, TileMill, OpenLayers, GeoExt

texts:
Would this work as a stack for serving web maps with the basic GIS viewer functionality (pan, zoom, layer switching, info popups)?

Apache Web Server
Mapnik to serve the geo-data (TileMill to Style the Mapnik data)
OpenLayers to display in browser
GeoExt or JQuery for added functionality

Do I have that right? Am I missing anything?
I've read through a bunch of the threads and other web map info, still not sure exactly where/how Mapnik fits in. It can replace GeoServer/MapServer in a stack?

And TileMill, not sure how the styled map gets added to the server.

Thoughts and advice very much appreciated!


texts:
I have some configuration or something like that where i can do to limit the size of log geoserver (like 1 Gb), or maybe delete entire log like 1 time of week?

Because my logs never stop to increase, like 5 GB, 7GB etc


texts: Point displacement between map (Cesium 3D) and raster (numpy array)

texts:
It is better to show than tell:

On this image, the layer below the colorful shape is a trivial heat map in raster format.
The colorful shape is the product of the skimage.graph.MCP_Geometric method.

In order to find the cell indexes for the start and end points, I am using the calculation shown in the second answer of this thread: How to extract value from a raster based on lat and lon? (I have used the first one also as a test, ending in the same results.).

The problem:

As you can see from the arrows on the resulting image, the start/end points seem to be displaced in accordance to the starting points of the corridor.
The points' coordinates are in srid 4326.

UPDATE:

The Point coordinates are taken by Cesium with 3d tile rendering. When checking the same rasters and the corresponding points in QGIS I see an even greater displacement:


texts: Breaking out single polygon (feature) of SpatialPolygonDataframe

texts:
I am trying to break out a single feature of a Spatialpolygon object.
That is e.g. plot a single polygon from a shapefile that contains multiple.
Any suggestions how to do this?

>#Import shapefile, generate spatialpolygondataframe
>vectormap<-readOGR(dsn="QGIS export", layer="section4")
>
>OGR data source with driver: ESRI Shapefile
>Source: "QGIS export", layer: "section4"
>with 4 features
>It has 6 fields
>
>
># Change vectormap CRS to Equirectangular Carree projection
>vectormap2 <- spTransform(vectormap, CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 >+x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

> #Summary of vectormap2
>Object of class SpatialPolygonsDataFrame
>Coordinates:
>       min      max
>x  1467755  5288689
>y -6018287 -3018986
>Is projected: TRUE
>proj4string :........
>Data attributes:
>  Layer                       SubClasses ExtendedEn Linetype EntityHand Text
> Layer 1:4   AcDbEntity:AcDb2dPolyline:3    NA's:4     NA's:4   41:1    NA's:4
>             AcDbEntity:AcDbSpline    :1                        5B:1
>                                                                7B:1
>                                                                94:1
> coordinates(vectormap2)
>     [,1]     [,2]
>0 4122892 -3791580
>1 4661889 -4500891
>2 4500413 -4079798
>3 3352295 -4392646

Happy for any input!

Side Note.
I can also get the shapefile in a format giving the features in separate layers.

> summary(vectormap)
>Object of class SpatialPolygonsDataFrame
>Coordinates:
>           min      max
>x  -0.01383442 457.4256
>y 125.27349310 484.3477
>Is projected: FALSE
>proj4string :............
>Data attributes:
>         Layer                           SubClasses ExtendedEn Linetype >EntityHand   Text
> big        :1   AcDbEntity:AcDbBlockReference:4    NA's:4     NA's:4   42:1       >NA's:4
> fullsection:1                                                          48:1
> medium     :1                                                          67:1
> small      :1
84:1

However when I do this, I cant change the projection.

vectormap2 <- spTransform(vectormap, CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
non finite transformation detected:
[1] 187.4585 468.4831      Inf      Inf
Error in .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  failure in Polygons 1 Polygon 1 points 1
In addition: Warning message:
In .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  :
  1808 projected point(s) not finite


texts: How to add Attachment to an object in ArcGIS Online?

texts:
I created attachments for some objects of a theme (as a File Geodatabase Feature Class) in  ArcGIS 10.0:  images (JPG) and doc files.

Then moved the same shp theme in my web-map of ArcGIS Online, and want to add in any way these attachments for these objects there as well.
So, my main goal -  find a way to see this attached jpg and doc files when I choose concrete object in the web-map of ArcGIS Online.

Also I’d like to find opportunity to place in the web-map of ArcGIS Online  - URL in attribute table for any concrete object.

Could you advise me please – how it’s possible?


texts: Leaflet displaying oddly curved TMS tiles

texts:
I'm trying to display some data using arctic ice extent data (projected in EPSG:3413, NSIDC Polar Stereographic North). I have processed the distributed images using gdal2tiles with the following command:
gdal2tiles -p raster -s EPSG:3413 -z '0-5' -e -n tmp.vrt
(tmp.vrt is the RGB-ified file I created from the TIFF using gdal_translate)

When I try to display the tiles with Leaflet (or with OpenLayers, for that matter), at zoom=0 the page has a line of oddly curved versions of my tiles, flat along the top but looking like they've had big circular bites taken out of them along the bottom. Zooming in causes the proper tiles to mostly appear, but there is still the line of curved ones cutting across it, obscuring a large area of the actual images.

I assume this is some kind of projection issue, and I have tried correcting it using a Proj4Leaflet CRS, but all that does is leave the oddly-shaped tiles and badly rearrange those few tiles it displays sensibly.

I'm fairly new to this, and can't figure out if the issue is in the way I've processed the tiles with GDAL, or in the way I'm trying to display them with Leaflet.

Here is a screenshot of the display at zoom=0:

And at zoom=2

Here's the level-0 image the page is trying to read from. All other tiles are just rectangular slices of this image: 0.png

And the relevant code:

  var map = new L.map('map').setView([90,0], 0);

  L.tileLayer('data/20170112/{z}/{x}/{y}.png', {
    tms: true,
    nowrap: true,
    attribution: "Map data courtesy University of Bremen",
    maxZoom: 5
  }).addTo(map)

In case this is the problem, here's my entire "processing chain", such as it is:

$ wget http://www.iup.uni-bremen.de:8084/amsr2data/asi_daygrid_swath/n3125/2017/jan/Arctic3125/asi-AMSR2-n3125-20170112-v5.tif
$ gdal_translate -of vrt -expand rgb asi-AMSR2-n3125-20170112-v5.tif tmp.vrt
$ gdal2tiles -p raster -s EPSG:3413 -z '0-5' -e -n tmp.vrt 20170112


texts: How to use rsgislib python binds in windows environment

texts:
I am right now would like to use rsgislib to create a simplify remote sensing image processing application based on python.
However I have no idea about how to use rsgislib in windows.

Perhaps someone would like to share about using rsgislib in windows or a substitute python library for rsgislib.

rsgislib: The Remote Sensing and GIS software library (RSGISLib) is a collection of tools for processing remote sensing and GIS datasets. The tools are accessed using Python bindings or an XML interface.
http://www.rsgislib.org/


texts: Using spatial regression on interpolated surfaces?

texts:
I have two raster surfaces where I want to test to see If any relationship (A depending on B) can be found and if so if that relationship is statisticaly significant.

What are key issues to consider and look out for?

How would you propose that I go forth?

Both surface A and B are interpolated surfaces which are indexed (internally, that is to say they do not share the same 'unit' or range, their distributions/histograms are however similar but not identical). The question I want to answer is if their geographical distribution correlates.

I'm especially concerned with the risk of having induced spatial autocorrelation. I already now that one of the raster datasets is spatially autocorrelated and suspect that the other would yield similar results since they both show quite clear geographical patterns.

Hence what should I be wary of when doing a regression analysis on two interpolated surfaces?

I'm also seeking good references and further reading in the subject.

I don't consider myself an expert at spatial statistics.

I got access to ArcMap with all relevant licenses/extensions. I'm not familiar with R nor with advanced scripting.


texts: Automatically panning the map in ArcGIS when using the Time slider

texts:
I have a point dataset with time-stamps at 1 minute intervals showing a driving route.

I have enabled the Time properties of the layer and I am able to use the Time slider to see the points move along the route.

What I would like to do is to zoom in closer to the map so that I can't see the full extent of the route, and when I press the "play" button on the time slider, for the map to automatically pan with the point as it moves in time.

Is this possible using just the desktop functionality?

I am doing this on ArcInfo 10.1, but I could do it in ArcInfo 10.2.1 if there are new features on there.

I do also have a gps enabled time-stamped video of the same driving route if this makes the process any easier?  But I guess linking the video to the route in Arc is somethign for later!


texts:
I have a set of polygon layers with thousands of polygons some of which are adjacent (Either rings or simply adjacent). I would like to simplify this polygons by merging them.

Is there a way to do this automatically?


texts: There are multiple rows in my .tab file in mapinfo. Is there any way to select nth row in mapbasic window of mapinfo?


texts: How to reduce 2 image bands simultaneously using grouped sum region reducer in GEE

texts:
I am trying to reduce a region to simultaneously obtain forest loss and forest cover by year from the Hansen dataset. I have tried doing this using a grouped reducer, however, I get the following error:   "Dictionary (Error)
ReduceRegion.AggregationContainer: Need 2 bands for Reducer.group(Reducer.sum), <Image<[forest_area, forest_loss, lossyear]>> has 3."
Is there an inherent reason why this is not possible or am I doing something wrong? From one of the GEE tutorials, it seemed like it should be possible to do this but I have found no examples.
This is the immediate code I am using with the first 3 variables being modified bands from the original dataset.
var protectedLoss = dataset.select(['loss']).clip(protectedAreas);
var protectedLossYear = dataset.select(['lossyear']).clip(protectedAreas);
var protectedForestLoss = blank.where(treeCover.gt(forestDef).and(protectedLoss.eq(1)),1).clip(protectedAreas);
var protectedChart = protectedDataset.reduceRegion({reducer: ee.Reducer.sum().group({
groupField: 2}), geometry: protectedAreas,  scale: 30,  maxPixels: 1e12`

Apologies if this is a basic question - I am new to GEE.


texts: Store elements of Imap into BLOB fields using IPersistStream

texts:
I have an IMap (include elements: text, line) and I want to store it into a field of a table in SQL, so I can load it again later. I used below codes, but I got an the error :

public class Test
{
  IStream stream;
  IMap myMap;
  void Save()
  {
    IPersistStream perStream;
    IMemoryBlobStream memoryStream = new MemoryBlobStreamClass();
    stream = memoryStream as IStream;
    perStream = (IPersistStream ) myMap ;
    perStream.Save(stream, 0);
  }
  void Load()
  {
     IMemoryBlobStream memoryStream = new MemoryBlobStreamClass();
     IPersistStream pstm = new MapClass();
     pstm.Load(stream); **// throw error: Error HRESULT E_FAIL has been returned from a call to a COM component.**

     myMap = (IMap)pstm;
  }
}

I found out one way to save elements to an binary file and can load it again later as below codes:

void Save (AxPageLayoutControl pageControl)
{
                IPersistStream perStream;
                IMemoryBlobStream memoryStream = new MemoryBlobStreamClass();

                IGraphicsContainer graphicsContainer = pageControl.ActiveView.GraphicsContainer;
                graphicsContainer.Reset();
                IElement pElement = graphicsContainer.Next();
                while (pElement != null)
                {
                    if (pElement is IMapFrame)
                    {
                        perStream = (IPersistStream)pElement ;
                        perStream.Save(memoryStream, 0);

                        memoryStream.SaveToFile("C:/image.bmp");
                        break;

                    }
                    pElement = graphicsContainer.Next();
                }
}

void Load (AxPageLayoutControl pageControl)
{
            IPersistStream perStream;
            IMemoryBlobStream memoryStream = new MemoryBlobStreamClass();
            IGraphicsContainer graphicsContainer = pageControl.ActiveView.GraphicsContainer;
            graphicsContainer.Reset();
            IElement pElement = graphicsContainer.Next();
            while (pElement != null)
            {
                if (pElement is IMapFrame)
                {

                    memoryStream.LoadFromFile("C:/image.bmp");
                    perStream = (IPersistStream)pElement;
                    perStream.Load(memoryStream);
                    IMapFrame frame = new MapFrameClass();
                    frame   = (IMapFrame)perStream;
                    pElement = (IElement)frame;
                    break;

                }
                pElement = graphicsContainer.Next();
            }
            pageControl.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);
}

But, I am finding a way to convert IPersistStream or IMemoryBlobStream to byte[] format, so I can store in SQL Server and load later.

Does anyone know how to do this?


texts: Using advanced settings of GDAL_rasterize in QGIS script

texts:
I’m working on a script (that later should be used in a processing model) that uses gdal_rasterize to write different attribute values (columns) in different bands of one raster file.

I know that this I possible with gdal_rasterize using the -b and -burn flags.
Anyway, I can’t manage to integrate this commands in my QGIS-Python script. Am I missing something with the syntax?

Here is the GDAL synopsis:

Usage: gdal_rasterize [-b band]* [-i] [-at]
   {[-burn value]* | [-a attribute_name] | [-3d]} [-add]
   [-l layername]* [-where expression] [-sql select_statement]
   [-dialect dialect] [-of format] [-a_srs srs_def]
   [-co "NAME=VALUE"]* [-a_nodata value] [-init value]*
   [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]
   [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/
         CInt16/CInt32/CFloat32/CFloat64}] [-q]
   <src_datasource> <dst_filename>

And here is my script so far:

 ##Rasterize_GDAL_KRK=name
 ##Input=vector
 ##Output=output raster
 ##Extent=extent
 ##cellsize=number 25

import os, processing, osgeo.gdal

processing.runalg('gdalogr:rasterize',
                                        {"INPUT":Input,
                                        "FIELD":"KRK_1",
                                        "DIMENSIONS":1,
                                        "WIDTH":cellsize,
                                        "HEIGHT":cellsize,
                                        "RAST_EXT":Extent,
                                        "TFW":0,
                                        "RTYPE":0,
                                        "NO_DATA":None,
                                        "COMPRESS":0,
                                        "JPEGCOMPRESSION":1,
                                        "ZLEVEL":1,
                                        "PREDICTOR":1,
                                        "TILED":False,
                                        "BIGTIFF":3,
                                        "EXTRA": None,
                                        "OUTPUT":Output})


texts:
Can I set the title of an SLD rule using a feature's field?

I have a very dynamic dataset where the value of a polygon changes after each update, therefore I need to be able to compare using one field and label the style rule with another field to create a legible legend.

Below I have a rule using the SLD syntax. I am able to use a PropertyName tag to compare two values, is it possible to similarly use a PropertyName tag (or equivalent) to change the Description's Title tag?

    <se:Rule>
      <se:Name>0</se:Name>
      <se:Description>
        <se:Title>0</se:Title>
      </se:Description>
      <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
        <ogc:PropertyIsEqualTo>
          <ogc:PropertyName>myFieldName</ogc:PropertyName>
          <ogc:Literal>1</ogc:Literal>
        </ogc:PropertyIsEqualTo>
      </ogc:Filter>
      <se:PolygonSymbolizer>
        <se:Fill>
          <se:SvgParameter name="fill">#820000</se:SvgParameter>
        </se:Fill>
        <se:Stroke>
          <se:SvgParameter name="stroke">#000001</se:SvgParameter>
          <se:SvgParameter name="stroke-width">1</se:SvgParameter>
          <se:SvgParameter name="stroke-linejoin">bevel</se:SvgParameter>
        </se:Stroke>
      </se:PolygonSymbolizer>
    </se:Rule>

In GeoServer simply adding a PropertyName tag returns the following:

  Element 'se:Title' is a simple type, so it must have no element
  information item [children].

I've tried adding a function within the tag and get the same result. Adding a PropertyName tag within the Description tag returns:

  line 11: cvc-complex-type.2.4.a: Invalid content was found starting
  with element 'ogc:PropertyName'. One of
  '{"http://www.opengis.net/se":Title,
  "http://www.opengis.net/se":Abstract}' is expected.


texts: Routing by means of OpenStreetMap with Qt and Maps.me?

texts:
I'm following the steps on this page:

https://github.com/mapsme/omim/blob/master/docs/MAPS.md
https://github.com/mapsme/omim/blob/master/docs/CONTRIBUTING.md

When I run this line "omim/tools/unix/generate_mwm.sh source.pbf any_random_string_asdf" I don't have an Error, then run Qt and Maps.me application appears but when I click on the route nothing happens.

What do I have to do? Should I execute something before?


texts: Seeking process in ModelBuilder that can choose among several shapefiles in folder?

texts:
I created a model in ArcGIS 10 which selects those parcels that intersect my polygon shapefile.

It will then remove all r.o.w parcels.

I then export the table it produces to an excel file. The shapefiles are in a folder.

Is there a process which will direct ModelBuilder to said folder then prompt me as to which shapefile I want to run without using python?


texts: Add Bing basemap to ArcMap with ArcGIS Portal

texts: I am aware that with a Bing maps key, one can add a Bing basemap to ArcMap by adding data from ArcGIS Online and choosing a Bing basemap. My issue however is that we use Portal and not ArcGIS Online. We have a Bing maps key and I am looking for a way to add a Bing basemap to ArcMap.


texts:
I'm trying to read in several rasters as array using gdal in a python script, but I keep getting NoneType returned and can't figure out why.

The process is this:

Use gdalbuildvrt (http://www.gdal.org/gdalbuildvrt.html) to read in several rasters into gdal's virtual format.
Loop through the rasters (aka the layers) of the .vrt, reading them into a temporary array and adding the temp array to the cumulative (or total) array.
Write the total array to the raster.

import os
import gdal
from osgeo.gdal import *
import numpy as np
gdal.UseExceptions()

# intxtlist = a text file with the rasters to be put in the vrt; in this case there are 3 files
VRTstack = '/path/tempVRT.vrt'
os.system('gdalbuildvrt -separate -vrtnodata 0 -input_file_list %s %s' % (intxtlist, VRTstack))

vrtDS = gdal.Open(VRTstack)
Nlayers = vrtDS.RasterCount # Nlayers = 3 because there were 3 input rasters
Nsamp = vrtDS.RasterXSize
Nline = vrtDS.RasterYSize

totalArr = np.zeros((Nline, Nsamp), dtype = 'uint16') # preallocate an array for the final output total
for N in range(0, Nlayers):
    tempArr = vrtDS.GetRasterBand(N+1).ReadAsArray() # tempArray is now Nth tiff of stack as array; N+1 b/c gdal starts at 1
    print type(tempArr)
    totalArr += tempArr  # add the temp array to the final array
## code here to write totalArr as raster

The problem I'm getting is with the second N in the for loop. When N=0, "print type(tempArr)" returns: type 'numpy.ndarray', which I would expect. When N=1, however, type(tempArr) returns type 'NoneType', which I can't understand. The entire script then fails at line "totalArr += tempArr" because you can't add NoneType to a numpy int array.

What's weird is that no matter how I loop through, the "second" N always fails. For example, "for N in [1, 0, 2]:" returns an array for N=1 and N=2 but returns NoneType when N=0 (even though N=0 didn't fail when going from 0 to Nlayers (aka 0, 1, 2)) and "for N in [2, 1, 0]:" works for N=2 and N=0 but not N=1.

I'm pretty sure there's nothing wrong with the vrt itself, as 'gdalinfo %s' % VRTstack' returns accurate information on all three bands.

Also worth mentioning is that this script worked with one set of images (call it S1), but not another (S2). The only real difference between each set is that S1 has a coarser resolution, and its rasters are 32-bit signed type, while S2's raster type is 8-bit unsigned. I thought that maybe it had something to do with the fact that I was trying to read in unit8 rasters and add it to a uint16, so I tried making totalArray type uint 8 as well, but it still failed. It later occurred to me that this wouldn't change anything, because tempArr is NoneType before it even gets added to totalArray.

When using 5 images instead of just 3, NoneType is always/only returned for the 2nd and 4th N in the loop. That is, when looping through N=0 to N=4, numpy arrays are returned for N=0,2,4, while NoneType is returned for N=1 and N=3. Similarily, when the loop is "for N in [1, 0, 3, 2, 4]," arrays are correctly returned for N=1,3,4 but NoneType is returned for N=0 and N=2. Bizarre.


texts: Creating and populating inset maps for densely-labelled areas using ArcPy?

texts:
I would like to know if it is possible to automate the example in the picture below, unlike the picture below i'm using vector data.

I would like for a fixed scale generate map, but when labels are too dense, we could have some little maps (inset maps) on the side of the Main map.

i have found some tips using arcpy.

but the problem is how to detect densely labelled areas


texts: Sentinel-2 Hub & satellite imagery data format

texts:
I acquire the raw JP2 imagery files from the Sentinel-2 project. Let's say I've selected a tile for a given sampling period and band (e.g. B01). Then, after using Python tools provided from here, I compare my JP2 files to the TIFF files (corresponding to the exact same spatial area corresponding to the JP2 file) provided via get calls to the Sentinel Hub AWS S3 source. These two files - regardless of sampling period, location, or (usually) band are different in the following ways: 1] array dimensions (e.g. JP2=1830x1830 vs. TIFF=10980x10980) and 2] measurement values (e.g. JP2=[[1000,1001,1023,...]] vs. TIFF=[[0.780,0.888,0.901]]).

This means two things. First, some kind of interpolation is going on from the Sentinel-2 raw data format to resample a coarser array to one comparable to more higher resolutions bands like B02, B03, etc. Also, the array value themselves corresponding to a given TIFF image have also undergone some form of normalization/standardization/regularization/equalization.

What does Sentinel Hub do - or perhaps the Python Sentinel-2 toolkit - to the images to make this large difference happen compared to the supposed raw JP2 format? That is, why the radical difference between the JP2 and TIFF files?


texts: How to add a node/vertex in a QGIS edit session?

texts:
I'm doing some polygon (boundary) cleanup, and am wondering if it is possible to add a node/vertex to a set of selected polygons or even a single polygon.

By the way, the QGIS topology editing tools are incredible - extremely simple but very accurate and take so much of the work away from doing such manual work.

QGIS gives the ability to remove selected nodes - amazing that the topology is preserved and even improved when cleaning up excess nodes. Great stuff!


texts: Implement security to qGIS Server (for windows)

texts:
I have few map services running on QGIS Server(v 2.2.X on Windows 7 32 bit) and anyone can connect or access those services (WMS & WFS).

Is there any way to implement security in QGIS Server? So that I can set, assign username & password to those map services.

I researched on my side and following few links in which provides few steps such as apply password to specific folder, creating .htaccess file, also researched official site of QGIS Server but cannot find any valid document.

Any help will be great!!! Thanks in advance !!!


texts:
We have a stack of open source software's which we use to build a GIS.
We would like to align our application to OGC as much as possible. I visited the site and read about the standards.

Now my question is how to start with making the application OGC compliant. There are lot of documents which talk about standards but I couldn't find much on how to get started. Specifically given that am building on top of Open Source software's and Implementing is very subjective word.

Can you guide me to some resources which can get me started?

Ideally we would like address Web Mobile and Desktop, We have a existing solution for Desktop but are planning to go ahead with Web. We are aiming for a GIS product.

Here is the list of Open Source Products: GDAL/OGR, GeoServer, OpenLayers, Cesium


texts: Leaflet L.Print (Mapfish) cannot set map scale

texts:
I am trying to print a map using the MapFish Print extension of Geoserver in Leaflet (L.print). I configured the .yaml file to suit my needs but the scale of the map doesn't change no matter my settings.

On the printed map it always displays at 1:2000000.

Any idea how to enforce the map scale on the print tool?

//Create a PRINT provider
    var printProvider = L.print.provider({
    capabilities: printConfig,
    method: 'GET',
    autoLoad: true,
    dpi: 300,
    legends: true,
    outputFormat: 'pdf',
    //layout:"A4 Landscape",
    pages: [{center: [46, 45.5],
            scale: 10000,
            dpi: 300,
            geodetic: true}],
    customParams: {
              mapTitle: 'title here',
              comment: 'here we can include a comment regarding the map - standard text'
          } });

//Create Print Control
var printControl = L.control.print({
provider: printProvider
});
map.addControl(printControl);`


texts: How to map OSM XML to screen coordinates?

texts:
I'm creating a simple application that allows visualizing OSM XML data.

I've been following JOSM as a guide, and I'm trying to understand how you can map the long, lat of the XML data to x, y screen coordinates.

I see that there's a link to slippy map tiles that gives an equation, to map to a tileset once given a zoom factor, long and lat. This may be related to my question but I don't think it gives the proper answer.

I'm good with graphics programming and I understand how to use the bounding box to get the screen coordinates if everything was given in an x,y coordinate information. But, everything is in the lon and lat, which means there has to be a projection.

How does JOSM or another application go from lon, lat to x, y on a flat screen simply using the OSM xml data(like JOSM)?


texts: What should I start with Web GIS MapServer?

texts:
I'm not native English Speaker and I'm a graduated student who wants to have a masters degree in software engineering.  I have a thesis relative with "Build a Web GIS for air pollution in Viet Nam's capital city". But I've never any experiences with GIS. I've searched some Google, this shows me that I should use MapServer, OpenLayers and PostGIS SQL.

Now, I think user will query from web browser and based in OpenLayers function JavaScript, It'll transfer data to MapServer (I think can use PHP MapScript here to handle query), after that MapServer can get data from PostGIS SQL if it needs (after all, MapServer will return a image that fits all user's query).

Please, correct me if I'm wrong or teach me some knowledge.


texts: ArcMap closing when running ArcPy script?

texts:
In the below script the processed file will be added to the table of content in ArcGIS, since my data is too long it cause the program to be closed due to opening too many files, how I can run the script without opening the processed file in ArcGIS.

import  arcpy
ws = arcpy.env.workspace = r'D:\BRB Snow Cover\Sample\Terra'
Terra_Rasters = arcpy.ListRasters("MOD*")
matches = 1
for item1 in Terra_Rasters:
    for item2 in Terra_Rasters:
        item1_num = int(item1[31:38])
        item2_num = int(item2[31:38])
        if (item2_num - 1) == item1_num:
            print item1
            print item2
            print "Match"
            output = "D:\\BRB Snow Cover\\Sample\\output1""\\" + str(item2) + '.tif'
            arcpy.gp.Con_sa(item2, item1, output, item2, "\"Value\" =8")
        else:
            print "No"
            matches = matches + 1


texts:
I am dealing with some large datasets and would like to minimise the amount of data sent over the network. I am using OpenLayers' BBox strategy to only transfer the data that is needed but I am wondering whether I can simplify the output so that I can show lower zoom levels without having to transfer data with unnecessary accuracy/detail.

I am using GeoServer as a WFS server as part of the OpenGeo bundle.


texts: Access-Control-Allow-Origin error in web map

texts:
I am having trouble getting a WFS layer working using Geoserver, OpenLayers and GeoExt.

I get the following error:

XMLHttpRequest cannot load http://10.0.0.155/geoserver/wfs. Origin http://10.0.0.155:8080 is not allowed by Access-Control-Allow-Origin.

Does anyone know how I get around this? I can't find anything in the examples for geoserver, openlayers or geoext.

I don't get any errors when using WMS on my maps.

Many thanks,

James


texts:
I'm basically cribbing code from Esri and other online sources to convert a CSV to a NetCDF. I get an HDF error on writing values to the variables. X and Y values are UTM coordinates if that matters. How do I get this to run correctly? I've tried using both x_ran and x_val as the write for x_var, and both throw the HDF error. System is Windows, Python37. Code below:
infile= 'Path\to\in\file.csv'
indata = np.loadtxt(infile,delimiter=',',skiprows=1)

outfile = 'Path\to\out\file.nc'
outdata = Dataset(outfile, 'w', format='NETCDF4')
matrix = indata
matrix = matrix[matrix[:,1].argsort()]
matrix = matrix[matrix[:,2].argsort(kind = 'mergesort')]
matrix = matrix[matrix[:,3].argsort(kind = 'mergesort')]

x_val = matrix[:,0]
y_val = matrix[:,1]
z_val = matrix[:,2]
D_val = matrix[:,3]

x_ran = len(x_val)
y_ran = len(y_val)
z_uni = list(set(z_val))
z_uni.sort()
z_ran =  len(z_uni)

outdata.createDimension('X',x_ran)
outdata.createDimension('Y',y_ran)
outdata.createDimension('Z',z_ran)

xvar = outdata.createVariable("X", "d", ("X",))
xvar.units = "meters"
xvar.standard_name = "Projected X"

yvar = outdata.createVariable("Y", "d", ("Y",))
yvar.units = "meters"
yvar.standard_name = "Projected Y",

zvar = outdata.createVariable("Z", "d", ("Z",))
zvar.units = "meter"
zvar.standard_name = "Elev"
zvar.positive = "up"

D = outdata.createVariable("D", "d", ("Z","Y","X",))
D.units = "Num"
D.standard_name = "D"

xvar[:] = x_ran
yvar[:] = y_ran
zvar[:] = z_ran
D[:,:,:] = data_matrix
outdata.close()


texts: Merge and Delete with priority overlapping polygons QGIS

texts:
I have two shapefiles with polygons: one is regions for a country and the other is the protected areas inside that country.
I have merged these two polygons with Union, but the regions and protected areas are currently overlapping. I would like the area where they overlap to be only the protected areas (when merging, I would like to assign a higher ranking to protected areas). However, I don't want the entire region underlying protected area to be deleted, JUST where they overlap, so one new map is created with protected areas and regions without gaps. they do not have any attributes in common so I'm not sure how to do this in QGIS.
What it looks like (protected area in grey overlapping green regions):

What I want (one layer with merged protected areas no overlap with regions):


texts: Show popup only with a specific zoom level

texts:
On my map I've a point vector that is showed as heatmap down to one a specific zoom level, after this I show the points normally.

var geoBlogHeatmap = new ol.layer.Heatmap({
  title: 'Heatmap',
  source: sourceGeoBlog,
  opacity: 0.8,
  renderMode: 'vector',
  gradient: ['#253494', '#2c7fb8', '#41b6c4', '#a1dab4', '#ffffcc'],
  blur: 20,
  radius: 20,
  shadow: 250,
  minZoom: 8,
  maxZoom: 13,
});

var geoBlog = new ol.layer.Vector({
  source: sourceGeoBlog,
  style: styleFunction,
  minZoom: 12,
  maxZoom: 18,
});

I've used the official Popup example for show a popup after a click on a feature.

/**
 * Elements that make up the popup.
 */
var container = document.getElementById('popup');
var content = document.getElementById('popup-content');
var closer = document.getElementById('popup-closer');

/**
 * Create an overlay to anchor the popup to the map.
 */
var overlay = new ol.Overlay({
  element: container,
  autoPan: true,
  autoPanAnimation: {
    duration: 250
  },
});

/**
 * Add a click handler to hide the popup.
 * @return {boolean} Don't follow the href.
 */
 closer.onclick = function() {
  overlay.setPosition(undefined);
  closer.blur();
  return false;
};

var feature_onClick;
map.on('click', function(evt) {

    feature_onClick = map.forEachFeatureAtPixel(evt.pixel, function(feature, layer) {
        // console.log(feature);
        return feature;
      });

  if (feature_onClick) {
    // console.log(feature_onClick.getProperties().pk);
    overlay.setPosition(evt.coordinate);
     content.innerHTML =
     '
'+ feature_onClick.getProperties().title + '
'
     + '
'
     + '
' + feature_onClick.getProperties().description + '
'
     + 'Approfondisci'
     + '
'
     + ' Pubblicato: ' + feature_onClick.getProperties().time_since_publication + ''
     }
});

I would like to see the popup only when I can see the points, then from minZoom: 12. How I can do this?


texts: Extract SpPolygon or SpLine vertices in RSpatial

texts:
I aim to calculate the concave Hull from all the points of my database.

I'm trying to extract the vertices:
-of a SpPolygonsDataFrame
-of a SpLinesDataFrame
both are imported by readOGR()

I tried different methods but I'm only extracting centroids.

require(rgdal)
poly<- readOGR(".", "sample")
coordinates(as(as(poly, "SpatialLines"), "SpatialPoints"))


texts: How can I create contours of 0.1 or 0.2 instead of 1?

texts: The default value is 1 and I have salinity values from 33 to 35 and therefore only contours 34 and 35 can be created with the contour interval 1... But I want it smaller so I can see contours 34.1, 34.2, .... 34.9, 35.0.....


texts: Large File Sizes and Print Composer not printing

texts:
I have been working on a heavily detailed map using qgis for the last two years. I am runnig QGIS Essen 2.14.1 on a windows 7 home premium 64 bit with 4gb ram. The file size of the map is around 150mb. i was nearing the the end of the project until this monolithic problem has arisen and i am not able to print the map using the composer. I get a 1cm band at the top of the page which looks great, nothing on the rest of the page and printer loads another sheet and sits there for about 10 minutes or so and then gifts me a blank page.  I initially thought it was a question of upgrading RAM but it seems, from what i have so far found, that there is a bug within qgis when trying to print high res. or  high dpi maps.

Printing at a low DPI is not an option.

The issuses i am having a very similar to a previous question:- How can I fix Print Composer errors due to huge file size?

However the 'workaround' suggested there  doesnt seem to be applicable with Essen. Is there a way through or around this this using just QGIS...


texts:
QGIS is running on my Windows PC.

I have a CSV file on my PC that is remotely updated with values of a raspberry pi's GPS position and some wifi-status readings.

Is it possible to create a vector layer to map out the raspberry pi's position automatically along with the wifi-status readings?


texts:
I have a single-band, georeferenced .tif file that I would like to display on top of a Leaflet.js map, using the ImageOverlay function.  Unfortunately, Leaflet seems only to accept PNG or JPG files, and I haven't been able to produce a PNG file that Leaflet will display.

The original .tif file is in EPSG:102022 and the pixel values are represented by floating-point 64-bit numbers (Float64).  Here is a snippet from gdalinfo:

Origin = (-1253000.000000000000000,-3095000.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Band 1 Block=2433x1 Type=Float64, ColorInterp=Gray
  Min=0.000 Max=28.907
  Minimum=0.000, Maximum=28.907, Mean=0.574, StdDev=1.419
  NoData Value=-1.69999999999999994e+308
  Metadata:
    STATISTICS_MAXIMUM=28.906996389906
    STATISTICS_MEAN=0.57436288189673
    STATISTICS_MINIMUM=0.0004002906752851
    STATISTICS_STDDEV=1.4194311264403

I transformed the original image from EPSG:102022 to EPSG:4326 using gdalwarp, and that seems to have worked fine.  However, when I try to convert the reprojected tif image to a PNG image, I get an error message that the data type has to be changed from Float64 to an unsigned 8 or 16-bit integer.  I tried

gdal_translate -ot UInt16 -of PNG <infile> <outfile>

but for some reason it appears to lose a lot of resolution (when I look at in in QGIS) and the resulting PNG file will not display in Leaflet.  I also tried an online converter which produced a PNG image that looked good in QGIS but wouldn't display in Leaflet.  I can see the bounding box in the correct position on the Leaflet map, but there is just a question mark inside.  My Leaflet code is standard, except for substituting a local file path for the URL:

var southWest = L.latLng(-37.0328788,12.7036709),
  northEast = L.latLng( -26.8167313,36.5804906),
  imageBounds = L.latLngBounds(southWest, northEast);
var imageUrl = '/<local_path_on_my_computer>/image.png';

L.imageOverlay(imageUrl, imageBounds).addTo(map);

I've tried lots of variations.  I'd also like to be able to control the color scale in the png file and I'm not sure how to do that.  I've seen some complicated suggestions using other software to produce tiled layers, but I'm hoping that there's a simpler solution, since this is just a single image.

I'd be grateful for any strategy suggestions.


texts: Very Large Coordinates when Reading Shapefile using Python

texts:
When using Python anad fiona package to read a shapefile of Boston, I am getting very large values in the tens of thousands. Why is this?

# Import Boston shapefile
shapefilename = 'ZIPCODES_NT_POLY'
shp = fiona.open(shapefilename + '.shp')
coords = shp.bounds
shp.close()

Then running coords gives the output

(33869.92130000144, 777617.2998000011, 330800.31099999696, 959741.1853)


texts:
I'm working with jupyter and I have been trying to open TIFF images that I have download from GEE. I have trying to open the images in Jupyter environment.

This is my code:

from PIL import Image
import numpy

img = Image.open("NDVI.tif")

but I get the following error message:

  OSError                                   Traceback (most recent call
  last)  in
  ----> 1 img = Image.open("NDVI.tif")

  ~\Anaconda3\lib\site-packages\PIL\Image.py in open(fp, mode)    2820
  for message in accept_warnings:    2821         warnings.warn(message)
  -> 2822     raise IOError("cannot identify image file %r" % (filename if filename else fp))    2823     2824

  OSError: cannot identify image file 'NDVI.tif'

I have tried to change the raster nae to end with TIFF and tiff but didn't help.
How can I solve this?

EDIT: This is where my jupyter notebook is and where the tif are:


texts: Mapbox.js: Bind dynamic popups with featureLayer?

texts:
I am building a multiple chloropleth map using v2.2.1 of Mapbox.js with featureLayer. Now  I would like to add tooltips with custom content on hover event.

This code is working well to set up the featurelayer and add styles to each layer:

_this.orgLayer = L.mapbox.featureLayer(null, {
            pointToLayer: function (feature, latlng) {
                return L.circleMarker(latlng, {
                    radius: 8,
                    color: "#000",
                    weight: 1,
                    opacity: 1,
                    fillOpacity: 0.97
                });
            }
        })
        .loadURL(boundsUrl)
        .addTo(map)
        .on('ready', joinDataAndSetupMap);
...
// once data joined and map set up, add styles to each layer.
// this gets called at initialization and again when the map is updated.
_this.orgLayer.eachLayer(function(layer) {
        var val = layer.feature.properties.data[month];
        var style = _this.getStyle(val);
        layer.setStyle(style);
    });

But I'm not sure how to bind popups on hover. Could anyone suggest the best way to do this?

The documentation mentions a popupOptions option, but I can't see how to use those to specify dynamic content, or bind the popups to mouseover not click.

Maybe I need to do something inside the eachLayer loop? If so, what?

UPDATE: I've got this, it works but I don't think it's the most efficient way, because it has to re-add the events each time the layer is regenerated. It would be better to add them once at the beginning, wouldn't it:

    _this.orgLayer.eachLayer(function(layer) {
        var style = _this.getStyle(monthData[ratio],
                                   quintiles,
                                   layer.feature.properties.name,
                                   _this.activeNames);
        layer.setStyle(style);
        var closeTooltip;
        layer.on('mousemove', function(e) {
            var layer = e.target;
            var html = utils.constructTooltip(_this.options,
                layer.feature.properties.name,
                monthData.date,
                monthData[y],
                monthData[x],
                monthData[ratio]);
            this.popup.setLatLng(e.latlng);
            this.popup.setContent(html);
            if (!this._map) this.popup.openOn(this.map);
            window.clearTimeout(closeTooltip);
            if (!L.Browser.ie && !L.Browser.opera) {
                layer.bringToFront();
            }
        }, _this);
        layer.on('mouseout', function(e) {
            var _this = this;
            closeTooltip = window.setTimeout(function() {
               _this.map.closePopup();
            }, 100);
        }, _this);
        layer.on('click', function(e) {
            this.map.fitBounds(e.target.getBounds());
        }, _this);
    });

Also, it's quite jerky on mouseover. Is there a way to debounce the events?


texts: Digitizing and adding attribute values to the points in QGIS

texts:
I want to build a digital terrain model from nautical charts. I need to capture the points where depth is indicated in the chart, and for each point, I must attribute the depth value. One by one, to be sure I capture all points and give the right value. And left a marker on the points to indicate which was already digitized. Below an example of a nautical chart. Blue is water, light brown is land, numbers are depth.

I found tools to digitize points, lines and plots. But didn't find how to attribute values manually to the digitized points.


texts: WFS Layer - capabilities document contained no layers

texts:
I use Quantum GIS 1.8.0 on Windows 7 Professional, 64 bit.  I have problem with WFS Layer. I tried tree WFS:

1) Download services - Parcels (Czech Republic) http://services.cuzk.cz/wfs/inspire-cp-wfs.asp?service=WFS&request=GetCapabilities

2) Geonames (Czech Republic) http://geoportal.cuzk.cz/WFS_GN/WFService.aspx?SID=27082012_03561414729376000_PM_021892838861&service=WFS&request=GetCapabilities

3) Example geoserv.weichand.de:8080/geoserver/wfs?service=WFS&request=GetCapabilities

Example #3 is OK.

Example #2 shows Error "The layer geoportal.cuzk.cz/WFS_GN/WFService.aspx?SID=27082012_03561414729376000_PM_021892838861&SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=gmgml:GEONAMES&SRSNAME=EPSG:102067 is not a valid layer and can not be added to the map"

Example #1 shows Error "capabilities document contained no layers."

I don't know, where are the differences between services. Could you help me?

Thanks Petr


texts: JSON Error? - Unfolding library - Bug or something else?

texts:
While trying to make use of GeoJSONReader from the Unfolding Library, while compiling I got the following error:

A JSONObject text must begin with '{'
JSON Error: Expected a ':' after a key
org.json.JSONException: Expected a ',' or '}' at 28 [character 29 line 1]
    at org.json.JSONTokener.syntaxError(Unknown Source)
    at org.json.JSONObject.<init>(Unknown Source)
    at org.json.JSONObject.<init>(Unknown Source)
    at de.fhpotsdam.unfolding.data.GeoJSONReader.loadDataFromJSON(Unknown Source)
....
A JSONObject text must end with '}'
JSONObject["features"] not found.
JSONObject["features"] is not a JSONArray.
Exception in thread "Animation Thread" java.lang.RuntimeException: java.lang.NullPointerException

And it goes on...
First, I thought it had to do with the JSON itself so I tested via several JSON web validators and it looked fine. I even compared it to another JSON that I knew was working when I used it. And actually,there is no problem with it as it can be displayed on geojson.io, as in the SS.

Now, I'm not sure whether I did something wrong but I'm pretty sure this is the way to load up a GeoJSON/JSON.

List<Feature> busLine = GeoJSONReader.loadDataFromJSON(this, "mures/linia4_linie.geo.json");

Link to the json gist

Any ideas?

Specs: using Unfolding library from here in Eclipse Mars, JDK 1.8.0_91

I've used it more than several times already but never had this issue before.


texts: Display only markers which are near user's current location in Leaflet?

texts: I am using Leaflet js to detect current user location from the browser. I have bunch of markers being displayed in a given area. I would like to show only markers that are closest to the user's current location.


texts: Problem with identify task in ArcGIS JavaScript API 2.8

texts:
I found interesting issue while using the “identify task” tool to display information at a point.

Situation: I am displaying records by Type (A, B, C etc). When I display Type "B" records and I click on anywhere on map, it shows attributes of other records (either "A", "C", or "D")

How is this possible?

Please see the screenshots for further illustration. I also have added the code to JS Fiddle (I changed the server name).

Is this a bug, showing attributes of a hidden record?

Please check this and let me know. Thanks in advance : )

The first screenshot shows all features and an "identify" popup. The second shows only "C"-type records with an identify popup at a location of a different record type.

Update:

I am using single layer and using the combo box to display particular layer. I researched on Google and ESRI forums

Looks like its bug and its ID  NIM042441 also it can be resolve according to ESRI forum so working on the same issue..

So thanks to all : )


texts: Cannot generate subsets of feature class with arcpy (ArcGIS library in Python 2.7)

texts:
I'm having a hard time here on processing GIS data in Python, using library arcpy.
I've been trying to generate independent features from a feature class based on a field of the attribute table which is a unique code representing productive forest units, but I can't get it done.
I've already done this in other situations, but this time I don't know what I am missing.
Here is the code and the error I get:
# coding utf-8
import arcpy

arcpy.env.overwriteOutput = True
ws = r'D:\Projeto_VANT\SIG\proc_parc.gdb'
arcpy.env.workspace = ws

talhoes = r'copy_talhoes'

estados = ('SP', 'MG')
florestas = ('PROPRIA', 'PARCERIA')

arcpy.MakeFeatureLayer_management(talhoes,
                                          'talhoes_layer',
                                          """ "ESTADO" IN {} AND "FLORESTA" IN {} """.format(estados, florestas),
                                          ws)
arcpy.FeatureClassToFeatureClass_conversion(in_features = 'talhoes_layer',
                                            out_path = ws,
                                            out_name = 'talhoes1')
talhoes1 = r'talhoes1'

arcpy.AddField_management(talhoes1, 'CONCAT_T', 'TEXT')
arcpy.CalculateField_management(talhoes1, 'CONCAT_T', """ [ESTADO] & "_" & [CODIGO] & "_" & [TALHAO] """, 'VB')

with arcpy.da.SearchCursor(talhoes1, ['CONCAT_T', 'AREA']) as tal_cursor:
    for x in tal_cursor:
        print(x[0] + " " + str(x[1])) # This print is just to check if the cursor works and it does!
        arcpy.MakeFeatureLayer_management(x,
                                          'teste',
                                          """ CONCAT_T = '{}' """.format(str(x[0]))) # Apparently the problem is here!
        arcpy.CopyFeatures_management('teste',
                                      'Layer{}'.format(x[0]))

Here is the error:
Traceback (most recent call last):
  File "D:/ArcPy_Classes/Scripts/sampling_sig.py", line 32, in <module>
    """ CONCAT_T = '{}' """.format(str(x[0])))
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 6965, in MakeFeatureLayer
    raise e
RuntimeError: Object: Error in executing tool

Can anyone help me out?


texts: Empty or strange cost path output in ArcGIS 10.1

texts:
I´m trying to map green corridors between different stands of hardwood, using Spatial Analyst Cost path tool. I think my cost distance and backlink raster are OK, but the cost path output normally gives a default message, saying in the Results that its empty.
In the map view, all that exists are a few raster cells, on and around the input feature/raster data, sometimes in a gradient from green (good, low costs) in the north and red (high costs) in the south.
Attached some pictures (upper: cost distance, middle: backlink, lower (cost path) to show what it looks like. This specific cost path result was however different, but still useless.

Anyone got solution for this?


texts: Is is possible to have a sld with a rotated point symbolizer along a line in geoserver?

texts:
Geoserver 2.1.3

I would like to symbolize a line with a Square but use the <Rotation>45</Rotation> tags to make the symbol look like a diamond.

My sld is below and unfortunately the symbols are not rotated...Is this not possible or is my SLD wrong?

thanks,

Rob

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
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">
  <NamedLayer>
    <Name>line-symbol rotation</Name>
    <UserStyle>
      <Title>line-symbol rotation</Title>
      <FeatureTypeStyle>
        <Rule>
          <Name>line-symbol rotation</Name>
          <LineSymbolizer >
            <Stroke>
              <GraphicStroke>
                <Graphic>
           <Mark>
             <WellKnownName>square</WellKnownName>
             <Fill>
               <CssParameter name="fill">#009900</CssParameter>
             </Fill>
           </Mark>
           <Size>12</Size>
           <Rotation>45</Rotation>
         </Graphic>
              </GraphicStroke>
            </Stroke>
          </LineSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: GrassGIS process for DTM of UAV point cloud

texts:
Using grassGIS I am attempting to produce a DTM (remove buildings and vegetation) from the las point cloud. Unfortunately the point cloud exported from dronedeploy.com does not distinguish first and last pulses so the method based on GrassGIS v.lidar.egdection does not work.

I need an alternative method/tool particularly for removing buildings?
Any recommendations to transform the z-values relative to the drone take-off location to absolute ones (e.g.: sea-level,..)?


texts:
In the course of a monitoring exercise for small scale farmers, we captured the coordinates of fields using a GPS tracking in the SurveyCTO app. We used automatic capturing of the coordinates every 5 seconds to compensate for low GPS accuracy. Therefore, the coordinate list for every field is very long. With these coordinates, we would like to create polygons in QGIS to e.g. determine the size of fields. Unfortunately, the coordinates can only be extracted in one cell of an excel file. This looks like this:

Now my question is: Does anyone have a suggestion how we can import such a list into QGIS in a way that a polygon is created for each of these coordinate lists?


texts: Identifying projection for shapefile with unknown coordinate system, no z values and no measures?

texts:
I am attempting to create a GIS Project with multiple Brazilian cities. I have been collecting, and am interested in, shapefiles that delineate the neighborhood boundaries of said cities.

One city's shape file that I have acquired is for Curitiba, Paraná (PR), Brazil. I obtained the city-neighborhood shape file from here: Institute of Research and Urban Planning: Curitiba. You will find the file under "Divisa de Bairros" (Neighborhood Divisions).

I am using Quantum GIS to overlay these city shape files over Google Maps, because ultimately I may need to incorporate street or highway information in to my analysis. I am ambivalent ultimately about which source I use, and am open to solutions using Open Streets if possible.

I have read the GIS help files and followed the process to view the metadata in ArcMap to see what the coordinate system and/or projection system is for this shapefile. However, on reviewing the metadata it says that the coordinates have "No Z Values" and "No measures" and there is not defined coordinate system.

My question is: How can I fix this coordinate/measurement issue so that I can correctly project the location of the Curitiba city shapefile on either the Google Maps or Open Streets projections so that I may move forward with my project?

Note: I have already contacted the creator(s) of the files and asked if they may provide any additional information. I will post or answer my own question as need be.


texts: How to add legend to Mosaic Dataset (ArcGIS Image Service Layer)

texts:
I have snow cover data of MODIS (46 images for an year) I put them in a mosaic dataset (with in a file geodatabase) and published them as an Image Service from ArcCatalog and added a time slider.

I am using the ArcGIS Javascript api to make a little web app. Now the last step I'm stuck at is making a legend. (plus I'm a beginner at javascript!)

I read over here ( https://developers.arcgis.com/javascript/jsapi/legend-amd.html) that the legend class supports some specific layer types. I believe the work I have done above is not supported by it? Because I have tried the normal legend making functions but to no use. Can anybody guide me how to make the legend?


texts: Counting average of all polygons instead of every polygon in PostGIS?

texts:
I'm facing a problem with counting avg of all polygons in pgsql query. Instead of that I receive avg(of an area) of every polygon separately. My syntax is looking like this:

CREATE OR REPLACE FUNCTION function(geometry, integer)
RETURNS geometry AS $$

DECLARE

c_density integer := avg(area(geom, 10));

BEGIN
RETURN(
SELECT ST_CollectionExtract(ST_Collect(points),1) as points FROM
(SELECT ST_Collect(function2($1, c_density), ST_Centroid($1)) as points) as points
);
END;
$$
LANGUAGE plpgsql;

I want c_density to be constant. For example I have dataset with areas of:
    10 20 30 40 50 60 70 80 90 100
And I want to c_density to be 55 for every polygon. Instead of that it is like 10, 20, 30...etc. for every next one.

Function area is simply making square root from ST_Area for polygons.


texts: Misplaced Google Maps with OpenLayers plugin in QGIS?

texts:
I was using Google Maps with OpenLayers plugin in QGIS with my vector layer with projection EPSG:4326 with google map projection EPSG:3857 and "on the fly" transformation enabled.

It was matching all the vector layer with Google Maps but now it's not working.

How can I reproject Google Maps so that it matches with my vector layer as before?


texts: How to generate a CSV file from Overpass with ID and tag value?

texts:
I am trying to use http://overpass-turbo.eu ... Something as

[out:csv(::id, wikidata)];
way["wikidata"];
out center;

PS: the answer can restrict area to a country by a kind of area["ISO3166-1"="BR"][boundary=administrative] filter... Or more specific administrative area.

A table with the pair of osm_id and wikidata_id  of each tagged way of the map... And I need also to limit (use a SQL-like LIMIT clause) to test some other variants before request full query.

NOTE 1 - Other solutions

Perhaps taginfo.openstreetmap.org also can generate this CSV, but I not see how to.

NOTE 2 - Overpass answer at comments

I am testing regegular expression but was ignored

[out:csv(::id, wikidata)][timeout:25];
(
  way["wikidata"~"Q[0-9]+"]({{bbox}});
);
out body;
>;
out skel qt;

The answer was  posted by @mmd as comment http://overpass-turbo.eu/s/A84

[out:csv(::id, wikidata)][timeout:25];
way["wikidata"~"Q[0-9]+"]({{bbox}});
out 10;

where out 10 is the answer for the "LIMIT clause" problem for tests or sampling...

NOTE 3 - Other solution for non-limited CSV

Overpass will not run without a "limit parameter", so @mmd suggested other solution, that is not Overpass... I agree.

PS:  the answer can also use terminal and some tool as the suggested by @mmd,  as osmium-tags-filter is perhaps a better solution.


texts: Is there touch screen support in QGIS for Windows?

texts:
I have a new 2in1 PC (laptop and tablet).

It would be very nice to navigate in QGIS using its touch screen I.e. zoom in and zoom out using double finger touch, like on a mobile phone, etc.

Maybe it is already available?


texts: Emitting signal when project property is changed using PyQGIS

texts:
QGIS API provides QgsExpressionContextUtils to set project properties.
But in the class, there isn't any signal.
I hope to connect to a signal whenever some of the properties were changed.
So is there a way to implement this?


texts: Filtering image collection - only one image from a certain date

texts:
My code generates an image collection. In the collection, there are multiple images from the same date (probably from different tiles). Is there a way to filter it and get only one image from a certain date? It's not important which image stays in the collection.
my code:
var image_collection = "COPERNICUS/S2_SR";//image collection.
var first_date       = '2019-06-01'; // First date in image collection
var last_date        = '2019-10-10' ;// Last date in image collection

var fc = ee.FeatureCollection(user_account+folder+shpfile_name).sort('Spatial_ID');
var S2 = ee.ImageCollection(image_collection)
        .filterDate(first_date, last_date)
        .filterBounds(fc)
        .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 18);

//-----------------------------------------------------------------------------------
print(S2);
print(fc)  ;


texts: CartoDB: Adding another view of layer using Layer Actions and Toggle Buttons

texts:
I would like to add multiple sublayers of the same table attached to a single layeraction in cartodb js. I written my webpage based on the following cartodb tutorial: http://docs.cartodb.com/tutorials/toggle_map_view.html. Essentially I would like to attach more than a single query to each button, something where you would populate sublayers[0] and [1] perhaps? Does anyone now how to do this?

var LayerActions = {
  all: function(){
    sublayers[0].setSQL("SELECT * FROM ne_10m_populated_places_simple");
    return true;
  },
  capitals: function(){
    sublayers[0].setSQL("SELECT * FROM ne_10m_populated_places_simple WHERE featurecla = 'Admin-0 capital'");
    return true;
  },
  megacities: function(){
    sublayers[0].setSQL("SELECT * FROM ne_10m_populated_places_simple WHERE megacity = 1");
    return true;
  }
}


texts: Leaflet - tell whether an object a layer or a group

texts:
I need to close popups on an object, but I don't know whether it is a layer or a geoJSON (and hence a feature group). If the former, I can do layer.closePopup(), but if the latter, I have to do layer.eachLayer().closePopup() (I think).

How do I test whether an object in Leaflet is a layer or a group? I'm getting some joy with if (layer.childNodes), but it doesn't always seem to work, and I'd rather have a Leaflet solution than a DOM one, as I think it would be more reliable.


texts: Choosing Geographical Coordinate System and Projected Coordinate System for India?

texts:
I want to create spatial data for entire India.

Usually I use WGS84 as GCS and UTM zones for PCS in the places in India, but when I have to create the spatial data for entire India, I get confused which geographical coordinate systems and projected coordinate system will be suitable.

What are suitable ones?


texts: Flat Earth approximation: forward and reverse formulas

texts:
For very small distances (< 100 m) I need to go back and forth between "polar coordinates" (i.e. bearing and distance) and "equirectangular projection".

I'm not sure whether I understand them correctly. Here what I've done:

double refLat, refLon; // reference (initial) coordinates in radians
double bearing, distance; // target definition

double x = distance * cos(bearing);
double y = distance * sin(bearing);

double targetLon = refLon + x;
double targetLat = refLat + y;

now I need to calculate the bearing and the distance from my current position to the target:

double currLat, currLon; // current position
double targetLat, targetLon; // target position

double meanLat = (currLat + targetLat) / 2;
double dLon = (targetLon - currLon) * cos(meanLat);
double dLat = (targetLat - currLat);

double distance = sqrt(pow(dLon, 2), pow(dLat, 2)) * 6371.0;
double bearing = atan2(dLat, dLon);

I'm not sure how to take in account the cos(meanLat) and Earth's radius in the first function.


texts: Seeking API provider for getting coordinates automatically?

texts:
I need to get the coordinates of some accommodations automatically (and not by extracting them manually from, e.g., Google Maps) and then import them in QGIS. I've already seen that TripAdvisor has an API to download data, but a key is required and it doesn't seem easy to obtain one, especially for academic research.

Do you know if there are other providers that offer this service for free and without the need of a key, with the same quality of contents like TripAdvisor or Maps?


texts: Unsupervised pixel classification in ArcGIS Pro 2.4.0: ISO cluster will not let me classify image into three classes

texts:
I am using ArcGIS Pro 2.4.0, and am attempting to classify a field into a schema with three categories: Leaves, Shadows, and Soil. Here is a sample of the original RGB image:

Through supervised classification (SVM), I was able to create training samples and classify the image into the the desired categories:

Note: Green = leaves, black = shadows, brown = soil. Despite not having the accuracy statistics generated yet, it looks very accurate.

When I try to do the same thing with an unsupervised pixel-based classification (ISO is the only option on ArcGIS Pro that I am aware of), it will not let me divide it into three classes. Instead, it only gives me two:

The only setting I changed from the default ISO cluster settings was the maximum number of classes. I changed that from 5 to 3:

Does anybody know why this would happen? I don't have as much experience with unsupervised classifications as I do with supervised, so I'm not sure why I am not able to get 3 classes. Do I need to change some additional settings?


texts:
Bonjour,

Peut-on avec Quantum GIS redresser une photo aérienne avec des points caractéristiques connus en 3D (XYZ) ?

-- Translation:

Hello,

Is it possible to rectify an aerial photo with known control points in 3d, in QGIS?


texts:
I cannot make ogr2ogr.exe work. I'm using GDAL 1.11.

I keep getting the following error message:

  ERROR 1: Did not get at least 2 values or invalid number of set of
  coordinates 709993.98 6181946.65 32.36 709990.65
  6181946.07
  32.36 709992.47 6181935.8 32.36 709995.78 6181936.38 32.36 709993.98 6181946.65
  32.36 ERROR 1: Invalid exterior ring

I run the following command:

  C:\Python34\Lib\site-packages\osgeo\ogr2ogr.exe -skipfailures -nlt
  MULTILINESTRING25D -dim 3 -f "DXF" -a_srs "EPSG:25832" -spat 711000
  6168300 713000 6168600 -geomfield posList bygning.DXF Bygning.gml

I'm working on a relatively big file 4.2GB.

UPDATE: Example of data to import

<gml:surfaceProperty>
    <gml:Surface srsName="EPSG:25832" srsDimension="3">
        <gml:patches>
            <gml:PolygonPatch>
                <gml:exterior>
                    <gml:LinearRing>
                        <gml:posList>
                            710820.22 6180637.52 36 710818.14 6180633.83 36 710825.47 6180629.71 36 710827.55 6180633.4 38.25 710828.6 6180637.3 38.25 710825.18 6180638.36 38.25 710824.22 6180635.27 36 710820.22 6180637.52 36
                        </gml:posList>
                    </gml:LinearRing>
                </gml:exterior>
            </gml:PolygonPatch>
        </gml:patches>
    </gml:Surface>
</gml:surfaceProperty>

Does the GML3 driver in org2ogr support surfaces?

UPDATE 2 - making it simple

I tried making a gml consisting of only one feature, same namespace etc. as before. I then output the feature to WKT.

Output WKT

"POLYGON ((711383.39 6182528.23,36.49 711379.31,6182528.35
36.48,711379.06 6182520.08,36.48 711383.15,6182519.97 36.48,711383.37 6182527.56,36.49 711383.39,6182528.23 36.49))"

Input GML
http://pastebin.com/76AGLWH0

Cleary the commas is placed in the wrong positions. The z-coordinate of the first point is parsed as the x-coordinate of the second feature. If I add -dim 3 or -nlt GEOMETRY, it adds zeroes.

BR
Thomas


texts: Using Sum() or something similar with CalculateField_management and FieldLists

texts:
I need some help with using variables with CalculateField_management tool. I've tried several different scenarios but don't seem to have any luck.

What I have is a list of fields containing water consumption data for various months.

In my tables I have fields with the following field names: [JUN_12], [JUL_12], [AUG_12], etc... Each table spans 2 years June 2012 - May 2013, the next table will go from June 2013 - May 2014, etc...

These tables are actually views that resides in our billing system. Eventually, I will have another table that will have [JUN_13], [JUL_13] (incremented a year).

I'm trying to avoid hard coding each column name in my code to make it easier to update the code once a year. As these tables are created throughout the years, I'd like to only have to update the yr1 and yr2 parameters and have the script run. This should be the only change from one year / table to the next...

I am fairly new to python and appreciate any help that I can get.

Here's an example of what I have:

#Variables...
yr1 = "08"
yr2 = "09"

sep = "_"

y1 = sep + yr1
y2 = sep + yr2

m1 = "JAN"
m2 = "FEB"
m3 = "MAR"
m4 = "APR"
m5 = "MAY"
m6 = "JUN"
m7 = "JUL"
m8 = "AUG"
m9 = "SEP"
m10 = "OCT"
m11 = "NOV"
m12 = "DEC"

fclass = EnergyGDB + "\\" + PZSPJoin
f1 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m6,y1))][0]
f2 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m7,y1))][0]
f3 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m8,y1))][0]
f4 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m9,y1))][0]
f5 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m10,y1))][0]
f6 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m11,y1))][0]
f7 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m12,y1))][0]
f8 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m1,y2))][0]
f9 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m2,y2))][0]
f10 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m3,y2))][0]
f11 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m4,y2))][0]
f12 = [str(f.name) for f in arcpy.ListFields(fclass, "{0}{1}".format(m5,y2))][0]

fieldList = [f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12]
fields = ",".join("!{0}!".format(x) for x in fieldList)

#Calculation:
arcpy.CalculateField_management(PZSPJoin, "TotalCCF", "sum([{0}])".format(fields), "PYTHON", "")

There are a few other things that are happening in between the calculation and the variables, however this is where I'm having the issue. I'm trying to add all of the fields and have a total in a new field.

If I replace the variables with the actual field names, the script works just fine; but as soon as I introduce variables, it breaks.

Any help would be greatly appreciated. Thanks!


texts: Editing data with Server Object Extension?

texts:
Is it possible to use an interface such as IWorkspaceEdit within an Server Object Extension to create new features or update data in a geodatabse?

Perhaps this can only be done through a managed feature service, but I don't want to expose editing to the client, but rather control it more precisely and programatically in an SOE.

I tried to do this, but I am getting an exception when I try to call ITable:CreateRow() .

For what it is worth, I was able to get an IWorkspaceEdit interface from a map services table, and start editing.


texts: Combining Python and QGIS out of console

texts:
I was wondering if it's possible to use Python and QGIS out of QGIS environment (namely, avoiding the use of the console). I'm pretty new to QGIS, I've been always using ArcGIS 10.1. What I am looking for is a thing similar to that "arcpy" module, since I work with tons of data at a time. Moreover, I'd like to use a good IDE, any suggestion?
I'm using qgis 1.8, python 2.7.2 32bit, windows 7 64bit.
P.S.: I couldn't really come up with a solution, neither following the instructions found in http://www.qgis.org/pyqgis-cookbook/intro.html#python-applications. Hope someone will enlight me!


texts: Adding an OpenLayers 3 layer parameter to change external layer controls

texts:
I have this JSFiddle that has a couple of layers and then a box after it where my opacity slider would be if all the correct external resources were loaded in.

In the 'Cities layer 2' layer parameters I have changed the parameters as so:

var vectorLayer2 = new ol.layer.Vector({
    source: vectorSource2,
    name: 'Cities Layer 2',
    visible: true,
    slider: false
});

I was hoping by doing so when I initialise my layer tree I could query the layer and if the 'slider' parameter was set to true then it would show the slider, and if it was false then it would hide the slider but I'm not sure exactly how to write this?

This also might not be the best way to go around this so I welcome any other alternative solutions.


texts: How to get CRS of a raster layer in pyqgis?

texts:
I have this problem accessing to the coordinate reference system(crs) of a raster layer map.
I wonder if such thing is possible?
I have used this xml from GDAL open map service to load in a raster layer(crs id 900913) to my stand-alone pyqgis application (http://www.gdal.org/frmt_wms_googlemaps_tms.xml)
It displays the map just fine, but when I tried to draw an icon at a specific lat long position (crs id 4326), the icon appeared weirdly.
I used the QgsCoordinateTransform class to convert the coordinates.
Please advise me what I should do to get it right? I am very new to pyqgis.
Thanks a ton!


texts: How to find markers in selected dragbox openlayers 3?

texts:
Here is a part of code,

draw = new ol.interaction.DragBox({
            /* dragbox interaction is active only if alt key is pressed */
            condition: ol.events.condition.altKeyOnly,
            /* style the box */
            style: new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: [0, 0, 255, 1]
                })
            })
        });
        /* add the DragBox interaction to the map */
        map.addInteraction(draw);

        draw.on('boxend', function (evt) {}

I googled but didn't find how will I do. I created a dragbox and added to map , I have some markers on the map. How I learn that dragbox contains markers or not ? Thanks for responses.


texts:
I want to change the a raster's CRS, because i want to fit it with the others layer and rasters.

When i load it, appears a text "The CRS hasn't been defined: will be used as CRS EPSG 4326: WGS 84".
So the raster appears shifted respect to the others.

How can I solve this problem?


texts: Difference between haversine formula & geodesic distance for calculating distance

texts:
Background

I have Vehicle driving data which includes Lat/Long points of a particular route. I drive in that route for 2 times.
The driving data in that route is around 20 minutes. (Approximately for 2 times = 40 mins)

Distance Calculation

Before I have been using haversine formula to calculate distance between every point between route 1 & route 2.
recently I came across geopy library which uses geodesic distance function to calculate distance.

Question/Requirement

So my question is, which one produces better results either haversine  or geodesic distance?


texts: How to identify types of polygons based on line segments

texts:
I have a mission that i have to identify several types of polygons. I got stuck at two types of categories that I don't know how to identify. I will post every situation here, maybe someone has an idea.
I have a layer with polygon segments (line type) which has a calculated field where the front street is equal to 1 and the others to 0. I also have the field with the number of faces.
I have to identify on this layer with line segments that also have the group field with identifier for each polygon. Then it's easy for me to turn this layer back into a polygon layer, but I don't know how to identify 2 types of plots.
Case 1 - angle plots - at the crossroads of two roads.

Here i have to identify the consecutive segments that have the value 1 in the " front street" field.
Case 2 - polygons with 4 sides, which has two sides overlooking two different roads.

Here i have to find the polygons that have 2 parallel street segments.


texts: Is it possible to convert a jpg file into lyr

texts:
I am a beginner with GIS 10.0. IS there any way to convert a jpg file into a layer?
I am working on a problem that needs to have the entire map of my city as a layer but I just have a jpg file of the map.


texts:
How do I use the ogr SetAttributeFilter with uppercases?

Code example:

lyr = shp_obj.GetLayer()
sel_str = "%s = '%s'" % (field, value)
lyr.SetAttributeFilter(sel_str)

Can I just do:

 sel_str = "UPPER(%s) = '%s'" % (field, value)

?


texts: Getting matrix from polygon attribute table that represents vertex couples using ArcGIS Desktop?

texts:
How can I get a matrix from a polygon attribute table that represents polygon vertex couples?

What I want is something like this :

A = [ (A1,B1) , (A1,C1) , (B1,C1)
      (A2,B2) , (A2,C2) , (B2,C2)
      ....    ,  ....   ,    ....
      ...     ,   ....  ,  ....  ]      ]


texts: How to programmatically interrupt QGIS initialization and startup

texts:
I'm wondering how to interrupt the initialization of QGIS 2.0 on Windows XP before start up. I've put a dialog in one of my plugins that shows up on start up. However, I'd like to make a quit button that closes QGIS if desired by user.

I've tried two ways without success:

This gives me an error for loading the plugin due to sys.exit(), so it doesn't successfully close QGIS:

result = self.dlg.exec_()
if result == 1:
    self.doStuff()
else:
    sys.exit()

This causes QGIS to crash--which successfully closes the program, but really not in the way I'd like it to:

result = self.dlg.exec_()
if result == 1:
    self.doStuff()
else:
    QgsApplication.exitQgis()

I've been wondering about this second situation. In the case that I don't find a better solution, what risks may I be taking if I let QGIS crash repeatedly on purpose?

I've taken a look at this link, but it only interrupts the start up until a user succeeds to authenticate. I'm trying to interrupt start up by choice of user who clicks a "Quit" button before QGIS even finishes setting up.


texts:
I have a set of locations/points and I am trying to compute the nearest distance of these points from a Polyline. The points and the polyline lie in China and the Geographic Coordinate System is WGS 1984. I have some questions regarding this.

Since my task involves computing distances, I am trying to use an equidistant projection. The following pdf suggests that the best equidistant projection for China is the Polyconic projection.

http://rcg.gvc.gu.se/data/ChinaPrecip/_notes/projection.pdf

The ESRI Polyconic projection page indicates that the parameters I will have to set include False Easting, False Northing, Central Meridian and Latitude of Origin.

http://resources.arcgis.com/en/help/main/10.1/index.html#//003r0000003t000000

What would be the best value for these parameters over China?

Also, I previously thought that the selection of central meridian and latitude of origin does not impact distortions. The ESRI link above, however, under the "Properties" and "Distance" section says "The scale along each parallel and along the central meridian of the projection is accurate. Distortion increases along the meridians as the distance from the central meridian increases.
" This seems to imply that the selection of central meridian impacts distortions or am I understanding this incorrectly?


texts: GeoServer - Ingest not supported for mosaics

texts:
GeoServer 2.18.0, Java 1.8 on Windows.
I have succeeded in creating an ImageMosaic store using the REST API.
This store is located as a single directory in the file system. The data is GeoTIFF.
I would like to add new tiles to this store. To do this I add a targetStore JSON object to the REST API import object. GeoServer returns an INIT_ERROR with a message that "ingest not supported for mosaics".
This seems to indicate that it is not possible to extend the mosaic once created. Is this the case, or maybe I am calling the API incorrectly?


texts: Testing if ArcPy script will work with lower license level?

texts:
I was looking at the question Exporting feature class to SQLite database using ArcPy with ArcGIS Pro? which states at the end

  The problem is that I cannot use feature class to feature class
  because in basic licence I cannot copy feature from geodatabase to
  sqlite database.

I couldn't find documentation to back this up, so I thought I'd test it.  My thinking was that I would be able to change my license level down to Basic and try it out, however the ArcGIS Administrator wouldn't let me do that because I don't have a Basic license (I have Advanced license only).

So I am unable to confirm that the tools won't work with a Basic license because I can't limit the license I have to Basic level.

This got me thinking.  How do you test that an arcpy script will work at a lower license level than the one you're currently using if you don't have a lower license available to test with?

I read that you can import the lower level before importing arcpy, but this fails for me.

>>> import arcview
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcview.py", line 18, in <module>
    gp.setProduct("ArcView")
RuntimeError: ERROR 999999: Error executing function.
>>> import arceditor
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arceditor.py", line 18, in <module>
    gp.setProduct("ArcEditor")
RuntimeError: ERROR 999999: Error executing function.
>>> import arcpy
>>> arcpy.ProductInfo()
u'ArcInfo'
>>>


texts: Getting ArcMap Python Parser to wrap text displayed in MapTips?

texts:
I have an attribute field [NOTE] that contains very long strings - often 255 characters - usually consisting of multiple sentences.

I'd like to be able to use ArcMap's MapTips with some form of wrapping, because without it the MapTip text usually runs off of the map pane, obscuring important parts of the note.

Is it possible to do this in the Advanced box using the Python parser? I

tried variations of a python function textwrap.wrap([NOTE]) to no avail.


texts:
The following code works very well in pyqgis 2.0 with a line shapefile, but when executed with pyqgis 3.0 the instruction measureLine return 0 and with no errors signaled.
What is the new instruction in pyqgis 3.0?

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

def dist_seg():
  layer=iface.activeLayer()
  L=[]
  features = layer.getFeatures()
  d = QgsDistanceArea()
  for i,element in enumerate (features):
    m = d.measureLine(element.geometry().asPolyline())
    L.append ((i,round(m,1)))
  return L


texts: Export Raster Data to TIFF, BMP, or PNG File using ArcPy?

texts: I am writing a Python script to export a Raster Dataset to a TIFF, BMP, or PNG file, as I would manually by right-clicking on the layer, choosing Data | Export Data, and then make my choices in the Export Raster Data window.  Can this be duplicated with ArcPy, and if so, how?


texts:
I am georeferencing some tiff images in QGIS using the same settings for the goereferencer's parameters.
A couple of images were referenced fine but for others the georeferenced output raster is weirdly incomplete and huge in file size.
details of files processed fine:

initial file size: just over 99 Mb
referenced file size: between 600-800 Mb

Details of problematic files:

initial file size: just over 290-350 Mb
referenced file size: > 4Gb

Main Problem:
Output file size a part, when I visualize the georeferenced raster of the problematic files I see that the image is not complete: I can see it has been rotated (as usual) but is has also cut.
On a side note the referenced image is not automatically added to the project but i need to add it.
I'm working with QGIS 3.10.6 on a Windows 10 machine.
the settings I'm using in the georeference are:

Transformation type: Elmerth
Resampling method: Nearest Neighbour
CRS: EPSG 32619
compression LZW
all images referenced using 2 points (please note that this protocol worked ok with previous images, in this session as well as in the past).

Here some screenshot to see what I mean:
Setting used in all Image refrerencing. In particular the screenshot is from one problematic image (you can see it whole image was imported in the groreferencing tool).

output:

on the left the problematic image on the right the ok one.
Any idea of what is causing this and/or workarounds?


texts: Changing point data coordinate system without moving points in QGIS?

texts:
I made a mistake regarding coordinate systems and lining up two layers.

I have a raster with elevations in EPSG:2949. I created a point data shapefile and entered about 600 hundred points manually lining up with where I want them on my raster. After adding all those and trying the Point Sampling tool I found that my point data shapefile has EPSG:4326 as it's coordinate system so Point Sampling will not get me an output.

Switching my coordinate system moves my point elsewhere as expected so basically I'm trying to get my two layers (raster and point data) on the same coordinate system while keeping the points where I drew them originally.

Is that possible?


texts:
I am trying to implement some TMS specs as per http://wiki.osgeo.org/wiki/Tile_Map_Service_Specification. I am particularly interested in the service discovery part ( and ).

The page is very clear wrt status of the specs, ie experimental. I am wondering if it is worth implementing it at all in fact.

Bottom line is that I am struggling to find a server that would actually comply to the specs above mentioned.

Does anyone know of such a server? + state of the specs? should I just stick to the good old url pattern and drop all the server service discovery stuff?

Thanks for any pointers,

Thomas


texts: Intersecting to vector layers when there is no common attribute

texts:
I am trying to join two vector layers into one, the first is a road noise barrier shapefile, the second is a vectorised DEM.
The road barriers have a relative height attribute for each section of road, one barrier passes through several DEM regions while maintaining the same relative height.

When running the intersect operation I get the following error: GEOS geoprocessing error: One or more input features have invalid geometry.

How can I fix this error or get around the problem another way to find the height of the barrier after subtracting the Local DEM.

Changing the CRS didn't solve the problem


texts: How to parse polygons from GML in OpenLayers?

texts:
I'm reading a GML file from an URL. This file contains polygons. I would like to draw the polygons on a map, where each should have a different color. But all examples I found on the web so far seemed to be outdated and didn't work. Any Hints?

Here the parsing (that's the part that already works):

    var theParser = new OpenLayers.Format.GML();
    theParser.internalProjection = null;
    theParser.externalProjection = null;
    theParser.extractStyles      = false;
    theParser.extractAttributes  = true;
    var features = theParser.read(response.responseText);
    if(features)
    {
         for(var i=0;i<features.length;i++)
        {
            var feature = features[i]
       //Each feature is an instance of OpenLayers.Feature.Vector
           //Each has an attribte geometry, which is of type
           //OpenLayers.Geometry.MultiPolygon and which contains an
           //OpenLayers.Geometry.Polygon
        }
}


texts: QGIS installation from OSGeoW4: entry point not found in qgis_core.dll

texts:
I have installed QGIS through OSGeo4W installation. My operation system is windows 7 64x. I try to start the program QGIS Desktop from start menu this error:

Do I need to do something? Installation was successful and re installing doesn't change anything.


texts: Use a field other than the category field to set marker-width in CartoDB's Torque Category

texts:
I have rabies data that contains a date, animal type and county.  The centroid of the county polygon was used to generate a geometric point, the date used for Torque's Time Column, and the animal type for Torque's Category column.  I would like to use the count of each type of animal for a given date county combination to set the marker-width using CartoCSS.

It seems that for Torque the value of a field is not accessible to CartoCSS as [field] although you can access the numeric value of your category using [value].  Is there a way to pass the value of another field to Torque either using sql to generate count(animal type) grouped by animal type, time and county or by having Torque directly calculate the number of category records?