texts: Saving Shapefile to File Geodatabase in QGIS 2.2?

texts:
I searched the StackExchange forum for an answer too my question and apologize if I missed it.

I'm using QGIS 2.2 64-bit and able to read feature classes from an ESRI file geodatabase without a problem.  The problem occurs when I attempt to save a layer to a file geodatabase by right-clicking on the layer and selecting "Save As."  I have "ESRI FileGDB" selected as the desired format and enter a name for the exported feature class.  The result is a file geodatabase created with the name of the feature class and an error message

  Export to vector file failed.

  Error: creation of layer failed (OGR error:Error: Failed at creating
  FeatureDataset  (General function failure.))

Any idea why the feature class export to a file geodatabase fails?


texts: MySQL performance issue with ST_Contains not using spatial index

texts:
We are having what seems to be a fairly large mysql performance issue on trying to run a fairly simple update statement. We have a table(1.8mil) with houses that contains a Lat+Long geometry point column(geo), and then a table(6k) that has a list of schools with a boundary geometry polygon column(boundary). We have spatial indexes on both, we are trying to set the school's id, that contains the point, to the house table with the update. The update is taking 1 hour and 47 minutes to update 1.6mil records. In other systems I have used in my paste experience, something like that would take just a few minutes. Any recommendations?
I have posted this same question in the DBA SE site as well, as it is very much a GIS & DBA question.
CREATE TABLE houses (
  ID int PRIMARY KEY NOT NULL,
  Latitude float DEFAULT NULL,
  Longitude float DEFAULT NULL,
  geo point GENERATED ALWAYS AS (st_srid(point(ifnull(`Longitude`,0),ifnull(`Latitude`, 0)),4326)) STORED NOT NULL,
  SPATIAL INDEX spidx_houses(geo)
)
ENGINE = INNODB,
CHARACTER SET utf8mb4,
COLLATE utf8mb4_0900_ai_ci;

CREATE TABLE schoolBound (
  ID int PRIMARY KEY NOT NULL,
  BOUNDARY GEOMETRY NOT NULL,
  reference VARCHAR(200) DEFAULT NULL,
  type bigint DEFAULT NULL,
  INDEX idx_reference(reference),
  INDEX idx_type(type),
  SPATIAL INDEX spidx_schoolBound(BOUNDARY)
)
ENGINE = INNODB,
CHARACTER SET utf8mb4,
COLLATE utf8mb4_0900_ai_ci;

-- type 4 means it's a elementary
Update houses hs
    INNER JOIN schoolBound AS sb ON ST_Contains(sb.boundary, hs.geo) AND sb.type = 4
SET hs.elementary_nces_code = sb.reference

The explain seems to show that it is not going to use the spatial index for schoolBound.
+----+-------------+-------+------------+------+---------------+------+---------+------+---------+----------+------------------------------------------------+
| id | select_type | table | partitions | type | possible_keys | key  | key_len | ref  | rows    | filtered | Extra                                          |
+----+-------------+-------+------------+------+---------------+------+---------+------+---------+----------+------------------------------------------------+
|  1 | SIMPLE      | sb    | NULL       | ALL  | NULL          | NULL | NULL    | NULL |    6078 |    10.00 | Using where                                    |
|  1 | UPDATE      | hs    | NULL       | ALL  | spidx_houses  | NULL | NULL    | NULL | 1856567 |   100.00 | Range checked for each record (index map: 0x4) |
+----+-------------+-------+------------+------+---------------+------+---------+------+---------+----------+------------------------------------------------+


texts:
Today I found that ArcGis can make colored hillshade like on piture below.

Here is link to their blog. I want to do that in QGIS. I know about layer coloring feature, and it's close for perfect, but it can use just one color. So maybe there is some plugin for QGIS with that I can use color scheme (color ramp) instead one color.

I'm done it! Here is my result:

Not so good like a result from ArcGis blog, but I'm on the right way!
What I've done?
Firstly, I've created a new hillshade layer, with Hillshade tool, as noted in the answer. Then, I've choosed pseudocolor style in the properties of this layer. Of course I've applied color ramp similar to one from the blog. Finally, I've repeated all three layers from the blog, and made them all 50% transparent. Done.


texts:
I am still beginner in Remote Sensing.

I have asked a GRASS programmer to prepare for me Landsat reflectance values from raw DN values (1 scene per year from 1987 to 2013). I have proposed him (following litterature and p.ex. Is there a way to automate the preprocessing of LANDSAT data?) to run step by step next tools to obtain reflectance values:

i.landsat.toar (addon for GRASS 6, included in GRASS 7) - convert DN
to top of atmosphere radiance
i.atcorr - correct top of atmosphere to   surface reflectance

These surface reflectance values were subsequently corrected by topographic normalisation DOS3 by i.topo.corr because of my mountainous AOI.

How is it possible to establish/check if calculated reflectance values are right?

The reason why I am asking is because if the reflectance values (right) are more appropriate for vegetation study.

Why do the raw Landsat data (left) values sees to be better visually distinctive ?

Is it only a question of appearance?


texts: QGIS Raster Merge Outputs Mysterious Bands

texts: I am using the Raster>Miscellaneous>Merge tool in QGIS to merge two different GeoTiff(s) into one layer. They are two different agricultural fields. When I merge them, it seems that the rgba bands are not retaining their values. Is there something I am missing? I understand the black is just the no data value I set. Not worried about that. I need the third image to have the first and second image with that rgba structure.


texts: Sampling values from a .vrt raster using PostGIS points

texts:
I have about 1 mill 2D points in a table in a PostgreSQL database.

I also have a DEM (Digital Elevation Model) in raster. It’s a lot of .tiff files, bound together by a .vrt file, so they function as one raster data set, at least when seen from QGIS.

I’m looking for the best way to attach to each point, the height in the DEM, at the points location. I would prefer to have that point in a new field in the point table (i.e. as a new attribute to the point feature).

After looking at the logical (PostGIS) possibilities, ST_Value  [http://postgis.net/docs/RT_ST_Value.html] seems to be the relevant function. But I can’t see any ways to point the ST_Value() to a .vrt file. And since the DEM is about 1Tb I would like to avoid pulling it inside the PostgreSQL database. I was briefly considering FDW (Foreign Data Wrapper) but have no experience in that direction.

Googling around I found GDAL’s gdallocationinfo  [http://www.gdal.org/gdallocationinfo.html]. So I generated a text file with the 1 mill x y coordinates of the points, and copy-pasted that list into a batch-file with lines like this.

gdallocationinfo -geoloc -valonly my_dem.vrt 444302.48 6157429.94 > collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444322.773687121 6156609.71326993 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444514.567093169 6209236.39142913 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444516.550364779 6209265.20736336 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444523.113049805 6209178.98623131 >> collect_dtm_values.values
gdallocationinfo -geoloc -valonly my_dem.vrt 444528.996558214 6209137.25582758 >> collect_dtm_values.values
...

Though this seems to work, it will estimated take 3 days to run the 1 mill lines batch file (it’s running).

A part from not being an elegant solution, it also leaves me with just a simple text file, holding one number per line. I would now have to figure out a good way to connect the values back onto the points, and the only connection is the fact that their order of appearance in the output text file, is the same as the order of appearance in the input batch file.

I will have to do it again and again, so I’m willing to invest in a better solution :-)

So – the question is: What would be the good way to sample points in a .vrt raster, given a set of points in a PostgreSQL table?

I work in PostgreSQL/PostGIS because I find it fast and convenient, so a solution in this environment would be preferable, and in any case, a solution that runs considerable faster, than what I have now, would be nice.


texts: Restrict Layer Style getting saved in QGIS project table in PostgreSQL

texts:
As postgres admin , I have created a layer and styled the same in QGIS AND saved the style in postgres. When another user uses the same project(user have been given the grant to access the project) and changes the styling(color of the feature) and saves the project , the default style gets changed.
How do I stop the user from saving the project with the style?


texts: QGIS: Smoothing two polygon sets where polygons are nearly adjacent

texts:
I have two polygon sets, one representing cities, the other representing school districts. I am trying to create a relationship between cities and districts, but I have a problem where the cities and districts share a boundary....the boundary line is not consistent and some points technically cross over the boundary, and as a result my ruby script is creating incorrect relationships when a city and school are adjacent.

Is there a way to smooth the datasets so that adjacent polygons share borders that are identical in terms of the lat/lng points along the border?

I have each dataset in a .csv file with lat/lng pairs and a descriptor attribute.


texts: Re-ordering fields permanently using ArcGIS Make Query Table tool?

texts:
In the Answer to Reordering fields (permanently) in file geodatabase using ArcGIS Desktop? it was stated that the Make Query Table tool could be used to re-order fields in tables and feature classes permanently.

However, when I open the dialog for that tool I cannot see anywhere that enables re-ordering of fields.  All I can see is that fields can be carried through or dropped using checkboxes but in the same order.

Assuming that it is possible can someone describe the steps in more detail, please?


texts: ArcGIS - how to join attributes of overlapping polygons

texts:
Given two tables each with unique attribute sets, how do I merge the attributes of one layer onto the other, or create a new layer containing the attributes from both?

The image shows my two layers: the white cells underlie the red cells. I would like the red cells to "sample" the attribute data form the white cells. I've tried to do a spatial join (right clicked first layer) but I get an empty table.

UPDATE

I tried intersect (green layer) but it generates a larger area than the red cells layer.

UPDATE 2 - using centroids


texts: Generating start and end points for LineStrings

texts:
I want to create start and endpoints in QGIS for a Linestring shapefile.
I tried using the "QChainage Plugin" but was not useful.
Can someone suggest a tool?


texts: How can I change the NoData value in QGIS 2.x using the python console

texts:
I'm trying to add a NoData value to a single band raster image in QGIS 2.0 using the python console but can't find the right bit of code to do it.

Here's what I see on the console:

>>> layer = qgis.utils.iface.activeLayer()
>>> layer
<qgis.core.QgsRasterLayer object at 0x06AFFDF8>
>>> layer.renderer()
<qgis.core.QgsPalettedRasterRenderer object at 0x1874E0C0>

Any help here would be great!


texts: Can Leaflet method .getBounds() be used on a layer to determine if a point is within it's bounds?

texts:
I am trying to determine if a set of points are within the bounds of a grid layer (specifically what rectangular cells are those points in) so I can highlight just those grid cells (or add a popup for that grid cell). Using MongodB or a spatial database query is not an option; the solution needs to be within the realm of javascript and Leaflet. I've looked at the leaflet-pip plugin, but I don't see why using .getBounds().contains() to do this won't work using the grid layer instead of the whole view. In the code below, I have created an empty var for the grid (var = L.geoJson();) which should act as a global variable once it is filled through the $.getJSON request (I am not including that part). I am then using it within my ajax request to get the point data. I am currently getting a leaflet-src.js:1873 Uncaught TypeError: Cannot read property 'lat' of undefined error, but I think I am just doing something simple wrong. Any ideas of how I could correct this to get it to work?

$.ajax({
    type: 'GET',
    url: "<url>",
    timeout: 8000,
    crossDomain: true,
    dataType: 'json',
    data: {},
    headers: {
        "VMUser": "user"
    },
    //Renders data to view
    success: function (data) {
      console.log(data);
        data.forEach(function (value) {
            var jobName = value.name;
            var createdBy = value.createdby;
            var latitude = value.latitude;
            var longitude = value.longitude;
            var id = value.jobid;
            var status = value.status;
            var provider = value.provider;

            var marker = L.marker([value.latitude,value.longitude]).addTo(map);

            var circle = L.circle([value.latitude, value.longitude], 500, {
             color: 'red',
             fillColor: '#f03',
             fillOpacity: 1.5 }).addTo(map);

            *if(grid.getBounds().contains(marker.getLatLng())){console.log('marker detected')}*

            grid.bindPopup("Job Name: "+jobName+"
Created By: "+createdBy+"
Latitude: "+latitude+"
Longitude: "+longitude+"
Status: "+status+"
Provider: "+provider+".");
            grid.on('mouseover', function (e) {
            this.openPopup();
            });
            grid.on('mouseout', function (e) {
            this.closePopup();
            });

        });
    },

  });
});


texts: Export multiple Layers at same time to XLSX

texts:
How to export multiple/selected layers to one or more Excel-Sheet?

I do not found an internal function or a plugin for this.


texts: Hiding features of a layer at different scales

texts:
I'm using ArcMap 10.5.1
I have a map with journeys made by patients across the country as lines features. All journeys are in the same layer. Some of these journeys go right across the country, while some are very localised. When I create maps of small districts I dont want to see the journey lines (features) which extend across these districts as they obscure whats going on within the district. I am aware of two potential solutions;

create separate layers out of the features I don't want to include. In my instance, the journeys I do want are not distinguishable from the features I dont want, so to select them by attributes of another feature class is not possible. I would need to do this manually, which would be laborious.
Convert to graphics. This would allow me to delete individual lines, but as I create different map images, this would be clumsy.

I was wondering if there is a way to specify that only lines (within the layer) completely visible within a given scale/map extent, be drawn?


texts: Temporally grouped reduction in Earth Engine: confusing behavior

texts:
I am trying to create two Image Collections in Earth Engine, each representing monthly composites of Sentinel 1 observations in my ROI. The first contains VV images, the second contains VH images.
I've adapted the code posted by Nicholas Clinton in this example.
The code works beautifully for the VH cross-polarized collection, but the same code applied to VV results in an error. The attached image shows the issue: When I create the initial image collection, before mapping the function to create a new collection of monthly composites, everything goes smoothly. Printing the first image in both collections shows that the three expected bands are present (VV, VH, and 'angle'). But when I map the monthly-composite function over both collections, I get an error for VV saying that no such band exists, while VH works fine and the first monthly composite is added to the map, as expected.
I assume I'm making an obvious mistake with my code but for the life of me I can't figure out what it is!
Here is the relevant code in full:
var S1VV = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .filterDate('2014-01-01', '2019-12-31')
        .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
        .map(function(image) {
          var edge = image.lt(-28.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage);
        });
print(S1VV.first())

var S1VH = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .filterDate('2014-01-01', '2019-12-31')
        .filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'))
        .map(function(image) {
          var edge = image.lt(-30.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage);
        });
print(S1VH.first())

// Group by month, and then reduce within groups by mean();
// the result is an ImageCollection with one image for each
// month.
var months = ee.List.sequence(1, 12);
var S1VV_monthly = ee.ImageCollection.fromImages(
        months.map(function (m){
          return S1VV.filter(ee.Filter.calendarRange(m, m, 'month'))
                      .select(1).select('VV').mean()
                      .set('month', m);
        }));
print(S1VV_monthly);

var S1VH_monthly = ee.ImageCollection.fromImages(
      months.map(function (m) {
        return S1_VH.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .select(1).select('VH').mean()
                    .set('month', m);
}));
print(S1VH_monthly);

Map.addLayer(S1VH_monthly.first(), {min: -18, max: -12}, 'VH');


texts: Drawing polygons with specific heights and widths in ArcGIS Desktop?

texts:
I have a dataset that defines my plot corners. Within each plot, there are a number of subplots, each with specific dimensions. I need to create polygon shapefiles for each of the subplots.
Is there a way to do this in the ArcGIS GUI?
In other words, I currently have this -

and I need to generate this -

Any ideas?


texts: Using variables from one Python script in other Python scripts?

texts:
I want to set up a single python file containing variables that have the locations of all my data sources.  This would then be used by all my other scripts, then as a data source changes, I only have to edit the one file.

My Data_sources.py would look something like this:

BC_BEC = "Database Connections\\BC.sde\\FOREST_BC.BC_Eco\\BC.BC_BEC"
BC_TFL = "Database Connections\\Forest_BC.sde\\BC.BC_Admin\\BC.TFL_Boundary"

and so on...

How do I get all my other scripts to load or access these variables in the Data_sources.py?


texts: Compare line with next line in attribute table shapefile

texts:
I am making a python code for ArcGis 9.3. The aim is to compare the value of my field "ID_BV" at line l with value of field "ID_BV" at line l+ 1. And depending of the test answer i populate another field or not.

The problem look very simple : i don't know how can i acces to line l+1. I try to put twice :

rowNext = rows.Next()

rowNext = rows.Next()

With row my line and rowNext the next line.

But it doesn,t work. I also try a lot of other thing like that but it doesn't work.

The entire code is as follow :

import sys, string, os, arcgisscripting, glob
gp = arcgisscripting.create()

# Set a toolbox
gp.AddToolbox("C:/Program Files (x86)/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")

MyShape = "C:\Users\marie.kremer\Desktop\Fabrication_Hydrogrammes\CalculMoyenne\\Test_Inter.shp"

gp.addfield (MyShape, "Moyenne", "double")

#Initialisation de variables
Somme = 3
#Position du curseur dans Mysghape
rows = gp.UpdateCursor(MyShape)
row = rows.Next()
rowNext = rows.Next()
rows.updateRow(rowNext)
rowNext = rows.Next()

gp.AddMessage(Somme)

while row:
    ID_BV = row.GetValue("ID_BV")
    ID_BV_Next = rowNext.GetValue("ID_BV")

    gp.AddMessage("ID_BV")
    gp.AddMessage(ID_BV)
    gp.AddMessage("ID_BV_Next")
    gp.AddMessage(ID_BV_Next)

    if ID_BV == ID_BV_Next :

        gp.AddMessage("Entrée dans le IF")

        Somme = Somme + row.GetValue("AfoisF10p6")

    else :

        gp.AddMessage("Entrée dans le Else")
        Somme = Somme + row.GetValue("AfoisF10p6")
        row.SetValue("Moyenne", Somme)
        Somme = 2

    rows.updateRow(row)
    rows.updateRow(rowNext)
    rowNext = rows.Next()
    row = rows.Next()

Thanks a lot ;)


texts: Anyone know of a map service (WMS, WFS, Esri, other) for the US National Grid?


texts: Geocoding Addresses - Error when selecting Address Locator to use

texts:
I'm running into some issues with ArcMap 10.3.1 when trying to use the geocoding tool to geocode a table of addresses from an .xls file. I've added my table into the Table of Contents, right clicked and chose "Geocode Addresses...". From here, the "Choose an Address Locator to use" dialog box appears. Our group uses a specific address location "Composite_Prod", which I've added by clicking the "Add" button and selecting it from our GIS Server (See screenshot. Everyone else in the office has no problem with this address locator). But, when I click "OK", the error below displays:

"Error: The Microsoft Jet database engine cannot find the input table or query 'GDB_Items'. Make sure it exists and that its name is spelled correctly"

I have no idea where "GDB_Items" came from and at which point this table was selected as an input.


texts: How to get SOAP API Access on an ArcGIS Online Map?

texts:
I published a tiled map service to ArcGIS Online.

This web site reported that we would be able to access the ArcGIS Online maps via the SOAP API from our applications.

I check my first published map and under Supported Interfaces, I have REST but not SOAP options.

Does anyone know how I can can get SOAP access to my ArcGIS Online Map?

Thanks.


texts: Extract cells above certain values of raster file

texts:
I have a very large raster layer and with the values for noise emission. I'd like to extract all cells which are above a certain threshold value (e.g. 50 db) and convert them to an ESRI Shapefile so that I can intersect then with a line shapefile at the next step. The files should be as small as possible. My TIF file is 600MB.
I'd prefer to solve the issue using R.


texts: How to work with Spatial Autocorrelation (Morans I) on equal weighted point

texts:
I have a point feature class. All points are the same, it's not weighted.

I want to use Spatial Autocorrelation (Morans I) tool but it only takes a field for weight. I tried fooling it by adding a field with the value 1 for all points but it comes back with an error.

How can I use Moran I in ArcGIS on a point feature class?


texts: Survey123 Connect-fieldLength doesn't seem to work

texts:
I'm working on my first webform using Survey123 Connect and am struggling to figure out what I'm doing wrong.
I'd like to change the length of some of the entry fields so if I am asking them to enter their age, it doesn't take up the entire width of the page as it does now.
I tried using the bind::esri:fieldLength, assuming that entering a number in that field would limit the entry to that number of characters. Is there something else I need to do to put that into effect?

With those values entered, all the entry fields still take up the entire width of the page. I'm sure I'm missing something very simple.


texts: Unable to download VIIRS images with GIBS API and gdal_translate

texts:
I'm trying to download some VIIRS images from GIBS API. I'm using gdal:

gdal_translate -projwin 19.07 42.16 28.39 34.5 -of "GTiff" "viirs.xml" "VIIRS_SNPP_DayNightBand_ENCC.tif"

The viirs.xml file contains:

<GDAL_WMS>
    <Service name="TMS"><ServerUrl>https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/VIIRS_SNPP_DayNightBand_ENCC/default/2017-08-21/500m/${z}/${y}/${x}.png</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-180.0</UpperLeftX>
        <UpperLeftY>90</UpperLeftY>
        <LowerRightX>396.0</LowerRightX>
        <LowerRightY>-198</LowerRightY>
        <TileLevel>7</TileLevel>
        <TileCountX>2</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:4326</Projection>
    <BlockSizeX>512</BlockSizeX>
    <BlockSizeY>512</BlockSizeY>
    <BandsCount>4</BandsCount>
</GDAL_WMS>

But, it returns the following error:

Input file size is 131072, 655360
ERROR 1: GDALWMS: Incorrect bands count 1 in downloaded block, expected 4.
ERROR 1: GDALWMS: ReadBlockFromFile (https://gibs.earthdata.nasa.gov/wmts/epsg4326/best/VIIRS_SNPP_DayNightBand_ENCC/default/2017-08-21/500m/7/21/88.png) failed.
ERROR 1: viirs.xml, band 1: IReadBlock failed at X offset 88, Y offset 21
ERROR 1: GetBlockRef failed at X block offset 88, Y block offset 21


texts: Python script for dark object subtraction on Landsat 8 Imagery

texts: Does anyone know of any python scripts that will allow me to conduct such a process? I am limited to open source software and have minor experience in Python.


texts: Preventing ArcMap from labeling when creating grids?

texts:
I am using the Create Polygon Grid Wizard to make a large number of grids.

Is there anyway I can stop ArcMap from drawing the grid labels as they are created?

It seems to take much much longer to create the grids because the labels are constantly being drawn. I still want them created but not drawn during the creation process.


texts:
I'm a beginner with cartodb.
Can I overlap layers with cartodb maps from different user accounts?
I mean, I've tried this code and doesn't work, from Rolling Stones Case Study. First, add 'saleiva' user basemap, and then, my other tables, user 'carmenamaya'

var baseLayerDef = {
   sql: "SELECT * FROM rolling_basemap",
   cartocss: "#rolling_basemap { polygon-fill:#333333; polygon-opacity: 0.7; line-opacity:1; line-color: #000; line-width: .3; [feature='Urban Area']{ polygon-fill:#000; polygon-opacity:0; line-width: 0; } }"
};

cartodb.createLayer(map, {
  user_name: 'saleiva',
  type: 'cartodb',
  sublayers: [
    baseLayerDef,
 //   pointsLayerDef,
 //   linesLayerDef
  ]
}).on('done', function(layer) {
   layer created, add it to the map
  map.addLayer(layer);

var pointsLayerDef = {
  sql: "SELECT *, date as date_proc, ST_asGeoJson(the_geom) as geom FROM carmen_amaya",
  cartocss: "#carmen_amaya::oth { marker-fill: #000; marker-opacity: .3; marker-width: 17; marker-allow-overlap: true; } #carmen_amaya { marker-fill: #FFF; marker-opacity: 1; marker-width: 5; marker-line-width: 0; marker-placement: point; marker-type: ellipse; marker-allow-overlap: true; }",
  interactivity: 'geom,city,cartodb_id,date_proc'
};

var linesLayerDef = {
  sql: "select * from carmen_amaya_concerts",
  cartocss: "#carmen_amaya_concerts{ line-width: 1; line-color: #FFF; line-opacity: 0.8; }",
};

cartodb.createLayer(map, {
  user_name: 'carmenamaya',
  type: 'cartodb',
  sublayers: [
//        baseLayerDef,
    pointsLayerDef,
    linesLayerDef
  ]
}).on('done', function(layer) {
  //     layer created, add it to the map
      map.addLayer(layer);


texts: How can I clip multiple geotiff tiles with a shapefile without greatly increasing their size?

texts:
I'm using gdalwarp to clip a number of geotiffs (tiles of a larger map) to a polygon described by a shapefile (which is about the size of the overall map, i.e. much larger than individual tiles) using:

gdalwarp -of GTiff -cutline shape.shp -crop_to_cutline tile.tif cropped_tile.tif

(where tile.tif actually iterates through all the map tiles in the full code)

This works fine, however the resulting tiles (cropped_tile*.tif) now cover the same geographic area as the shapefile, and consequently have much larger file sizes (2.5M -> 64M).  Is there any way to clip the areas of the tiles which lie outside the polygon in the shapefile without making each tile cover the same area as the shapefile itself?


texts: Keep the same image scale after mosaicking

texts:
I am trying to keep the same original image scale after creating a mosaic with the same images. But when I try to add in a map this image I can´t configure the visualized parameters.

The problem is the image is never shown, it seems to be all masked.

Is there any way to keep the image scale and can be possible to see the image after?

I used this function:

var composite_free_mosaic_nova = composite_free_mosaic.resample('bilinear').reproject({
  crs: projecao.projection().crs(),
  scale: 30
});

After the mosaicking the nominal scale was: 111319.49079327357. And I try to keep in 30 meters.

My code link is this one:

https://code.earthengine.google.com/bd6aec5ffd9f98cb2d0284e5c1ece421


texts: Add a cursor/icon in OpenLayers to mark precisely

texts:
We are using OpenLayers in a mobile application developed in Ionic. When we trace any feature (Point/Line/Polygon) by using finger precision is not accurate and user is not able to exactly mark the specific point.

Is there any way that we can add icon (plus/cross) and adjust the point to the center of the icon, Now if user tap finger within a tolerance range then the point at the center of the icon is considered.


texts:
I am looking for a quick way of defining several variables at once in GEE. My code currently is written like this:
var a = 0;
var b = 0;
var c = 0;
var d = 0;

which takes up space uselessly and is quite repetitive.
I found that in javascript it is possible to write it as:
var [a,b,c,d] = [0,1,2,3]

but trying that in GEE only produces an error:

(SyntaxError: Unexpected token (1:5)).


texts: Collecting geometric network distances and updating field with these distances for whole feature class?

texts:
I have a geometric network representing a golf course irrigation system. The pump house is the only source and each irrigation head is a sink. I want to update a field in the heads fc ([Dist_from_Pump]) based on each head's distance from the pumphouse. I suspect this will be a Python solution, which I'm fine with, but how can I iterate over all of the heads, place a flag at each head, and then trace the network for a distance?

Or is there a different/better way?

I'm thinking something like:

    for head in heads:
      dist = [trace distance from head to pump]
      # then use an update cursor to add these distances to the attributes table. But I don't know how to get these distances from arcpy.

I'm using ArcMap 10.5


texts: Resolving relation does not exist in PostGIS?

texts:
I have a view, named PointsForGpxExport, that contains points I want to export to a GPX file. The ogr2ogr command I'm trying to use for the export is included below. When I run the command, a GPX file is created. The file is a shell of a GPX file - headers and root element but none of the data from the view. An error is reported.

The error reported is

ERROR 1: ERROR: relation "pointsforgpxexport" does not exist
LINE 1: DECLARE executeSQLCursor CURSOR for SELECT * FROM PointsForG...

The command:

ogr2ogr-f GPX c:\temp\points.gpx PG:"host=localhost port=5432
dbname=SpatialPlayground schemas=public user=postgres password=password"
-sql "SELECT * FROM PointsForGpxExport"

When I run SELECT * FROM PointsForGpxExport inside the pgAdmin tool, an error occurs ERROR: relation "pointsforgpxexport" does not exist I resolved this error by including the schema name in the query as in SELECT * FROM "public"."PointsForGpxExport"; This error is identical to the error I'm getting when I run ogr2ogr and I thought it would be resolved by using schemas=public in the connection string.

What am I doing wrong and how do I fix this?

PostgreSQL 9.3, PostGIS 2.1.6, GDAL 1.11.2 released 2015/02/10


texts: Create cross sections in QGIS and save as DXF file

texts: Is there any way to create cross sections of rasters in QGIS and save it as DXF file? I need to continue working with it in AutoCAD LT.


texts: Convert invalid multilinestring to linestring with recursive ST_Union PostGIS

texts:
intro

I explain my geoprocess:

start of a very long linestring
Segmented linestring, 20 meter segments with ST_Segmentize, save the order of segments, ex 0 is start segment and 200 end segment
I need to select certain segments within the linestring; selecting by order, ex: select segments between 10 and 20
with the selected segements I need to join them in 1 single
linestring

problem

the union of the selected segments does not generate a valid segment according to the answer

Converting MultiLinestring to Linestring with PostGIS?

I usually use st_union to join the segments, and it works correctly while these are valid; the problem is that I have linestring not valid; this generate multilinestring

this is the case:

posible solution???

I think the best solution could be a recursive st_union that links all the segments in a certain order

generating step by step, singing linestring to generate a valid and simple linestring, this approach would be correct?

or can there be any tool in postgis that realizes this?


texts: arcpy.da.NumPyArrayToTable() RuntimeError: Unspecified error

texts:
Using arcpy for ArcGIS Desktop 10.7.1...

In order to append data from a Pandas data frame to an ArcGIS table, I am attempting to create a temporary in-memory table from the Pandas data frame as follows:

tmpTable = "in_memory/tmpTable"
labResults = pandas.read_csv(csvPath, parse_dates=[3], dayfirst=True)

######  If this ONE next line is commented out and the following line is UNCOMMENTED, the error does NOT occur
labResults['ResultQualifier'] = np.where(labResults.Result.str[0].str.isdigit(), None, labResults.Result.str[0])
#labResults['ResultQualifier'] = np.where(labResults.Result.str[0].str.isdigit(), labResults.Result.str[0], labResults.Result.str[0])

labResults['ResultValue'] = pandas.to_numeric(np.where(labResults.Result.str[0].str.isdigit(), labResults.Result, np.where(labResults.Result.str[1].str.isdigit(),
labResults.Result.str[1:], labResults.Result.str[2:])))
del labResults['Result']
labResults.rename(columns={'ResultValue': 'Result'}, inplace=True)

arr = np.array(np.rec.fromrecords(labResults.values))
arr.dtype.names = tuple(labResults.dtypes.index.tolist())
arcpy.da.NumPyArrayToTable(arr, tmpTable)
######  ERROR OCCURS on this next line
arcpy.Append_management(tmpTable, resultsTable, "NO_TEST")

But the arcpy.da.NumPyArrayToTable() fails with the frustratingly generic error:

`RuntimeError: Unspecified error`

I can print the data frame by entering labResults on the Python console, and it prints out the expected data.  And labResults.dtypes produces the expected output also, so I'm confident that the data frame is all OK.

NB:  I get the same error whether tmpTable is a fGDB table (fGDB exists, table does not) or an in_memory table.

How can I get this to work as expected?

Note:  This code is based on:  https://my.usgs.gov/confluence/display/cdi/pandas.DataFrame+to+ArcGIS+Table where there is a full explanation of what it all does.

SAMPLE DATA:

"JobNo","LabNo","SampleID","DateSampled","Method","AnalyteCode","Units","Result"
90529,366330,"site1",15/1/19 10:15:00,"Glyphosate and Glufosinate in Water by LC-MS/MS","Glyphosate","µg/L","<5"
90529,366330,"site1",15/1/19 10:15:00,"Solids, Total Suspended dried at 104 ± 3 °C","TSS","mg/L","21"
90529,366331,"site2",15/1/19 10:35:00,"Glyphosate and Glufosinate in Water by LC-MS/MS","Glyphosate","µg/L","<5"
90529,366331,"site2",15/1/19 10:35:00,"Solids, Total Suspended dried at 104 ± 3 °C","TSS","mg/L","577"
90529,366332,"site3",15/1/19 11:00:00,"Glyphosate and Glufosinate in Water by LC-MS/MS","Glyphosate","µg/L","<5"

UPDATE:

I have found that if I change the labResults['ResultQualifier'] = line to never include NULLs (None's), then the error does NOT occur.  This line is there to create a new column in the data frame that includes either < if the first character of 'Result' is not a digit, or otherwise 'None' (NULL).  I'm not sure why NULL values would be causing such a problem.


texts: Leaflet: Linking a marker click event to an image

texts:
I'm breaking down my question into three parts.  I am a beginner to Leaflet and to writing code.  Currently, I have a JSON file with all of my markers linked to my map page. The idea is that every marker will have a different image so when the user clicks on the marker the thumbnail for the image will appear in the map, which will link to an image in the website.

Here is an example of a marker in the JSON file:

var mapshs4 = {
    "type": "FeatureCollection",

    "features": [

    { "type": "Feature", "id": 0, "properties": { "idcode": "A", "popupContent": "Here!" + '< img src="thumbnail">' }, "geometry": { "type": "Point", "coordinates": [ -89.2, 22.4 ] } }

Is there a way to resize the popup box to the image in the popup of a marker.
At the same time, I would like the click event to connect to a larger sized image embedded within the site.  My guess is that I would have to attach an ID to each image and to each marker to link them together.
Conversely, how do you make it so that when the user clicks on the image embedded in the site, it will go to that marker on the map?  (For an example that uses a label: http://www.theglobeandmail.com/life/home-and-garden/real-estate/exclusive-how-much-are-homes-in-your-toronto-neighbourhood/article4520994/)

Any help would be much appreciated.


texts: Select by attribute expression invalid in python

texts:
I'm not sure how my expression is failing. Any ideas?

inPoints = "points.shp"
i = 0
arcpy.MakeFeatureLayer_management(inPoints, "pts")
arcpy.SelectLayerByAttribute_management("pts", "NEW_SELECTION", '\"FID\" = i')


texts: Limiting number of features loaded by QGIS2Leaf?

texts:
I'm using the QGIS2Leaf plugin. As others have noted, Leaflet slows down when the number of features becomes too large.

https://stackoverflow.com/questions/15567378/efficient-way-of-accessing-large-data-sets-to-display-on-leaflet

I've worked out how to limit zoom levels, and have the map load zoomed right in. I'm now trying to work out how to just load the features within those bounds, and have the load refresh as the user pans to change those bounds.

Could anyone familiar with Leaflet/QGIS2Leaf tell me if this is possible, and further how to achieve it?


texts: Converting from QML to SLD on the command line

texts:
We have the problem of automating the conversion of MapInfo TAB styles to SLD. I've managed to achieve the first part of the process (TAB -> QML) using this tool:

https://github.com/NathanW2/MapInfo-to-QGIS-style-generator

I can load the resulting QML file into QGIS and export it as SLD manually. However, would I really need is a way of doing this via the command line.

Is this possible with QGIS, or does anyone know another tools that might help here?


texts: Python loop stops before running all files

texts:
I have a Python loop, the loop has to run commands on 200 plus files. But it randomly stops and does not continue producing the new output files.

here is the code:

from subprocess import Popen, PIPE
import os,glob

os.chdir(r'C:\Program Files\QGIS 3.0\bin')
grassbat = r'C:\Program Files\QGIS 3.0\bin\grass74.bat'  #not needed as the `os dir is changed`,

osgeobat = str(r'C:\Program Files\QGIS 3.0\OSGeo4W.bat')
input_las = r'E:\Rohan\Heimdall\data\LAZ\heimdall_final_demo_lidar'
output_directory = r'E:\Rohan\Heimdall\data\Raster\GRASS'
grass_location = r'E:\Rohan\Heimdall\grassdata\lidar\PERMANENT' # r'  reads a path as a raw string
print(grass_location)

# Loop to transform LAS to raster tif.
def las_to_GTIFF_lidardata(input_las, grass_location, output_directory):
    las_file_list = glob.glob(input_las + '\\*.las')  # lists all files in folder that (*.las)
    print(las_file_list)
    for file in las_file_list:
        tif_file_name = os.path.basename(file)[:-4] + '.tif'
        out_tif = os.path.join(output_directory, tif_file_name)

        firstcommand = str(osgeobat + " grass74 -text -f E:\Rohan\Heimdall\grassdata\lidar\PERMANENT --exec r.in.lidar -e -n -o --overwrite --verbose input=" + file +" output=lidar_in resolution=2")
        secondcommand = str(osgeobat + " grass74 -text -f E:\Rohan\Heimdall\grassdata\lidar\PERMANENT --exec r.out.gdal --overwrite --verbose -m input=lidar_in@PERMANENT output=" + out_tif + " format=GTiff type=Float32 nodata=-9999")
    firstcommandlst = firstcommand.split(" ")
    secondcommandlst = secondcommand.split(" ")

    las2raster_commands = []
    las2raster_commands.append(firstcommandlst)
    las2raster_commands.append(secondcommandlst)

    for command in las2raster_commands:
        print("command: ")
        print(command)
        process = Popen(command, stdout=PIPE, stderr=PIPE)
        stdout = process.communicate()
        print(stdout)
return True

las_to_GTIFF_lidardata(input_las=input_las, grass_location=grass_location, output_directory=output_directory)


texts: Calculating median elevation for local authority areas using OS Terrain 50

texts:
I haven't been able to find anything that could comprehensively help me with my problem.

I want to calculate median elevation for UK local authority areas. So far I've taken the OS Terrain 50 data  which is available in three formats [ESRI ® Shape (Contours): 1.1 Gb], [GML (Contours): 1.1 Gb], [ASCII Grid and GML (Grid): 162 Mb].

I think I need to intersect the height data with the local authority boundaries, but I'm struggling. There are lots of related questions on OS Terrain 50 such as here, here and here but some of them focus on loading in the data in different formats on other operating systems. I'm on Windows 7 and attempting to do this in QGIS. I have limited Python experience, but I suspect I'm going to have rectify that to solve this problem.

How do I load in the very large number of WinZipped files associated with the OS Terrain 50 and then take the intersection with the local authority boundaries, in order for me to calculate basic statistics on the associated elevation attributes?


texts:
I've tried the solution in this post but haven't been able to resolve so hoping that someone here might be able to help.
I'm attempting to create a buffer layer of 5km radius around a collection of pinpoints that I have in a separate layer.
The layer that has the pinpoints I want to buffer is in CRS EPSG:4326 - WGS 84, but my overall project map is EPSG:3857 WGS 84 / Pseudo Mercator.
If I change the layer with my pin points to 3857, all the pinpoints relocate to one pinpoint on a completely different continent. If I change the overall project to 4326, my base map (Open Street Map) distorts and stretches horribly.
Based on previous answers on older versions of QGIS I've also looked for On-The-Fly projections settings but they don't appear to be here - I've tried both Project > Properties > CRS and also Settings > Options > CRS as suggested in a duplicate thread, but it's not in either spot.
How can I get circular buffers from my pinpoints that doesn't horribly distort my base map, or make my pinpoints disappear into the sea?
I've found a workaround which is to make a radius circle using symbology, but this obviously isn't ideal because it means I can't use the added benefits of a buffer, like being able to count the number of pinpoints from another layer that fall within the buffer.
QGIS 3.12.3 Bucaresti Mac 10.15.5 Catalina


texts: Adding a second field to QgsGraduatedSymbolRenderer using PyQGIS

texts:
I am using PyQGIS within my work and I currently have a script that uses QgsGraduatedSymbolRenderer to symbolise my data. Previously, I have used a script that symbolises my data using CategorizedSymbolRenderer (see script below). However, this prints out too many values on my exported map, so I have decided to use Graduated instead, so that I can use min/max values to classify my data.
default_style = QgsStyle().defaultStyle()
color_ramp = default_style.colorRamp('Spectral')  # Spectral color ramp
color_ramp.invert()

layer = QgsProject().instance().mapLayersByName('Intensity')[0]

exp_string = '\'Intensity: \' || "Intensity" || \'\nAffected Population: \' || sum("LAC_Popu_1", "Intensity")'
unique_values = []
exp = QgsExpression(exp_string)
#print(exp.isValid())
context = QgsExpressionContext()
context.appendScope(QgsExpressionContextUtils.globalScope())
context.appendScope(QgsExpressionContextUtils.projectScope(QgsProject.instance()))
context.appendScope(QgsExpressionContextUtils.atlasScope(None))
context.appendScope(QgsExpressionContextUtils.layerScope(layer))
exp.prepare(context)
for feat in layer.getFeatures():
    context.setFeature(feat)
    val = exp.evaluate(context)
    if val not in unique_values:
        unique_values.append(val)
# print(unique_values)

categories = []
for value in sorted(unique_values):
    symbol = QgsSymbol.defaultSymbol(layer.geometryType())
    category = QgsRendererCategory(value, symbol, str(value))
    categories.append(category)
renderer = QgsCategorizedSymbolRenderer(exp_string, categories)
renderer.updateColorRamp(color_ramp)
layer.setRenderer(renderer)
layer.setOpacity(0.55)
layer.triggerRepaint()

My script that I have created so far retrieves intensity data from a field called 'Intensity' within my vlayer and uses Graduated Symbology to say if the intensity is between 1.0 - 1.99 for example, then a blue colour is given. This script runs perfectly fine to produce graduated intensity values with corresponding colour. However, I am struggling to get the 'Affected Population' factor within my label, like previously shown, where a value is underneath each intensity. The affected population field will have no involvement as to how the graduated symbology is laid out, it is just extra information for each intensity. This script is as follows:
Intensity_rules = (
    ('MMI Intensity 1.0 - 1.99', 0.0, 1.99, '#2b83ba'),
    ('MMI Intensity 2.0 - 2.99', 2.0, 2.99, '#71b4ae'),
    ('MMI Intensity 3.0 - 3.99', 3.0, 3.99, '#abdda4'),
    ('MMI Intensity 4.0 - 4.99', 4.0, 4.99, '#d4eeb1'),
    ('MMI Intensity 5.0 - 5.99', 5.0, 5.99, '#ffffbf'),
    ('MMI Intensity 6.0 - 6.99', 6.0, 6.99, '#fed58e'),
    ('MMI Intensity 7.0 - 7.99', 7.0, 7.99, '#fdae61'),
    ('MMI Intensity 8.0 - 8.99', 8.0, 8.99, '#ee7144'),
    ('MMI Intensity 9.0 - 9.99', 9.0, 9.99, '#d7191c'),
    ('MMI Intensity 10.0 - 10.99', 10.0, 10.99, '#d7191c'),
)

# create a category for each item in intensity
ranges = []
for label, lower, upper, color in Intensity_rules:
    symbol = QgsSymbol.defaultSymbol(vlayer.geometryType())
    symbol.setColor(QColor(color))
    rng = QgsRendererRange(lower, upper, symbol, label)
    ranges.append(rng)

#create the renderer and assign it to layer
expression = '"Intensity" ' #field name '\'Intensity: \' || "Intensity" || \'\nAffected Population: \' || sum("LAC_Popu_1", "Intensity")'
renderer = QgsGraduatedSymbolRenderer(expression, ranges)
vlayer.setRenderer(renderer)
vlayer.setOpacity(0.70)
vlayer.triggerRepaint()

I run zonal statistics on the intensity shapefile layer, along with a population raster, to produce an affected population column called "LAC_Popu_1". Is there a way I can add this column of produced statistics underneath each line of intensity like in this first script/photo where I used Categorized symbology e.g.
MMI Intensity 1.0 - 1.99: 'Blue colour'
Affected Population: 545432


texts: How to get DXF to show up in the right spot?

texts: I just followed the advice of this forum post in order to export the vector layers of a PDF for usage in ArcMap. The map vector layers in the PDF were spatially referenced. I was able to export the layers as a DXF using Inkscape and open the layers in ArcMap. I then went to the properties of the DXF in catalog and set the proper projection. However, when I brought the DXF into the map document, it was not showing up in the right location. It appears the the spatial information was somehow removed in the process of exporting into a DXF. Does anyone know of a solution to my situation? Thanks!


texts: Adjust image size in popup to fit inside popup container

texts:
I would like to display images in popups. It works but the problem is the images are too large:

I've added a class in the img tag in my js-file:

var picURL2 = 'http://www.lib.utexas.edu/maps/historical/newark_nj_1922.jpg';

shelter1.bindPopup("<img src='" + picURL2 + "'" + "class=popupImage" + "/>");

The class in my css file looks like this:

.popupImage {
    max-width: 20%;
    max-height: 20%;
}

But the image still streches beyond the popup container. The console doesn't produce any errors so I assume I've added the class inside the img tag correctly.


texts: Measuring elevation at points on DTM surface using R

texts:
Based on answers to Preventing reflective surface in rasterVis? I have been able to create natural looking digital terrain models using a combination of raster, rasterVis and rgl working on data held in ESRI ASCII grid files.

Next I want to measure the the height (elevation) of points on the surface of the digital terrain model. Ultimately I am interested in establishing the height of certain points on the hills compared to the floor of the valley nearby. Obviously, I could look at the data in the ASCII file - this is raw elevation data, after all - but it is impossible to relate the points in the plot to the data in the file. I would like to be able to visually inspect the plot, click on a point and read off the elevation, then click on another point and read that elevation.

While plotting the data in R has been very easy, I can't think of any obvious way to measure heights using an R-based solution.

I seek suggestions on how to approach this using R. Image showing what I'm after and the code used to produce it follow below. (NB I do not have permission to reproduce the original DTM data, hence the mathematical surface.)

# Create an ESRI ASCII raster grid and plot it.

out.file <- "your_file_name_here.asc"

num.cols <- 50
num.rows <- 50

# Mathematical surface taken from ?persp example
x <- seq(-10, 10, length.out = num.cols)
y <- x
rotsinc <- function(x,y) {
    sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }
    10 * sinc( sqrt(x^2+y^2) )
}

my.mat <- outer(x, y, rotsinc)

# Create header as per format specs
header.df <- data.frame(variable = c('ncols','nrows','xllcorner','yllcorner','cellsize','nodata_value'),
                        value = c(num.cols,num.rows,378923,4072345,2,-9999))

fileConn <- file(out.file)
write.table(header.df,
            quote = FALSE,
            row.names = FALSE,
            col.names = FALSE,
            sep = " ",
            file = fileConn)

# Append matrix to file
fileConn <- file(out.file, open = "a")
write.table(my.mat,
            row.names = FALSE,
            col.names = FALSE,
            file = fileConn)

close(fileConn)

# Now we can plot
library(raster)
library(rgdal)
library(rasterVis)

my.raster <- raster(readGDAL(out.file))
plot3D(my.raster,
       smooth = FALSE,
       specular = "black")


texts: Looking for dissolve algorithm for Javascript

texts: I have three adjacent polygons that together make up a larger shape. I need an algorithm that will dissolve the internal borders leaving only the outside border. An example would be a group of adjacent counties. I want to merge together all the county polygons and dissolve the internal borders so I'm left with just the outline of the state.


texts: Visualization libraries for rendering 3D GIS world on iOS/Android

texts:
Are there any good 3D libraries/frameworks for rendering 3D worlds in geospatial coordinates on mobile platforms?

Use case: I have some known 3D geospatial features (eg. buildings) with their corresponding lat,long,elevation values as well as an underlying DEM. I would like to render such a "scene" from a given camera viewpoint on either an iOS or an Android device.

I could use opengl-es but it would take some effort to convert the geospatial coordinates to cartesian coordinates. Is there a "geospatial-aware" soft of 3D library for mobile platforms that I could use? Something like the AGI insight3D libraries but for mobile.

Thanks!


texts: How Do I Capture when the Cluster Rendering is Completed

texts:
I have a graphics layer that draws clusters using the sample custom cluster from ESRI.  http://help.arcgis.com/en/webapi/silverlight/samples/start.htm#CustomClusterer

We have a Busy Indicator that we pop up during long calculations to let the user know that something is happening.  I would like to pop up this indicator during the clustering, and have it go away once the clusters have done drawing on the screen.  I can turn it on, but I am not able to capture a "done rendering" event to turn off the indicator.  I've tried capturing the  LayoutUpdated event on the Map, but it gets sent several times and am not able to determine which one is the graphics layers saying its done rendering.  I found the  OnClusteringCompleted event of the GraphicsClusterer and tried to override it by adding it to the custom clusterer but the code never gets there.

Any suggestions on how to accomplish what I want?


texts: Conversion of 3 degree UTM zone Easting Northing to WGS-84

texts:
In Nepal we use 3 degrees width zone UTM with 3 central meridians 81 E, 84 E and 87 E and datum: Everest-1830. I am trying to convert the Everest-projected coordinates to WGS-84's Latitude and Longitude. I applied reverse of projection for conversion of Everest projected Easting and Northing to Everest's Latitude and Longitude and then transformed to WGS-84's Latitude and Longitude

I used the formula in ogp publications 373-7-2 geomatics guidance note number 7 page number 50 and 54 to calculate latitude and longitude from Easting and Northing Projected coordinates and created code to do so.

But whenever I tested the conversion and input Easting, Northing and height values 617000, 3038250, 1230 respectively of zone 88 (Central meridian 84 degree E, datum : Everest-1830).
 I get output latitude 27.461834699447184 and longitude 83.86573999170062 (Datum: Wgs-84). But the expected output is latitude: 27.732 and longitude: 85.417.

I must be doing something wrong. Please help me


texts: I am trying to use Scanaerial in JOSM to locate water sources in my county.  I was able to install all of the required Python modules and ExtTools and load the scanaerial plugin, but now when I run it, it just scans the WMS layer endlessly and never gives me any output, at least for the hour I have let it run on an acre sized area..  All the documentation is on the OSM Wiki and the files are at this GitHub repo.  I am not sure what I am missing to make this plugin produce some output.


texts: Best Landsat-5 TM band combination for detecting fire scars

texts: I am looking for the best band combination with Landsat 5 imagery to detect fire scars.


texts: Convert Geo Coordinates System X/Y to lat &lon from netCDF4 file using python 3

texts:
I have netcdf4 file with some variables like: Geo Coordinates System X/Y.
I need to convert it to lat/lon. How can i do that with python?

This information from my netCDF4 file:

    import numpy as np
    from netCDF4 import Dataset

     path = '/home/data/Mmultic3kmNC4_msg04_201905080200.
        ...: nc'

    In [78]: f=netCDF4.Dataset(path,'r')

     f.variables.keys()
    Out[80]: odict_keys(['time', 'dtime', 'IR_016', 'commentaires', 'satellite', 'geos', 'ImageNavigation', 'GeosCoordinateSystem', 'Y', 'X', 'Albedo_to_Native_count_IR_016', 'IR_039', 'Temp_to_Native_count_IR_039', 'IR_087', 'Temp_to_Native_count_IR_087', 'IR_097', 'Temp_to_Native_count_IR_097', 'IR_108', 'Temp_to_Native_count_IR_108', 'IR_120', 'Temp_to_Native_count_IR_120', 'IR_134', 'Temp_to_Native_count_IR_134', 'VIS006', 'Albedo_to_Native_count_VIS006', 'VIS008', 'Albedo_to_Native_count_VIS008', 'WV_062', 'Temp_to_Native_count_WV_062', 'WV_073', 'Temp_to_Native_count_WV_073'])

    In [81]: IR=f.variables['IR_108']

    In [82]: IR
    Out[82]:
    <class 'netCDF4._netCDF4.Variable'>
    int16 IR_108(ny, nx)
        standard_name: toa_brightness_temperature
        long_name: 10.8 microns infra-red channel MSG Imager
        units: K
        scale_factor: 0.01
        add_offset: 273.15
        valid_min: -27315
        valid_max: 32767
        _FillValue: -32768
        available: yes
        comment: channel 9
        nuc: 931.122
        alpha: 0.9983
        beta: 0.6256
        CalibrationParametersOrigin: rttov
        _CoordinateSystems: GeosCoordinateSystem
        grid_mapping: GeosCoordinateSystem
        resolution: 3.0
        Offset : radiance = offset + CN x slope: -10.456819486590666
        Slope : radiance = offset + CN x slope: 0.2050356762076601
    unlimited dimensions:
    current shape = (3712, 3712)
    filling on

    print(IR.dimensions)
    ('ny', 'nx')
x=f.variables['X']

In [3]: x
Out[3]:
<class 'netCDF4._netCDF4.Variable'>
float32 X(nx)
    standard_name: projection_x_coordinate
    units: m
    _CoordinateAxisType: GeoX
unlimited dimensions:
current shape = (3712,)
filling on, default _FillValue of 9.969209968386869e+36 used

In [4]: y = f.variables['Y']

In [5]: y
Out[5]:
<class 'netCDF4._netCDF4.Variable'>
float32 Y(ny)
    standard_name: projection_y_coordinate
    units: m
    _CoordinateAxisType: GeoY
unlimited dimensions:
current shape = (3712,)
filling on, default _FillValue of 9.969209968386869e+36 used

y = f.variables['Y'][:]

In [12]: y
Out[12]:
masked_array(data=[ 5568754. ,  5565753.5,  5562753. , ..., -5559743. ,
                   -5562743. , -5565744. ],
             mask=False,
       fill_value=1e+20,
            dtype=float32)

link for file netCDF4

What do I need to use to convert X/Y to lat /lon coordinates?


texts: Resuming (continuing) digitizing polygon feature from any vertex in ArcMap?

texts:
When I digitize in ArcMap with a Wacom tablet and a pen I often accidentally double click and finish my sketch.  More often then not the accidentally finished sketch crosses itself.  I would like to be able to continue digitizing the existing sketch from the last vertex or any vertex I choose.  The OOTB "Continue Feature" tool does not work in this situation because it just adds a new feature part to the sketch and does not resolve the "crossing segment" artifact.

I also tried using the ESRI's "Modify by Continuing Sketch" tool (developer sample) which does what I need it to do but it always resumes digitizing from the "Last" (red) segment.  In theory this is correct unfortunately last vertex is often not the last digitized vertex - Editor seems to assign the last (red) vertex randomly.

My question is: As a workaround, I would like to be able to assign the "last vertex" to the actual last or any vertex in the sketch I choose.  Can this be done?

Figure: Example of an accidentally finished sketch (this is a polygon with a transparent fill) in which the last vertex was assigned to polygon boundary crosses itself not where I finished digitizing (which is indicated by an arrow).  This also happens with simple polygons - last vertex becomes another random vertex, not the actual last vertex I created in the sketch.


texts: Refreshing IDockableWindow with DGV using C#

texts:
I am having a IDockableWindow which will display a row of data which is a DATATABLE on DGV.
below is my code.

Page1.cs:

 private void toolStripMenuItem_Open_detail_page_Click(object sender, EventArgs e)
    {
        UID dockWinID = new UIDClass();
        dockWinID.Value = ThisAddIn.IDs.DetailsPage;

        IDockableWindow dockWindow = ArcMap.DockableWindowManager.GetDockableWindow(dockWinID);
        dockWindow.Show(true);

        DNRToolBar.WaterRightsDetailsPage WRDetailsPage = new DetailsPage(this.Hook);
        WRDetailsPage.WRDetailsPageLOAD(3);//passing in an integer directly for test.
    }

Page2.cs: - this is the default page which the arcGIS template creates once you choose to create a dock able window and below is my method in it

  public void WRDetailsPageLOAD(int ID)
    {
        WaterRightsInfoDBAccess dbAccess = new WaterRightsInfoDBAccess();
        DGV_WaterRightsDetailsPage.DataSource = dbAccess.WRDetailsPageData(ID);
    }

Once the Dockable window is shown with the DataTable when I try to call the click event again it does execute the WRDetailsPageLOAD(int ID) but the DGV in the dockwindow does not refresh.

http://forums.esri.com/Thread.asp?c=93&f=993&t=110408 this is exactly the same issue but I dint find it helpful


texts: Buffering only sides of polygons in QGIS?

texts:
My goal is to create a buffer from a polygon. But the edges are the problem. I can't really explain it in words, so I will post a screenshot. That's the way I want it to look. So the goal is, that the edges are not connected. The example on the screenshot is done with powerpoint, but I want to do it with QGIS, so I can do it for a lot of polygons at once.

Does anyone have an idea?


texts: IF/ELIF statement returning NULL for values that are not calculated with field calculator in Python

texts:
I am using field calculator to change several values of field.  I am receiving NULL values for those which are not in my function. What is the cause of this and a workaround?

In the code below; I would like to change only the values in my function.  If a value is not in the function; the code automatically sets the new value as NULL.

import arcpy

codeblock = """def replacestring(DAMAGING_PARTY):
        output = None
        if DAMAGING_PARTY == "Arizona Pipeline contact is Chad":
            output = "Arizona Pipeline Company"
        elif DAMAGING_PARTY == "James A Shirley Construction":
             output = "James A. Shirley Construction, Inc."
        elif DAMAGING_PARTY == "NPL Construction (Kevin Sena)":
            output = "NPL Construction"
        elif DAMAGING_PARTY == "NPL Construction DENNIS PARKHOUSE":
            output = "NPL Construction"
        elif DAMAGING_PARTY == "Rami / Benito Sustaita":
             output = "Rami Alabassi"
        elif DAMAGING_PARTY == "Shea Homes (Contractor Working for)":
             output = "Shea Homes"
        elif DAMAGING_PARTY == "Shea Homes/Rudy Andrade":
             output = "Shea Homes"
        elif DAMAGING_PARTY == "W A RASIC CONSTRUCTION/CODY ROLLER":
             output = "W A RASIC CONSTRUCTION"
        elif DAMAGING_PARTY == "WARASICK CONSTR COMPANY/GIANCARLO PAMDOSI":
             output = "W A RASIC CONSTRUCTION"
        elif DAMAGING_PARTY == "Young Construction; contact Rick Musmecci":
             output = "Young Construction"
        return output"""

ShortCodeExpression = "replacestring(!DAMAGING_PARTY!)"

table = r'table'

arcpy.CalculateField_management(table, "DAMAGING_PARTY", ShortCodeExpression, "PYTHON_9.3", codeblock)


texts: OpenStreetMap displaying very different data depending on zoom level following edits

texts:
Yesterday (19 hours ago) I made an edit to OpenStreetMap to reflect walkways around a new local train station. Some informal walkways have been blocked by a wall and some new walkways have been constructed around the station, bus stop and car park.

These edits are showing up on the map at openstreetmap.org at the zoom level map=18 and higher, but at map=17 and lower the old walkways are visible and my edits are not showing up. The images below illustrate:

I made the edits in JOSM for Windows. I did not notice any way of making different maps based on the zoom level. The OSM help suggests that the main Mapnik map should be updated every few minutes or in some cases a few hours, which leads me to think I may have done something wrong in my edits and would like to confirm this as I would like to be a good custodian of OpenStreetMap in my future contributions.

Have I done something wrong here?
How long does it take for Mapnik to update all the zoom level tiles?


texts: Superimposing high resolution NetCDF data on Google Earth?

texts: I have got output from Weather Research Forecasting Software(WRF) in NETCDF format. These are basically high resolution forecasts for precipitation and wind barbs(sea breeze circulations or slope valley circulations). I want to be able to display this on Google Earth or some similar GIS software. It should be interactive meaning I should be click to the highest possible resolution to tell the viewer in which neighborhood the rainfall occurred. It should also show the terrain etc if in mountainous area. The tools on this web page - http://www.nusculus.com/wtools allows me to see static geographical domain information on Google Earth by converting WRF output(in NETCDF format) into KML format. I want to be able to project meteorological data onto Google Earth. Does anyone know of any tool(could be in any language-java,python) out there that does this ?


texts:
I'm making a choropleth map based on this example

I would like to add links to the map, so that when a user clicks on a state they will be taken to a specific page.

For example, a user clicks on California and is taken to a page about California.

I saw a great example of this, but the site has since been taken down so I can't link to it.


texts: Is there a way to create a DEM that follow the point cloud so you don't need to clip the DEM after?

texts:
I would like to create a DEM from an XYZ file (gridded 0,25x0,25) with 0,25 cell. I'm wondering if there is a way so that my interpolation follows my point cloud, instead of covering a rectangle like in the below picture. I've tried the output option in QGIS called "Use layer extent", thinking it would use the boundaries of the point cloud and clip it automatically.

So I'm wondering if there is a way so you don't need to mask/clip the raster after the interpolation, so that it follows the edge of the point cloud.


texts: Help with Rastercalculator - Filter areas with a certain height difference

texts:
I'm struggeling a little bit with the rastercalculator tool of QGIS.
At the moment I am working with a DSM and a DTM. My aim is to filter all wind turbines out of the DSM and replace it with the corresponding DTM data.
Here im thinking about a calculation like:

If the height difference in "DSM" is like 50 m above the average in a certain area, replace it with pixels from "DTM"

So my question is, is this the right approach, or is there an easier way to deal with this problem?


texts: Use a SQL 2005 table as attribute data in QGIS

texts:
We have financial data in a SQL 2005 database which we would like to link to a cadastral shapefile in QGIS. The idea is to be able to display up-to-date outstanding balances for each plot.

How can we link QGIS to the SQL table?

We are using QGIS 2.0.


texts:
I have a line that exists between two latlon positions, and a point on a certain latlon position. What I want to know is what the nearest point on that line is (in respect to the other point). So a latlon position on that line.

I know how to do the basic math, but I just can't get my head around the latlong calculation. Can anyone help me or give me some pointers on how to do it?

Thanks


texts: Identifying value of closest non-NA pixel

texts:
I have a raster representing a large landscape, where pixel values indicate particular land classes.  I have been using the "distance" function in Package 'raster'  in R to calculate the distance from all NA pixels to the nearest pixel with a value.  In addition to the distance, I'm interested in knowing what the value was associated with the distance calculation.  For example, if I have two land class values, e.g. 1= "Urban" and 2= "Agriculture", I'd like to know if a pixel is closest to an "Urban" pixel or an "Agriculture" pixel.  I looked into the source code for the distance function, which converts the pixels to points and then calculates all relative distances but only returns the shortest distance.  See below for an example of the distance function:

    library(raster)
    r <- raster(ncol=36,nrow=18)
    r[] <- NA
    r[500] <- 1
    r[100] <-2
    dist <- distance(r) #calculates the distance to the nearest non-NA pixel

Is there a function in a different package that does this (I'd like to stay in the R environment)?


texts:
I am finished with my QGIS mapping project and now I want to save it and be able to view it from desktop only.

I am not wanting to publish this map to the internet and I do not want the other users on these particular computer to be able to modify/edit this map only view.

So the question is how would one go about doing this. Do I need a special "Map Viewer" to be able to do this?


texts: How to create a multi-story building infrastructure web map?

texts:
I have absolutely no experience working with GIS before.  During a recent meeting, our client was bouncing around their complaints about how detailed our annual reliability reports were and openly admitted to just reading the covering letter and the introduction before skipping 500-odd pages to the recommendation section.

They are looking for a means to visually represent the mechanical and electrical health of plant machinery which is not a major issue as we build a web portal with WPF dashboards, guages etc. which they seemed impressed with.  However, because each building is multi-story they wanted a 'Google Maps' type interface control to navigate between building levels and ultimately zoom from floor --> room --> plant machine --> machine component.  Again this is achievable through WPF controls but to get it to work we have to build each individual map from scratch.

After searching the net...  I ended up here.  Most of the application resources I have stumbled across discuss working with existing SHP files and I am searching for resource that focuses on creating SHP files from scratch.


texts: Importing and merging GeoJSON files in QGIS

texts:
I have got a dataset of GeoJSON files (26000 polygons). This is the first time for me to work on GeoJSON format and I need to import all rectangles into QGIS, and crucially, I need to dissolve any touching rectangles into a single polygon WITH A COUNT of how many rectangles were dissolved into that polygon.
I need to be able to sort the merged resulting polygons by a number of former layers because it indicates relative importance. First of all, I thought that I should convert the GeoJSON into shapefile, but the process is too slow for too many polygons.


texts: openlayers-3 linestring getLength not returning expected value

texts:
I am using getLength to retrieve the linestring length.

For the same segment:
1- when using google map measure tool, I get 228m
2- when using IGN geoportail measure tool, I get 228m
3- when I use e.feature.getGeometry().getLength() I get 330m

Here are the flat coordinates:
e.feature.getGeometry().getFlatCoordinates() :
[571382.4214041593, 5723486.068714521, 571593.8175605105, 5723741.65502785]

in 4326:
[5.132815622245775, 45.644023326845485, 5.134714626228319, 45.64562844964627]

When I check the coordinates position on either ol3 or google map, I get the same points. The difference must come from the calcul...

Did I miss something and I should not use the getLength method? Please give me some direction if you think this is not an issue.

I've done many tests and I get a 1.44 factor all the time between the distances. Maybe it's a unit problem. But I though the default unit in ol3 was meter.
I just divide the getLength result by 1.44 but I don't like it...


texts: How does QGIS Zonal Statistics handle partially overlapping pixels?

texts:
There are a couple Stackexchange posts asking this question, but neither of them were sufficiently answered:

How does Zonal Statistics work exactly?
How does Zonal statistics handle pixels that are on the border of the polygon?

I am writing Python scripts that perform essentially the same operations as Zonal Stats, and I have been using Zonal Stats to check my scripts' outputs.

I have noticed that when there are many cells in the polygon being queried, my results match perfectly, but there are discrepancies when there are only a few raster cells in the polygon. The following is an example of this. I am computing "count" for each polygon.

Exhibit 1:

The returned count is 0.45, which looks roughly correct if pixels can be treated fractionally rather than simply in/out.

Exhibit 2:

Similar result; count is 0.39.

Exhibit 3:

With a larger poylgon, now we see the count is exactly 2. Note that there are exactly 2 pixel centers contained in the polygon.

Exhibit 4:

Here, only one pixel center is within the polygon. The count is 1.008, though. If you exclude the polygon area from the pixel whose center is in the polygon, the remaining area is obviously greater than 0.008 of a pixel.

From these exhibits, it seems that Zonal Statistics treats these cases specially. It looks like to me that if there are less than 2 pixel centers within the polygon, it performs some kind of averaging such as polygon_area/raster_cell_area to return the count. I should note that I have 48 of these polygons, and all returned counts greater than or equal to 2 are integers; fractional counts are only returned for <2.

How are these fractional-pixel counts incorporated into the other statistics output by Zonal Stats?

How does Zonal Statistics work exactly? indicates strange behavior with max/min when using large raster cells, and my own testing returns mean values not congruous with simple "pixel center in/out of the polygon." To reiterate, I get exact agreement when the underlying raster contains many pixels within the polygon. I couldn't find any documentation about this, and it may be important when people have never had problems with Zonal Stats but use a coarse raster and remain unaware that there is a different behavior when there are only a few pixels within the polygon.


texts: Get the longest route from collection of polygons

texts:
const startPoint = [x,y];
const nodes = [[x,y], [x,y]....];

const polygons = [...];

const createConnectedPath = (startPoint, nodes, polygons) => {
 // The path should follow the center of each polygon.
 // The path should connect the startpoint with the node which have the longest path when following the path through the polygons.
 // The polygon is not required to be "connected" with another polygon. If so, connect the closest polygon with one another.
};

expectedOutput = [startpoint, [x,y][x,y][x,y]]; // A path

Is there any package that have this functionality out of the box? If not, what would be the correct way to approach this?

I'm thinking the first step would be to convert the polygon to lines. Then extend all lines until they intersect with another line. Then somehow connect each node with the closest line from the polygon. From there calculate the longest route and create a path of this route. Am I close? Or totally off?

It's a local coordinate system. Not GIS coordinates.

Totally new to this kind of logic; can anyone just point me in the right direction?

Python or JavaScript are preferred languages.


texts: Leaflet display only clicked element of a loaded map layer (polyline)

texts: I am creating a web mapping application that will return appropriate routes based upon ST_length query on a cartoDB dataset. So far, the code below successfully maps those routes with individual colours and clears any others that may have previously been selected:


texts:
I'm working a jQuery Mobile page with an ESRI jsapi map and a dialog which calls on a geocoding service to return an x/y which I plunk a graphic on the map and then center up on.

This site runs very nicely on my desktop..  However, on the phone the results do not populate the map.  I used the chrome debugging tools for android on my phone and my tablet to capture the error.

'Error: Invalid value for  attribute x="NaN"..

The madding thing is that the console and watch shows me the values for x and y are correct.

Is there a different place other than the point geometry for the mobile that I need to set the x and y coordinates?

Andy..

var map = new Map('ui-map-content', {
           basemap: "topo",
           autoResize:false,
           center: [-96.7968380, 32.7762360],
           zoom: 12,
           slider: true
       });

//Yep it was a success ('IknowwhereIam...)..
       function IknowwhereIam(data) {
           try {

               clearGraphics(map, pointGraphicLayer);
               var x, y;
               x = data.geometries[0].x;
               y = data.geometries[0].y;

               console.log("IknowwhereIam,  x: " + x + ",  y: " + y);

               var pt = new Point(x, y, map.spatialReference);

               var location = new Graphic(pt);
               console.log(location);
               pointGraphicLayer.add(location);

               map.centerAndZoom(pt, 16);

           }
           catch (e) {
               alert("problem inside the 'IKnowwhereIam' function: " + e.message);
           }
           finally {

           }
       }


texts:
I got a habitat map. The empty area is developed area and I want to fill it with grey color. How can I do that efficiently?

I tried to use clip but don't work.

I am using ArcMap 10.5.1 standard


texts: How to Convert 3D coordinates to 2D coordinates for R-Tree?

texts:
I am trying to implement a bounding box via a spatial index and get an error stating that I must either provide coordinates as either (minx, miny, maxx, maxy) or (x,y) for 2D indexes.

I've tried using the R-Tree example from this thread:  More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc

The first time I try to do a spatial index, I have the following value for a point ("Z" is the current value of "i" before my code breaks):

POINT Z (11.52951677300007 0.7729100360000416 -50000)

How can I gracefully get rid of the Z-coordinate to pass on as a 2D coordinate in my "for j" statement?  Code snippet is below:

    from rtree import index
    idx = index.Index() #Create an R-Tree index and store the features in it (bounding box)
    for pos, poly in enumerate(polygons):
        idx.insert(pos, shape(poly['geometry']).bounds)

        for i,pt in enumerate(points):  #iterate through the points
            point = shape(pt['geometry'])
            for j in idx.intersection(point.coords[0]):  #iterate through spatial index
                if point.within(shape(multi[j]['geometry'])):
                    discoveredBuildings.append(point['id'])


texts: The expression x(start_point($geometry)) doesn't return anything, how can I fix it?

texts:
I'm trying to get the coordinates of start and end of lines in a layer, but the expression x(start_point($geometry)) doesn't return anything, what's wrong?

Here's my code:

capa = "//atlantico2/ACTIVOS/0702 - EIA LINEA DE CONEXION ELECTRICA PA/04 Elaboracion/02 Ambiental/05 SIG/01 GDB/SHP/BIOTICO/ParcelasVirtuales_LN.shp"

layer = QgsVectorLayer(capa,'','ogr')
fields = layer.fields()

camposNuevos = ['X_Ini','Y_Ini','X_Fin','Y_Fin']
indices = []
for field in fields:
    for campo in camposNuevos:
        if field.name() == campo:
            indice = layer.fields().indexFromName(campo)
            indices.append(indice)
layer.dataProvider().deleteAttributes(indices)
layer.updateFields()

layer.dataProvider().addAttributes([QgsField('X_Ini',QVariant.Double),QgsField('Y_Ini',QVariant.Double),
QgsField('X_Fin',QVariant.Double),QgsField('Y_Fin',QVariant.Double)])
layer.updateFields()

exp = QgsExpression("'x(start_point($geometry))'")
context = QgsExpressionContext()
context.appendScopes(\
QgsExpressionContextUtils.globalProjectLayerScopes(layer))

with edit(layer):
    for f in layer.getFeatures():
        f['X_Ini'] = exp.evaluate(context)
        layer.updateFeature(f)

layer = iface.addVectorLayer(capa,'','ogr')

And this is the attribute table with the field needed in Zeros:

I'm using Qgis 3.4


texts: Deploying Projects ESRI.ArcGIS.ADF vs ESRI.ArcGIS.ADF.Local

texts:
I have a project which implements BaseTool and BaseToolbar.  For those of you who don't know already, the name of the dll containing these base classes has changed from ESRI.ArcGIS.ADF to ESRI.ArcGIS.ADF.Local between ArcGIS 9.x and ArcGIS 10.x.  However, the namespaces I need have remained the same.

To deploy this toolbar, I have been switching back and forth between dev machines running ArcGIS 9.x and 10.x, updating the reference, and then recompiling.  I don't have to change a single line of code.  Although this is a fairly easy process, I then have to create and manage separate installers for 9.x and 10.x versions.

I've been banging my head against the wall trying to come up with a more elegant solution.  Something along the lines of changing the reference at install time?  Does anyone have any suggestions?


texts: Finding out whether Leaflet is using Canvas or SVG?

texts:
Leaflet maps use either an SVG or Canvas renderer.

How do I figure out which one a created leaflat map is using?


texts: Python Automated Data Driven Pages - Labeling and Symbolizing Errors

texts:
I am trying to use a python script to automate the export of the data driven map document.

However when I use the python script to run the DDP export some of the labels are not exported and the symbols that have a rotation angle based on a attribute in the layer table is not rotating.

Other than the labels and the symbols the python script to automate the exports works fine.

When I run the DDP manually (through the MXD) all the labels and symbols are exported correctly.  Have anyone else encountered this problem?


texts:
I want to change the default styles for vector layers in my Leaflet version 0.7.3. I can't find anything about that in in the documentation, so I made this hack:

map.on('layeradd', function(e){
  if(typeof e.layer.setStyle === 'function'){
    e.layer.setStyle(myApp.Map.Styles['default']);
  }
});

I don't find that a sustainable way of setting a default style, since I may sometimes want e.g. to overrule the default style or some other scenario.

Any suggestions?


texts: determine the closest points on the network?

texts:
In ArcGIS, I have only a point layer relating to a hydrography network's pour points. Now I want to write code in VBA that selects a point from this layer, then highlight the closest points to this point. One of this points must be downstream of the desired point and another point must be upstream.
In this work, I must measure the distance between two points on the network, not a direct distance between two points. Also both points are related to a layer.
Also the program must be able to determine one of points in downstream and another points in upstream.
I don't know how to determine the closest points on the network.
Please help me again.


texts:
I am going to ask this question better:

What happened with NVIZ?

I do not know the exact version that nviz stopped working but it seems around version 2 of QGIS. There are two ways to get to nviz in QGIS 2.4

1) Processing toolbox - I double click on the nviz which is under visualization which is under GRASS commands and receive this error:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\gui\ProcessingToolbox.py", line 181, in executeAlgorithm
message = alg.checkBeforeOpeningParametersDialog()
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 511, in checkBeforeOpeningParametersDialog
msg = Grass7Utils.checkGrass7IsInstalled()
 File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 361, in checkGrass7IsInstalled
cmdpath = os.path.join(path, 'bin', 'r.out.gdal.exe')
  File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\ntpath.py", line 96, in join
assert len(path) > 0
TypeError: object of type 'NoneType' has no len()

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

QGIS version:
2.4.0-Chugiak Chugiak, 8fdd08a

Python path: ['C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\\processing', 'C:/PROGRA~1/QGISCH~1/apps/qgis/./python', u'C:/Users/Bobby Pescatore/.qgis2/python', u'C:/Users/Bobby Pescatore/.qgis2/python/plugins', 'C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGISCH~1\\bin\\python27.zip', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\DLLs', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\plat-win', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\lib-tk', 'C:\\PROGRA~1\\QGISCH~1\\bin', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\PIL', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\jinja2-2.7.2-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\markupsafe-0.23-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\pytz-2012j-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\six-1.3.0-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\wx-2.8-msw-unicode', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\xlrd-0.9.2-py2.7.egg', 'C:\\PROGRA~1\\QGISCH~1\\apps\\Python27\\lib\\site-packages\\xlwt-0.7.5-py2.7.egg', 'C:\\Users\\Bobby Pescatore\\.qgis2\\python\\plugins\\MetaSearch\\ext-libs', 'C:\\PROGRA~1\\QGISCH~1\\apps\\qgis\\python\\plugins\\fTools\\tools', 'C:/Users/Bobby Pescatore/Desktop/Batimetria']

2) From GRASS Tools - In the module list you can type nviz and if you try to initiate it from there you get this error:

Any help on this issue would be great. I get that the whole thing is open source so there is no dedicated support but I can't imagine that nobody has solved this yet.

For anybody who needs this functionality immediately, I have learned from the research I have done that if you work directly from GRASS and not through GIS this should work. I haven't tried yet but plan to tomorrow. Will update then.


texts: Rendering huge vector data whose geometry is fixed but attributes are changing every 2 minutes?

texts:
I want to create an application displaying road network with traffic details.

I am using OpenLayers 3 and GeoServer 2.14 for this purpose. My data (stored in PostgreSQL+ PostGIS) has around half a million geometries (road network) which is published in Geoserver through SQL view.
Geometry of the data is fixed but its attributes related to traffic are changing every 2 minutes. The approach we are currently using is slow. Is there any other option available to render data faster.

We have road network data with unique link id assigned to each road segment. There is one more column which stores speed data. According to the speed we display color of the road segment (red, orange, green etc.) on map using openlayers. This speed info updates every 2 min while road geometry is fixed. We have indexes built on the table. Currently we are fetching these details using SQL view using Geoserver where we are calling a function which fetches required data. For better performance we are also passing bbox (via view params) of tile into the function to fetch data of specific bbox only. This has improved the performance when data count is approx 2.5 lacs. Now the data count has just doubles and this approach is not working as expected.


texts: How to create a polygon with a circular hole in openlayers 3?

texts:
To highlight an area on my map I need to make a circular hole in a polygon. The rectangle polygon covers the whole map, the circle hole should act as a spot on my area to highlight.

But actually I only manage to make a rectangle hole, for example using Polygon's appendLinearRing function:

mypolygon.appendLinearRing(new ol.geom.LinearRing([
    [239850, 5069850],
    [239850, 5070150],
    [240150, 5070150],
    [240150, 5069850],
    [239850, 5069850]
]));

The problem is that I can't find out how to add "Circle" as optional geometry (opt_layout) to ol.geom.LinearRing. Any hints?


texts: It is possible to enable CORS in an old Windows GeoServer version (2.6.1) using Jetty 6.1.8?

texts: Of course I found some info about how I have to do this. But I cannot find the jetty-servlet package matching this old version.


texts: How to find/replace element in field with Field Calculator

texts:
I just want to perform the following operation done in Matlab into the Field Calculator using Python:

p(p==0) = 0.001 ; p(p == 1) = 0.999;
tVar = .5*log(p ./ (1 - p))

the operation consists of a Logit transformation:

Find 0's and 1's and replace with .001 and .999 respectively, and
Perform the second calculation with logarithms

I've tried the following for 1) in Python to no avail:

Pre-Logic Script Code:

def iLogit(a):
  return [.0001 if x==0 else x for x in a]

Any help please?


texts: GeoExt Legend Modification To Include Feature Count

texts:
GeoExt has Legend capability but is is quite basic which shows only feature types like in example http://api.geoext.org/1.0/examples/legendpanel.html. How could it be modified display legend with count(no of values)?

   I am using geoserver to serve map and using geoext. The solution I am looking for is to add this feature into legend panel (may be through modification of javascript source) and not make my own solution. Any path, direction code for this would be appreciated.


texts: GRASS7 r.water.outlet doesn't run correctly in QGIS 2.14.2

texts:
I am a new user of GIS software and I am currently facing an issue with the GRASS7 tool `r.water.outlet. I am a beginner with both Linux and GIS. I am using Ubuntu 14.04 with QGIS 2.14.2 and GRASS 7.0.3 I downloaded and installed the software by following the instructions of the official QGIS webpage.

I am using GRASS7 r.watershed to create a "Drainage" raster that I use as an input for r.water.outlet. The second input is the coordinates of vectorial points I want to use to delineate my basins (it would be even better if I could use several points to do so because I need 5 different basins).

This is my r.water.outlet interface:

First thing I noticed: I can't enter decimal numbers, so I lose accuracy for my input points, and the interface is different from the one in the second image for the tool r.water.outlet for my version of QGIS.

Then when I run the process (In my QGIS 2.14.2, the processing plugin version is 2.12.2), I obtain the log:

*Démarrage de l'algorithme r.water.outlet - Watershed basin creation program. ...
g.proj -c proj4="+proj=lcc +lat_1=44.1 +lat_0=44.1 +lon_0=2.337229166666667 +k_0=0.999877499 +x_0=600000 +y_0=3200000 +a=6378249.2 +b=6356514.999904194 +units=m +no_defs"
r.external input="/home/kyle/Bureau/ERM/Data_ERM/GISProjectFolder/Rasters/Drainage_Direction_Watershed.tif" band=1 output=tmp1462467882073 --overwrite -o
g.region -a n=1959000.0 s=1843500.0 e=1038600.0 w=923100.0 res=300
r.water.outlet input=tmp1462467882073 easting="989927" northing="1891138" output=output881924edb38d42189f35bcf5166a0da3 --overwrite
g.region raster=output881924edb38d42189f35bcf5166a0da3
r.out.gdal -c createopt="TFW=YES,COMPRESS=LZW" input=output881924edb38d42189f35bcf5166a0da3 output="/home/kyle/Bureau/ERM/Data_ERM/GISProjectFolder/test_water_outlet.tif"
DÃ©marrage du SIG GRASS ...
ExÃ©cution de '/home/kyle/.qgis2//processing/grass7_batch_job.sh' ...
ATTENTION: Datum non reconnu par GRASS et aucuns paramÃ¨tres trouvÃ©s
Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default
Projection information updated
ATTENTION: Datum non reconnu par GRASS et aucuns paramÃ¨tres trouvÃ©s
ATTENTION: Over-riding projection check
Reading band 1 of 1...
r.external terminÃ©. Link to raster map created.
Description:
Creates watershed basins from a drainage direction map.
Keywords:
raster, hydrology, watershed
Usage:
r.water.outlet input=name output=name coordinates=east,north
[--overwrite] [--help] [--verbose] [--quiet] [--ui]
Flags:
--o Autoriser la sortie Ã  Ã©craser les fichiers existants
--h Afficher le rÃ©sumÃ© d'utilisation
--v Sortie du module en mode bavard
--q Sortie du module en mode silence
--ui Force launching GUI dialog
Parameters:
input Name of input drainage direction map
output Name for output watershed basin map
coordinates Coordinates of outlet point
ERROR: r.water.outlet: Sorry, is not a valid parameter
ERROR: r.water.outlet: Sorry, is not a valid parameter
ERROR: Required parameter not set:
(Coordinates of outlet point)
ERREUR :Couche raster non trouvÃ©e
ERREUR :Raster map or group not found
ExÃ©cution de '/home/kyle/.qgis2//processing/grass7_batch_job.sh' terminÃ©e.
Nettoyage des fichiers temporaires ...
Converting outputs
Charger les couches de résultat
The following layers were not correctly generated.
Basin
You can check the log messages to find more information about the execution of the algorithm*

Does anyone have any idea what could be the problem?  If you need further information, do not hesitate to ask me.


texts:
Last couple of days i am working with creating thematic map using MapBasic. I am facing some difficulties while doing scripting and compilation.i am getting error due to Dynamic Window Id. does any code or solution to set window id? current window id i got while doing thematic map manually using mapinfo tools. if i used that window id in to my code it working, but only thing is that i have to open that window. if i closed that window and run my code with that window id it showing me incorrect window id.

anybody can tell me..how to solve this issue

thanks and regards


texts: OpenLayers map with Nominatim lookup function not working

texts:
I am trying to create a simple map with a search function, my minimal example is here: http://jsfiddle.net/TimLucas/vbaupe30/4/

For some reason the coords transformation after the script has received data from nominatim from EPSG 4326 to EPSG 3857 does not work well, even though it does work in an earlier instance in the script (right at the beginning).I have stepped through the function and the x - coord seems to be fine but the y - coord is wildly off.

Has anyone encountered a similar problem?


texts: Turn detached coordinate values into linestring in GeoPandas


texts: Grouping points based of user-provided distance in ArcGIS for Desktop?

texts:
I have a point layer and I need to group/aggregate the points within a 1km proximity into 1 single point to create a new point layer.

I was trying to use Aggregate point and then Feature to point, but it eliminates the single points I have in case there are no other points in  1km proximity, but I need to be able to save those single points too. It also groups only three or more points.


texts: Adding WMS group layer (only) to ArcMap?

texts:
I am trying to add a WMS group layer to ArcMap, without having the sublayers being listed in the TOC.
I will give a little background here: We are sourcing images from Nearmap.com, and their WMS data is served as follows:

NearMap
|
Australia

 Adelaide
-----2009-10-19
-----2009-10-29
 Perth
-----2010-10-12

As per Nearmap, the group layer link (Adelaide or Perth shown above) always points to the most recent image avaible, while the sub-layers present date specific images (one of them would of course be the latest one).
If this is true, it releives us from updating and re-publishing our maps every time a new image is available from Nearmap (which happens at least once a month).
When I try to add the group layer in Arcmap, it adds all the sub-layers along with it, and I have to choose a specific date image to publish in our web maps.

Does anyone here has experienced the same behaviour in other GIS software, open-source or otherwise?

I was informed by the Nearmap guys that "other" open source software have the "ability" to add group layer directly...any takers for that?


texts: Add a GeoJSON feature collection of linestrings to OpenLayers map

texts:
I have a GeoJSON object of line string coordinates, geobj
I am attempting to display them on my OL map object.
  var layer = new ol.layer.Vector({
    source: new ol.source.Vector({
      features: (new ol.format.GeoJSON()).readFeatures(geobj)
    }),
    style: new ol.style.Style({
      stroke: new ol.style.Stroke({
        color: "black",
        width: 10
      }),
    })
  });

  layer.setZIndex(110);
  map_obj.addLayer(layer);

But nothing appears.
I have double checked the coordinates are correct, by plotting them individually on the map.
EDIT: Ran some tests, it's clear the coordiantes are fine. The issue is to do with how I generate the geojson object, I've done this using a loop - something is not quite right with the final product.
var geojsonObject = {
    'type': 'FeatureCollection',
    'crs': {
      'type': 'name',
      'properties': {
        'name': "EPSG:27700",
      },
    },
    'features': []
  };

  $.getJSON(fn, function( data ) {

    data.features.forEach(calculateTraj);

    function calculateTraj(item){

      const MvtSpeed = item.properties.MvtSpeed;
      const MvtDirec = item.properties.MvtSpeed;

      const startLon = item.properties.LonTrajCellCG[0];
      const startLat = item.properties.LatTrajCellCG[0];

      const kmToLon = 1/(111.32 * Math.cos(startLat * (180/Math.PI)));
      let kmToLat   = 1/110.57;

      const distKMin15m = (MvtSpeed * 60 * 60)/1000;
      const alpha = (MvtDirec * (180/Math.PI));

      const deltaLon = (distKMin15m * Math.cos(alpha) * kmToLon);
      const deltaLat = (distKMin15m * Math.cos(alpha) * kmToLat);

      const endLon = parseFloat((startLon + deltaLon).toFixed(4));
      const endLat = parseFloat((startLat + deltaLat).toFixed(4));

      var t1 = ol.proj.fromLonLat([(startLon), (startLat)], ukngProj)
      var t2 = ol.proj.fromLonLat([(endLon), (endLat)], ukngProj)
      geojsonObject.features.push(
        {
          'type': 'Feature',
          'geometry': {
            "type": "LineString",
            "coordinates": [
              [t1],
              [t2] ],
          },
        });
    }
  });
  return geojsonObject;
}


texts: Clustering Geoserver with shared session

texts:
I have succeed in clustering geoserver on windows using three Apache tomcat instances and Apache server as load balancer. I also shared the data directory but I have two issues:

when I try to login in the geoserver through the load balancer, I failed to login although the username and password are correct and I didn't get any error message. And when I logged in geoserver through any instance of tomcat, the geoserver through the load balancer didn't login.
when I upload a layer on an instance of geoservers, the others didn't sync until I reload the configuration of each one.


texts:
How can I flip labels on my lines upside down? I am styling some OSM data and the cartocss is currently

#streets{
text-size: 14;
text-name: [NAME];
text-face-name: 'ATRotSemSanLigIta Regular';
text-fill: #000;
text-halo-fill: fadeout(white, 40%);
text-halo-radius: 2;
text-placement: line;
text-dy: 12;
}

I think its via text-orientation but cannot get it to fire