texts:
I use a compiled version of GDAL, including OGDI lib :
-   GDAL 2.0.1
-   OGDI 3.1.6

I can open DNC Charts and retrieve layers / features.
DNC Charts use VPF format.

Problem : I’m not able to retrieve notes which can be found in NOTES.RAT and NOTES.NJT, like explained in this VPF specification :

VECTOR PRODUCT FORMAT (VPF) products General Specification

3.16.4.2.1.1 Notes related attribute tables

Have you got ideas about this question ?


texts: Removing tree stems from a point cloud in R

texts:
I am studying the foliage density of forest plots below 2 metres in height using mobile terrestrial laser scanning. I want to filter out the tree stems from the point clouds so that the point density estimates do not include them. I have been using the lidR and TreeLS packages to process my point clouds. So far I have normalized the clouds and classified the tree stems using example code given in the TreeLS package documentation:
tls1 = readTLS("tls1.las") %>%
tlsNormalize %>%
tlsSample
map = treeMap(tls1, map.hough())
tls2 = treePoints(tls1, map, trp.crop(circle=FALSE))
tls3 = stemPoints(tls2, stm.hough(pixel_size = 0.03))

This successfully identifies the tree stems.
To filter the stems out, I have tried using the filter_poi() function:
tls4 <- filter_poi(tls3, Stem = "FALSE")

There is a column in the point cloud dataset entitled "Stem" with TRUE/FALSE values so I thought this would work, but I get this error:
Error in lasfilter_(las, lazyeval::dots_capture(...)) :
  `conditions` must be logical.


texts: How to create an Erdas batch list file (.bls) using Python?

texts: How can I generate the Erdas batch list file (.blf) using Python?  For example, I want to run a 5x5 median filter in the Batch Command Editor using a subset of images within a folder.  Erdas provides a way to load lists of files, although I am not sure how to generate these outside of Erdas.


texts: ArcGIS Pro Python Interpreter Setup with Visual Studio Code

texts:
I'm fairly new to Stack Exchange and not very well-versed in IDE's.  I'm struggling to set up the ArcGIS Pro Python interpreter with Visual Studio Code.  I have ArcGIS Pro 2.4.3 and VSC 1.44.2.  I'm trying to follow the advice on the below thread, but I cannot ask questions in that thread, only add answers.

Original Thread

So here's my question.  I cannot find the ""+ Custom..." and select "Configure" location where I add the Python interpreter, as mentioned in the thread above.  In VSC I can go ctrl+Shift+P and it shows me the current Python interpreters that the program already knows, but I cannot add a new location.  In File-->Preferences-->Settings under Workspace I can search for Python and there are 147 settings for Python.

I tried adding the C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\python.exe path under Python Path and adding C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3 under Pipenv Path but this hasn't changed anything in terms of the program understanding that I'm trying to add a new Python interpreter option.

There are so many settings for Python I don't even know where to start or if I am even changing the right settings.  It's not very intuitive.  Please help!


texts: How does multiuser work with versioned ArcGIS?

texts:
I am new at versioned editing and am confused a bit and have some questions. For example I have five editor users and versioned dataset or featureclass. So,

Should each editor create their own version under DEFUALT version? (named John, Mike, Bob,....)
Should we create only one version named Editors and all editor users use it? And in this case is there any possibility of locking each other?


texts: Viewing polyhedrals through the new 3D viewer QGIS 2.99 (3.0)

texts:
I've installed QGIS 2.99 to test its new features, particularly the 3D viewer, when I try the following I only see a square, not a 3D cube.

In PostGIS I make a square and extrude it

-- make a square
CREATE TABLE public.square1
(
    id serial,
    geom geometry,
    CONSTRAINT square1_pkey PRIMARY KEY (id)
);

INSERT INTO public.square1 (geom)
SELECT ST_MakePolygon(ST_GeomFromText('LINESTRING (0 0, 0 100, 100 100, 100 0, 0 0)')) AS geom;

I add the following sql via the db manager and add as a layer.

SELECT
ST_SetSRID(ST_Extrude(ST_Buffer(square1.geom, 0),0,0,10), 27700) AS geom
FROM Square1

But I only see 2D geometry, no 3D geometry.
I had to set a srid, QGIS would not allow no projection or ESPG:4326

I test the polyhedral for its volume, this is correct

ST_volume(ST_Extrude(ST_Buffer(square1.geom, 0),0,0,10) ) AS volume
FROM Square1

Where have I gone wrong?


texts: Registering database to ArcGIS for Server gives null error?

texts:
I'm trying to register a database to ArcGIS for Server 10.2.0 from the web server manager using sde file from ArcMap.  The connection works fine in ArcMap.
But when registering I'm only getting a "null" error.

I can register the database with SDE user but not with the user that has schemas with map services.

We are trying to do this after copying the database from production environment to a development environment.  In production everything works well.


texts: PostGIS: calculate intersection, then split intersection by attribute?

texts:
I'm using Postgres 9.6 and PostGIS 2.3.2. I want to find the intersection of two tables, then the area of the intersection for each of the attributes from the first table. My first table is called usgovt:

 category   | numeric                     |
 geom       | geometry(MultiPolygon,5070) |

And my second table is called forest:

geom       | geometry(MultiPolygon,5070) |

I want to calculate the area of the intersection that is in each category.

The tables are large, so I am keen to do this efficiently. This is my first attempt but it is VERY slow (been running for many hours):

SELECT g.category AS category,
         ST_Area(ST_Transform(g.geom, 4326)::geography) AS area,
         ST_Area(ST_Transform(ST_Intersection(g.geom, r.st_union), 4326)::geography) as overlap
    FROM (SELECT
         category,
         ST_Union(geom) AS geom
       FROM usgovt
       GROUP BY category) g, forest f
    WHERE
      ST_Intersects(g.geom, f.geom);

I think I'm doing it wrong.

What is the most efficient way to calculate the intersection once, THEN look at the area of the intersection in each category?


texts: Split points evenly to neighboring polygons?

texts:
I have a crime location data set, containing geocoded crime location points that have been geocoded to streets and street interactions. I also have a census tract shapefile. What I would like to do is to allocate these points to census tracts (polygons).

The problem is that many crime points fall right on or very close to the boundaries of tracts. I would like to split these points evenly to neighboring tracts.

For example, in the picture attached, there are 6 selected points right on or within 30 feet of the boundary of two tracts, I would like to count 3 points toward each tract.


texts: Combining information from raster and vector layers in QGIS?

texts:
in QGIS I have 1 raster-layer (elevation) and 1 vector-layer (administrative units).

Now I would like to be able to query one of the administrative units and see the mean and min/max of the elevation within its area.

Which QGIS tool can do that?


texts: find lat,long inside polygon with regular interval of time using mysql

texts:
I have table latest location there i have lat and long with device id.
 one more table with fencing name,there i have polygon points.
 so now how i can check point lies inside polygon.
 i need this with regular interval of time and storing that to new table.
 so how to find this with MySQL function,event or with procedure?


texts: ogr2ogr ST_Difference operation on OGRVRTDataSource

texts:
I'm trying to subtract from one shapefile the portions of another which overlap using the following:

land.shp: subunits.shp lakes.shp
  COOKIE=subunits CUTTER=lakes; \
  echo '<OGRVRTDataSource><OGRVRTLayer name="'$$COOKIE'"><SrcDataSource>'$$COOKIE'.shp</SrcDataSource></OGRVRTLayer><OGRVRTLayer name="'$$CUTTER'"><SrcDataSource>'$$CUTTER'.shp</SrcDataSource></OGRVRTLayer></OGRVRTDataSource>' | \
  ogr2ogr -f "ESRI Shapefile" \
    land.shp /vsistdin/ \
    -dialect sqlite \
    -sql "SELECT ST_Difference($$COOKIE.geometry,$$CUTTER.geometry) AS geometry FROM $$COOKIE,$$CUTTER"

lakes.shp: ne_10m_lakes.shp
  ogr2ogr \
    -t_srs EPSG:4326 \
    -clipdst -92.956317 49.113085 -74.952721 40.992926 \
    lakes.shp \
    ne_10m_lakes.shp \

subunits.shp: ne_10m_admin_0_map_subunits.shp
  ogr2ogr \
    -where "ADM0_A3 IN('CAN')" \
    -t_srs EPSG:4326 \
    subunits.shp \
    ne_10m_admin_0_map_subunits.shp

(Makefile, sorry about the long lines)

The result I'm getting is identical to the input, which is what should result if the 2nd file (CUTTER) is empty. It's not, however, and as a sanity check, when I reverse the order of the first line (COOKIE=lakes CUTTER=subunits) and leave everything else the same the result is what I'd expect: the US part of the lakes minus everything else. That suggests to me that the query itself is fine. Is there something else I might be missing?

Here are the contents of the .prj files. subunits.prj:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

lakes.prj:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]

The only difference being that lakes.prj uses floats for all numbers, though they're otherwise the same.

TL;DR: I want the area of the polygon around the Great Lakes to follow the shoreline, not the int'l. border.


texts: Can you dynamically set print options using ArcGIS JS API?

texts:
I'd like to create an option in my web application where a user can click on a print button and a list of options appears:

Page template (already available via esri.dijit.print)
Map Title
Notes or Custom Text

No. 1 is known to be available but I can't find any good resource on how to be able to dynamically change these inputs. Any ideas?

Code I'm working with:

    dojo.require("esri.dijit.Print");

function initPrinter() {

    var printTemplates = [{
        label: "Map Only",
        format: "png32",
        layout: "MAP_ONLY",
        exportOptions: {
          width: map.width,
          height: map.height,
          dpi: 96
        }
      }, {
        label: "Portrait",
        format: "PDF",
        layout: "Letter ANSI A Portrait",
        layoutOptions: {
          titleText: "my layout",
          authorText: "GIS Team",
          copyrightText: "Copyright 2013",
          scalebarUnit: "Miles"
        }
        }, {
        label: "Landscape",
        format: "PDF",
        layout: "Letter ANSI A Landscape",
        layoutOptions: {
            titleText: "my layout",
            authorText: GIS Team",
            copyrightText: "Copyright 2013",
            scalebarUnit: "Miles"
            }
      }];

    // print dijit
    var printer = new esri.dijit.Print({
        map: map,
        templates: printTemplates,
        url: "path/to/my/Export%20Web%20Map%20Task"
        }, dojo.byId("printButtonDiv"));

    printer.startup();

}


texts: How to use GRASS command "r.stream.distance" in Qgis 2.16.2?

texts:
I like to figure out the distance from each cell to the outlet along the flow path.
Maybe you can suggest some other opportunities/commands to solve my problem.
It would be nice, if someone can help me because i am a newcomer.


texts: How can a 9.2 user access a geodatabase created in 9.3?

texts: I have a colleague that uses 9.2 and has access to our server in which I have multiple geodatabases that were created in 9.3.  Is there any way for her to access the data?


texts: Clipping one shapefile file with 100 different shapefiles (Arcmap 10.1) and retrive 100 diferent outputs using ModelBuilder

texts:
I have tried without success to perform this task with ArcMap 10.1. I´m struggling with ModelBuilder.
I have one shapefile (mean precipitation in one country) and I need to clip this shapefile with 100 different shapefiles (watershed areas). I can´t use Union to create a single shapefile from the 100 different shapefiles, because some of them overlap.
Can this kind of task can be automated with Iterate Feature Classes inside ModelBuilder and, if so, where can I find an example?


texts: WMS layer not displaying in cesium globe

texts:
Im downloaded cesium 1.23 library from cesium site and added geoserver wms layer on the cesium globe, but its not displaying.

my query has mentioned below.

<!DOCTYPE html>

  <title>Hello World!</title>

How to add wms layer over the cesium globe?


texts: How to quantify the similarity of spatial patterns?

texts:
I am not a GISer, but I do believe a GIS (more specifically, spatial statistics) algorithm could help with the task I am carrying on. If posting a question that has already been posted (by me) on stackoverflow here does break the rule, I would delete this question.
I have several 10000*10000 grids (raster data) denoting the same tissue in mice, each representing the expression pattern of a specific type of gene (each grid-point has a value between 0 and 100). The task I am working on now is to quantify the similarity between each spatial expression patterns of genes. Since I want to take spatial information into account, I am thinking about either applying image similarity comparison algorithms from image processing (or computer vision), or spatial similarity comparison algorithms from GIS (spatial statistics). Among the two I believe spatial statistics might be more appropriate.
Sample images:

Please refer to my original question on stackoverflow for more details:
https://stackoverflow.com/questions/65912256/how-do-i-quantify-the-similarity-of-spatial-patterns


texts: Exporting a globe projection of a TIFF as a PDF

texts:
I'm trying to export a globe projection of the 'Natural Earth II with Shaded Relief, Water, and Drainages (large)' as a pdf but I keep getting the following error: Could not transform bounding box to target CRS

I'm using the following custom CRS:
+proj=ortho +lat_0=29 +lon_0=98 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs

I have also exported the 'Countries Map', the 'States, Provinces map' and the 'Graticules (10) map', all of which at least export using those CRS. The latter however does export some weird horizontal lines that should not be there.

I came across a few tutorials such as this one that seem to suggest that my problem has to do with improper clipping/projection of the layer. However, none of the solutions seem to work for me because the suggested Python scripts keep coming up with errors.

I'm working in QGIS 3.4.4.


texts: Create a new list from a raster list in R

texts:
I am very new to R, so this is tough for me to explain. But, I have a raster list made up of 3872 elements (or separate rasters, I am guessing) and I am hoping to make a new smaller raster list by pulling out some of the rasters in the larger list based on their name.

The data I am using are raster .bin files for sea ice concentration values from NSIDC. Here are my steps to uploading the .bin files into R and creating the raster list itself:

#making the raster
pixel <- 25000 #pixel dimension in meters for both x and y
xMin <- -3837500 #From NSIDC: ulxmap -3837500
xMax <- xMin + (pixel*304)
yMax <- 5837500 #From NSIDC: ulymap 5837500
yMin <- yMax - (pixel*448)
r <- raster(nrow=448, ncol=304, xmn=xMin, xmx=xMax, ymn=yMin, ymx=yMax)
projection(r) <- '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs'

# making the function to loop over the files
bin_raster <- function(ymd){
filename <- paste("data/","nt_", ymd,"_v1.1_n.bin", sep="")
con <- file(filename, 'rb')
x <- readBin(con, "raw", 300)
x <- readBin(con,"int", size=1, signed=FALSE, 150000)
close(con)
x[x==251] <- 1
x[x>251] <- 0 # set all non-ice pixels to 0% ice cover
x[x<251] <- x[x<251]/250 #scale ice cover to 1
rr <- setValues(r, x) # place result in raster

#changing file names to just the date
geoTiffFile<- paste('geoTiffs/', ymd,'.tif',sep='')
rr <- as.list(rr)
}

# writing the loop
# making a vector of all dates for every .bin file (extract all file names, then isolate date)
files <- list.files("data") # make a vector of all file names
head(files)
dates <- str_sub(files, 4, -12) # drop all text except for the date

# creating an empty list
raster_list <- rep(NA, length(dates))

# looping over all the files and filling the list
for(i in seq_along(dates)){
raster_list[i] <-   bin_raster(dates[[i]])}

# renaming the files so that only the date remains
list_names <- str_sub(dates, 0, -5)
names(raster_list) <- list_names

So now I have a raster list (raster_list) of all my files that are named by their date. How do I create a new list and choose specific dates to put into it?


texts: Plug in / API for getting sea routes with leaflet

texts:
I am developing a maritime routing application prototype with Leaflet.js for my company Deep Blue Globe and using the windy interface.

The problem is since windy uses Leaflet version 0.7.7 there are not many routing APIs compatible with it and the one which is compatible doesn't provide sea routes.

So if anyone could suggest me something, any alternative it would be really helpful. Getting sea routes is of major importance here.

We want an interactive interface like windy (any alternatives to windy also welcome).


texts: Openlayers - parse WFS features to create Attribute Table

texts:
I am trying to parse this WFS to create an Attribute Table with HTML having the fields as headers and their values as attributes.

I saved the GetFeature response in a variable (features) which I'd like to use for my purpose, but i have no idea how to go on.

I used the readFeatures() method which should return the features as an Array, but from the dev console I see a 0 length Array which I guess is not really what I need.

Here is my code (also the jsfiddle, please look in the console to see the result):

import WFS from 'ol/format/WFS';

export function testWFS() {
    var format, xmlDoc;

    var url = 'https://www.wondermap.it/cgi-bin/qgis_mapserv.fcgi?map=/home/ubuntu/qgis/projects/Demo_sci_WMS/demo_sci.qgs&SERVICE=WFS&REQUEST=Getfeature&TYPENAME=domini_sciabili,impianti_risalita,piste_sci';

    var xmlhttp = new XMLHttpRequest();
    xmlhttp.onreadystatechange = function() {
        if (xmlhttp.readyState === 4) {
            format = new WFS({
                featureNS: 'http://www.qgis.org/gml',
            });
            xmlDoc = xmlhttp.response;
            console.log(xmlDoc)
            var features = format.readFeatures({
                source:xmlDoc,
                dataProjection:'EPSG:32632',
            });
            console.log(features)

        }
    };

    xmlhttp.open("GET", url, true);
    xmlhttp.send();
}


texts: Spatial Join Polygon to Polygon with Centroids Inside features

texts: I have a parcel feature class with +5M features. I want to spatially join it to a feature class of regions to check if the regional code matches that in the parcel data. I've tried using "have their centroids within", but there are many cases where multi-part parcel features do not join correctly because the centroid is outside the features. Is there a way to spatially join with centroids inside features? The only other place I've seen this is when using Feature to Point, and checking "Inside"


texts: Using arcpy.da.SearchCursor to sort field in attribute table by ascending order

texts:
I am using ArcMap 10.6 with Python 2.7.14 with ArcPy.
I want to sort a field named Seq_ID in my attribute table by ascending order using ArcPy.
I believe I have to use arcpy.da.SearchCursor.
Here is my code snippet:
for row in arcpy.da.SearchCursor(points, fields, sql_clause=(None,'ORDER BY Seq_ID ASC'):
     arcpy.AddMessage("")
else:
    arcpy.AddMessage("")

However, after using this snippet above, my table is not sorted.


texts: Merging TIFF images based on any of the single images

texts: I have lot of tiff images for a small area. I need to merge all the images based on any one of the image(for color ref). Actually we merge the data in QGIS but all images all slightly variation in their RGB color.


texts: Adding jpeg. files and not displaying in colour

texts:
When I try to add a jpeg file into my project on Arc it displays only in black and white.

I have checked the file and it is definitely in colour.

Does anybody know how I can add such files in the project in colour?


texts: How can I create a Map Document (including layers, symbology info etc.) for DotSpatial just like we have MXD File in ArcGIS?


texts: GEE Tile error: Expected a homogeneous image collection, but an image with incompatible bands was encountered:

texts:
I'm having the problem of Tile error: Expected a homogeneous image collection, but an image with incompatible bands was encountered. Since number of total S-2 bands have changed during time, what can we do to solve this, when we don't even need to use that 2 new ones?
std all: Tile error: Expected a homogeneous image collection, but an image with incompatible bands was encountered:
  First image type: 21 bands ([B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12, AOT, WVP, SCL, TCI_R, TCI_G, TCI_B, QA10, QA20, QA60]).
Current image type: 23 bands ([B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12, AOT, WVP, SCL, TCI_R, TCI_G, TCI_B, MSK_CLDPRB, MSK_SNWPRB, QA10, QA20, QA60]).
          Image ID: null

var datasetall = ee.ImageCollection('COPERNICUS/S2_SR')
                 .filterDate('2017-05-01', '2020-09-30')
                  .filterBounds(polygon);

var medall = datasetall.median();
var resultall = medall.select('B4', 'B8');
var ndviall = resultall.normalizedDifference(['B8', 'B4']);

var palette = ['FF0000', 'ffa500', 'FFFF00', '4CA64C', '207401'];

var stdall = medall.reduce(ee.Reducer.stdDev());

Map.addLayer(stdall.clip(polygon), {min: 0, max: 0.2, palette: palette}, 'std all');

var ndvi_all = datasetall.map(
    function(img) {
         return img.normalizedDifference(['B8','B4'])
                  .rename('NDVI')
                  .copyProperties(img, ['system:time_start']);    });

print(ui.Chart.image.series(ndvi_all, polygon, ee.Reducer.stdDev(), 10, 'system:time_start' ));


texts: I want to cut a line shape into parts by polygon areas. No problem so far. But the polygons overlaps each other at the border. I need one line result for the whole area of every polygon. Any ideas for an automatic solution? Gdal ogr can't do it because of self intersecting polygons.


texts: ArcMap Python script format string padding to line up multiple lines of string with variable lengths

texts:
I'm using the advanced Python expression to build a label and am having trouble with the padding function in the format function as follows (rows is an arcpy SearchCursor output):

for row in rows:
     if row.DETECTFLAG=='N':
        texttoadd = ""
        width1= 60-len(row.CHEMICALNM)
        texttoadd = texttoadd + "{0}".format(row.CHEMICALNM) +"<BOL><ITA>{0:>{2}} {1}</ITA></BOL>".format(row.RRESTXT,row.LABQUALIFI,width1) + '\n'
        myText3 = myText3 + texttoadd
     elif row.FLAG == 'x' and row.DETECTFLAG=='Y':
        texttoadd = ""
        width2= 60-len(row.CHEMICALNM)
        texttoadd = texttoadd + "{0}".format(row.CHEMICALNM)+ "<BOL><ITA><CLR red='255'>{0:>{2}} {1}</CLR></ITA></BOL>".format(row.RRESTXT,row.LABQUALIFI,width2) + '\n'
        myText3 = myText3 + texttoadd
     else:
        texttoadd = ""
        width3= 60-len(row.CHEMICALNM)
        texttoadd = texttoadd + "{0}".format(row.CHEMICALNM)+"<BOL><ITA>{0:>{2}} {1}</ITA></BOL>".format(row.RRESTXT,row.LABQUALIFI,width3)  + '\n'
        myText3 = myText3 + texttoadd

the result of the myText3 variable is this:

I'm trying to use the format script to line up the column on the right side of the box.

Benzene              123245 J
Benzo(A)Pyrene       123 J
Silver               123456789
1,2-Dichlorobenzene  ND U

I have tried to set a variable to the right justify (>{2} instead of an actual number) but it doesn't seem to work.  My hunch is that the formatting tags are adding lengths to the string that need to be accounted for in the format function, but that's purely a shot in the dark right now.

Any tips on how to write the format function to pad the proper length?


texts:
I have the following project and I'm figuring out how to best implement it.  I would like the project to be as highly automated as possible.  The suggestion was made that I make and access table, link it to ArcGIS and then share the directory/map created through sharepoint.

Objective: Put together a list of facilities (company owned, company leased, government leased) on a map.

Criteria:

Develop a system that can be easily updated by non-GIS users.
Background data can be stored in an easily updated program (such as Microsoft Access) and dynamically linked to the mapping function.
Data sets should be accessible to various people in the company (via something like sharepoint).
Mapping function should allow for layers of data

My main question is the use of access and share point the best way of going about the project?

Also lets say I wanted to link my map to an online database (Orcale based) that's protected by a login screen is there any way to do that?


texts:
I want to create a viewshed around a LineFeature. The problem is, I have no observer points.
I would like to create them with a gap of 10m around the line feature. The LineFeature should be 1m over the terrain model and the observer points should have a height of 1.5m.
The height of the woods around the Linesegment are estimated with 20m, which hinders the view of the objects that lay behind it.
Ultimately I would like to calculate the total surface around the lineFeature from where at least a part of the Line is visible.
Anybody knows a away how I should approach that task?


texts: Combining two reclassified rasters using matrix of values with ArcGIS?

texts:
My purpose is to generate an accessibility raster for hiking trails using two reclassified rasters based on slope and land cover data. Both slope and land cover reclassified rasters have values from 1 to 5, so there will be 25 combinations in the matrix.

I need to assign an accessibility code for each combination of values (eg: value 1 for slopes + value 5 for land cover = 2 accessibility code, so that ArcGIS will generate pixel values of 2 everywhere the combination 1 for slopes and 5 for land cover is present).

In other words, I need a tool to automatically generate a code everywhere a condition is met (value X+value Y=Z value in the outcome).

Is there a way to achieve this?


texts: PostgreSQL 'forgot' PostGIS functionality

texts:
System details:

PostgreSQL 9.3.5 64bit on Windows 7 Ultimate;
PostGIS 2.1.3

Two days ago I re-organised my schema structure  - effectively moving a series of tables of the form vm[X]YYYYMM.[tablename] (X∈[admin, prop,add,reftab, plan]) to vmYYYYMM.[tablename]. There are no conflicting table names, so the reorganisation went off without a hitch.

This was done primarily to enable me to set schema when doing intra-period operations - tidying data, adding indexes, setting foreign keys etc - and to enable some more flexible scripting approaches to data import/indexing.

The change - creating the new schemas (schemae? TNIRN) and moving tables en-masse to their new homes - was the only significant change to the database.

Today I tried to perform a standard query of the type that I have performed a thousand times in the old schema-prefixed setup:

select a.lga,
       b.group3,
       sum(ST_Area(a.g1::geography))/10000 as area_ha
  from vmplan201506.plan_zone a
  join vmplan201506.plan_codelist b
    on a.zone_code = b.zone_code
 group by a.lga, b.group3
 order by a.lga, b.group3;

I ran that same query just now, within the vm201506 schema, as follows:

set schema 'vm201506';
select a.lga,
       b.group3,
       sum(ST_Area(a.g1::geography))/10000 as area_ha
  from plan_zone a
  join plan_codelist b
    on a.zone_code = b.zone_code
 group by a.lga, b.group3
 order by a.lga, b.group3;

Result: "ERROR:  type "geography" does not exist"

To test, I dropped the geography cast, and re-ran the query... this time 'ERROR: Function ST_Area() does not exist'. News to me...

If I add back the schema prefixes and set schema 'public', everything runs fine (ergo PostGIS exists, and both the geography cast and the area calc happen as expected).

Another oddity: select postgis_full_version() returns an error (even within schema public) - ERROR:  function postgis_full_version() does not exist... but look at the graphic below - it's clearly installed (and I can run spatial queries so long as I am operating in the public schema).

Anyone ever seen this sort of thing before? If so, was it remediable? (Note: restarting PostgreSQL didn't help).

I realise I could simply uninstall and reinstall PostGIS through Control Panel-> Programs and Features, but thought I would pose the question in case it's a symptom of something egregious. There's no indication of a pending disk failure.

UPDATE (no way to spend a Sunday afternoon)...
So I took this opportunity to

pg_dump all data (150GB... less than I thought it would be);
upgrade PostgreSQL to 9.4.4;
uninstall the old PostGIS;
install new 9.4.4-compatible PostGIS (2.1.8);
set up a fresh new database instance;
postgis-ify it;
psql \i the data

Everything went off without a hitch.

EXCEPT ... the problem still exists in a very slightly modified form.

In 9.3.5 I got an error - even with set schema 'public' for select postgis_full_version(). Now I get the output I expect - namely

POSTGIS="2.1.8 r13780" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER

Seems like progress, yeah? Not so fast.

If I set schema to anything else, select postgis_full_version() chucks an error- good old 'function does not exist'.

Talk about maddening.


texts: How to supply a weight matrix to focal_hpc in R

texts:
I'm using the function focal_hpc from spatial.tools in R to compute for each cell of a raster the average value (mean function) around windows of a specific size (3 cells x 3 cells).
I used to perform this analysis using the similar function focal, which is however slower than focal_hpc when using large raster files.
The advantage of focal is that the user can specify a matrix of weights. However, focal_hpc does not seem to provide this option.
I have been trying to use a custom function, such as:
f_mean = function(x,weight_pix){
  x[is.na(x)]<-0
  weighted.mean(x,weight_pix)
}

However, I keep getting this error, which I do not really understand:
Error in weight_pix * x : non-conformable arrays

Do you have any suggestion?


texts: Changing raster value based on point feature

texts:
I am working on some landslide modeling.

Within ArcMap I have a raster layer, representing the area of the county I am working in, and a point layer that shows landslide locations. I reclassified the county raster so all cells have a value of "0".

Is there a native GIS tool or script that allows me to change the value of the county raster based on landslide point location?

Whenever a landslide point (or multiple points) are located within one sell, the value should change from "0" to "1".

My goal is to have a county raster layer where cells with the presence of landslides are coded "1" and absence "0".

I guess I could also convert the point feature to raster and then use some form of map algebra?! Not sure about the expression needed though.

So far I was only able to find the "extract values to points" tool which accomplishes pretty much the opposite.

Any thoughts or ideas?


texts: Getting buildings multipatch from LiDAR data without building footprints?

texts:
I have a 3D LiDAR data and aerial images.

Is there a way to obtain the buildings multipatch (.shp), given that I don't actually have the buildings' footprints?


texts: Debugging ERROR 000539 from CalculateField in ArcPy?

texts:
I previously successfully used the sample code here:  How to calculate Field using part of filename in ModelBuilder?

# Import standard library modules
import arcpy, os, sys
from arcpy import env

# Allow for file overwrite
arcpy.env.overwriteOutput = True

# Set the workspace directory
env.workspace = r"C:\temp.gdb"

# Get the list of the featureclasses to process
fc_tables = arcpy.ListFeatureClasses()

# Loop through each file and perform the processing
for fc in fc_tables:
    print str("processing " + fc)

    # Define field name and expression
    field = "SID"
    expression = str(fc[:5]) #subsets first 5 characters of fc name

    # Create a new field with a new name
    arcpy.AddField_management(fc,field,"TEXT")

    # Calculate field here
    arcpy.CalculateField_management(fc, field, "expression", "PYTHON")

But now I'm running it again, with no changes, even with the same datasets I used before, and it no longer works.  I cannot see how or why it errors out.

Error message:

  Runtime error : ERROR 000539: Error running expression: expression : name 'expression' is not defined Failed to execute (CalculateField).


texts: Database Population in wrong Authentication Type

texts:
I've populated a database while connected to it as an editor with an Operating System Authentication Type.  I should have populated it while connected to the Database as the GIS administrator with a Database Authentication Type.

Consequently, After I populated it, I connected to the database as the GIS Administrator and the container is empty.

Seems like I need to push my version from one to the other and I am not sure if there is a way to do that and what that would be?

In other words: I have to push the editing Default version into the GIS Default Version.


texts:
We are using the latest version of degrib (2.02) http://www.nws.noaa.gov/mdl/degrib/tutorial.php with some grib files but I am having a hard time getting the "Input on stdin" feature to work. See "-in" #3: http://www.nws.noaa.gov/mdl/degrib/txtview.php?file=degrib.txt&dir=base

I am unable to get this to work doing something like cat gribfile.bin | degrib -opt1 -opt2 -etc From what I could tell, that should work, but it complains of "-opt1" not being a valid file (still trying to read the first arg). Has anyone gotten this to work?

The end goal is to be able to do a database probe (-DP) from streams, not sure if that will be possible though.

Thanks,

-Seth


texts: Image Resampling in Google Earth Engine?

texts:
I've wrote below code for landsat and sentinel mosaicing, however I don't know how can I resample sentinel-2 pixel size to landsat-8 data? Could anyone guide me please?

code link: https://code.earthengine.google.com/e65292e8cda25e846d2a9889a91a3b45

Map.centerObject(table)

var sen2mosaic = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2018-01-01','2019-01-01')
.filterBounds(table)
.map(function(img){
  var bands = img.select('B[3-4]','B8')
  .multiply(0.0001).clip(table)
   return bands
  .copyProperties(img,['system:time_start','system:time_end'])

})
.min();

Map.addLayer(sen2mosaic,[],'sen2mosaic',false)

var lan8mosaic = ee.ImageCollection('LANDSAT/LC08/C01/T1')
.filterDate('2018-01-01','2019-01-01')
.filterBounds(table)
.map(function(img){
  var bands = img.select('B[3-5]')
  .multiply(0.0001).clip(table)
  return bands
  .copyProperties(img,['system:time_start','system:time_end'])

})
.min();

Map.addLayer(lan8mosaic,[],'lan8mosaic',false)


texts: Writing Custom Map Tools using RectangleMapTool

texts:
I am creating a plugin with a preview window where the user must select a rectangle to define the study zone.
I see in this https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/canvas.html the Writing Custom Map Tools example but I do not know how to call it and give the possibility to the user to draw the rectangle.
I put this code in the .ui dialog python script:
class RectangleMapTool(QgsMapToolEmitPoint):
    def __init__(self, canvas):
        self.canvas = canvas
        QgsMapToolEmitPoint.__init__(self, self.canvas)
        self.rubberBand = QgsRubberBand(self.canvas, True)
        self.rubberBand.setColor(Qt.red)
        self.rubberBand.setWidth(1)
        self.reset()

    def reset(self):
        self.startPoint = self.endPoint = None
        self.isEmittingPoint = False
        self.rubberBand.reset(True)

    def canvasPressEvent(self, e):
        self.startPoint = self.toMapCoordinates(e.pos())
        self.endPoint = self.startPoint
        self.isEmittingPoint = True
        self.showRect(self.startPoint, self.endPoint)

    def canvasReleaseEvent(self, e):
        self.isEmittingPoint = False
        r = self.rectangle()
        if r is not None:
          print("Rectangle:", r.xMinimum(),
                r.yMinimum(), r.xMaximum(), r.yMaximum()
               )

    def canvasMoveEvent(self, e):
        if not self.isEmittingPoint:
          return

        self.endPoint = self.toMapCoordinates(e.pos())
        self.showRect(self.startPoint, self.endPoint)

    def showRect(self, startPoint, endPoint):
        self.rubberBand.reset(QGis.Polygon)
        if startPoint.x() == endPoint.x() or startPoint.y() == endPoint.y():
          return

        point1 = QgsPoint(startPoint.x(), startPoint.y())
        point2 = QgsPoint(startPoint.x(), endPoint.y())
        point3 = QgsPoint(endPoint.x(), endPoint.y())
        point4 = QgsPoint(endPoint.x(), startPoint.y())

        self.rubberBand.addPoint(point1, False)
        self.rubberBand.addPoint(point2, False)
        self.rubberBand.addPoint(point3, False)
        self.rubberBand.addPoint(point4, True)    # true to update canvas
        self.rubberBand.show()

    def rectangle(self):
        if self.startPoint is None or self.endPoint is None:
          return None
        elif (self.startPoint.x() == self.endPoint.x() or \
              self.startPoint.y() == self.endPoint.y()):
          return None

          return QgsRectangle(self.startPoint, self.endPoint)

    def deactivate(self):
        QgsMapTool.deactivate(self)
        self.deactivated.emit()


texts: How to Warp vector file based on control points in ArcMap?

texts:
I have a very comprehensive polyline file that, when compared with imagery, is obviously offset by a few hundred meters, but the offset varies from one area to another. It would be easy for me to create a set of control points to correct this, but I can't find the correct tool to do so. I'm using ArcMap, and have just about every ESRI extension known to man, if that helps.

I've entertained the idea of creating a raster and using the raster warp tool, then reconverting back to vector, but I don't want the inevitable loss of fidelity from multiple conversions. If anyone has another way to do this I'm all ears.


texts: GeoServer 2.17.1 WAR running on Tomcat 9 not finding GEOSERVER_DATA_DIR context param

texts:
I have GeoServer 2.17.1 installed as a WAR servlet in Tomcat 9. If the data files are located in
C:\Program Files\Apache Software Foundation\Tomcat 9.0\webapps\geoserver\data
everything works correctly. However, If the data files are located in
C:\Users\geoserveruser\Documents\GeoServerData
and the following context-param is added to
C:\Program Files\Apache Software Foundation\Tomcat 9.0\webapps\geoserver\WEB-INF\web.xml
<context-param>
   <param-name>GEOSERVER_DATA_DIR</param-name>
   <param-value>C:\Users\geoserveruser\Documents\GeoServerData</param-value>
</context-param>

the data directory is not found. In fact a new directory is created in
C:\Program Files\Apache Software Foundation\Tomcat 9.0\webapps\geoserver\data
when Tomcat is restarted. The GeoServer documentation
https://docs.geoserver.org/stable/en/user/datadirectory/setting.html#web-archive
indicates this is the correct way to configure the GeoServer data directory and the context param is provided in web.xml and just needs to be uncommented and edited.
Why is this not working as indicated in the documentation?
UPDATE:
Found these in the Tomcat log file
29 Jun 16:35:17 WARN [geoserver.platform] - Found Servlet context parameter GEOSERVER_DATA_DIR set to C:\Users\geoserveruser\Documents\GeoServerData , but this path does not exist

29 Jun 16:35:17 INFO [geoserver.platform] - Falling back to embedded data directory: C:\Program Files\Apache Software Foundation\Tomcat 9.0\webapps\geoserver\data

This explains why a new data directory is being created.
C:\Users\geoserveruser\Documents\GeoServerData
definitely exists and I verified the current (and all users) have read-write permissions for the directory.
Thoughts?
UPDATE: Attached are screenshots of the path in File Explorer and the security settings for the GeoServerData directory where there is Full Control. (note: this is a local dev environment so path is slightly different from question above, but equivalent).


texts: How to get sum of a subset as opposed to the last value? CARTO

texts:
When choosing a date subset the displayed value of my application is the lastest value as opposed to the sum on that subset.

In the example, 943.974 is the dosage for June 12th. The sum of dosage from June 10th to June 12th is what I need. That is 6306 in this case.

In my data I have:

Geom     Dosage     Date
123123    4132      6/10/2018
123123    1231      6/11/2018
123123     943      6/12/2018
342523    321       6/10/2018
342523    645       6/11/2018
342523    342       6/12/2018
...      ...        ...


texts:
I'd like to overlay a GeoTIFF that is in EPSG:32650, WGS84 on OpenStreetMap (zoom15).
But they don't match scales because resolution of GeoTIFF is 5m/pixel and OpenStreetMap is 4.773/pixel.
How can I get same size of map with GeoTIFF from OpenStreetMap and overlay to right position?
I know that QGIS can do that but I'd like to operate from Python.
And also, I can't resample GeoTiFF because of some reasons.
Here is the code to get OSM.
from staticmap import StaticMap

m = StaticMap(200, 200, url_template='http://a.tile.osm.org/{z}/{x}/{y}.png')

image = m.render(zoom=15,center=[22,140])
image.save('marker.png')


texts: How to remove marker added map with .addTo() function in leaflet?

texts: I am having trouble removing marker which is added with marker.addTo(map). The remove works correctly if added as layer map.addLayer(marker) and then map.removeLayer(marker). But removelayer is not working with .addTo().


texts: Anything similar to TMS or WMTS for WFS?

texts:
I've been working with OpenLayers for a while, but problem got a question with WFS.

Not sure if it's supported by WFS, or OpenLayers, I want to achieve something like this, mostly for cache purpose:

As user pans/zooms the map, OpenLayers keeps sending out WFS requests with BBOX. I want to make coordinates in such requests keep only 3 decimal digits, in a sense, aligning to 0.001 degrees. Do I have to tinkle with OpenLayers, or is it built-in?

EDIT

I'm thinking about tiling for WFS, like TMS and WMTS do. Just we don't render image here, instead we "render" XML files.


texts: Changing data source of ArcSDE layers to replace connection and names in thousands of MXDs?

texts:
I am facing a situation in wich our GIS data provider changed the path to their SDE Data server, and also changed the names of the the layers in their SDE Data.

Need to do changes on thousands of .MXDs:

the path to the new .sde connecting file
the layer names of SDE Data

I am not a programmer, but only a GIS mapper:). It is obvious that a script is needed to be batch-applied on all the thousands .mxds, which now have broken links to all the previous layers on SDE Data of our GIS data provider.

Do you have a code for a script, which can successfully change on all the maps - the path to the old .sde connecting file into a new path to the new .sde connecting file, & also to replace the old layers with the new layers which have different names now?


texts:
I am new to CARTO VL and i want to make a simple dashboard with an animation panel as shown in one CARTO VL guide.

The problem that i get is that when i try to animate it by the atribute of time of the GeoJSON, it doesn't work.

const viz = new carto.Viz(`
     @animation: animation(linear($TIME, time('2019-03-22T00:00:00'), time('2019-03-24T23:00:00')))

I don't know if it has to be with my GeoJSON time format.


texts:
I have several rasters that show the percentage of sediment (mud,sand, gravel) at different locations.  I followed along with the ESRI "Advanced labeling for the stretch renderer" guide but I am encountering a problem where my new labels are being overridden by the default labels.

I want a stretched raster that is comparable across all locations so I am using the Stretch Type "Minimum-Maximum" to change the values to be 0 - 100.  I would like my legend to be labeled at 0, 25, 50, 75 and 100. I used the advanced labeling to set the interval size to 25.  For a second the labels show correctly as I have indicated but after another second they revert back to High : 100 and Low : 0 without any additional labels between. See the screen shots below. What can I do to avoid my labels reverting to the default values?


texts: What programs would allow for the mapping of a geoid in 3d?

texts:
I would like to find a method for plotting an exaggerated geoid (e.g. topography):

I primarily use R for my work but would be interested in other open source software options (e.g. QGIS w/ GRASS). I found this wonderful tool for MATLAB (link) but haven't discovered other options.


texts: Layer Extent NOT Set When Bringing In ArcGIS Server map service?

texts:
I published a mapping service on ArcGIS Server 10.3, set the layer to 'Don't show layer when zoomed' 'in beyond' to 1:577,791, prior to publishing.

When I bring in the mapping service into ArcMap, the layer properties is set to 'Show layer at all scales'.

How do I get the mapping service to automatically show up as 'Don't show layer when zoomed' into my set extent? It can manually be set, but I want to avoid that.


texts: Output feature class is arbitrarily split into sections by Arc GIS

texts:
I’ve attached a picture illustrating my problem. In the picture areas 1, 2, 3 and 4 should all be joined as one polygon. Instead it looks like arc has split the polygon into four parts using a grid, see the lines.

My question is why has Arc GIS 10.3 (Advanced) done this? I know how I could fix this but I would like to prevent this from happening in the first place.

My starting data set was a shapefile containing points, I then ran a model which selected some of those points and buffered around them for 2 meters dissolving overlaps where they occur.

Model summary:

Iterate feature class = this imports my point feature class data

Make feature layer = converts a feature class to a layer

Select layer by attribute = selects all points with  “gridcode” > 0

Select layer by attribute = selects all points with  “p_value” < 0.05

Select layer by attribute = selects all points with  “value” > 1.96

Buffer = creates a 2m buffer around the selection and dissolves overlaps

END OF MODEL


texts:
I used QGIS and its qgis2web plugin to create a web map, using the Openlayers option (from Tom Chadwin and those who helped him!).  I got everything working great, except, the initial view goes to 0 long 0 lat, in the ocean, south of Ghana.  I have searched the web and found many suggestions for setting the initial view of an Openlayers map, and tried them as best I can (I am not a programmer and feel like a monkey copying and pasting code and hoping it works), but with one minor exception after making the edits I get a blank map so have to undo the edits. (The minor exception was I got the initial view to move slightly by changing the extent coordinates in qgis2web.js to coordinates for the extent of my QGIS project, not a solution.)

I have qgis2web set to "fit to layers extent" (if I try canvas extent, the result is worse, my layers are miles above Earth and not lined up with Earth).

I have tried edits to qgis2web.js and other files based upon what these pages say, https://www.qgistutorials.com/en/docs/web_mapping_with_qgis2web.html (step 26), https://stackoverflow.com/questions/3556966/how-to-set-default-view-in-openlayers-without-restricting-bounds, and more.

I am using QGIS 3.4 and according to Plugin Manager the current version of qgis2web, on a Windows machine with Windows 10.  All of my layers are EPSG:3857 - WGS 84 / Pseudo-Mercator - Projected.  I have 3 x, y, z, tile service base maps, also the same CRS.

Here is a the unedited code in qgis2web.js where I think the fix should be applied (but am not sure).

var map = new ol.Map({
controls: ol.control.defaults({attribution:false}).extend([
    expandedAttribution,new measureControl()
]),
target: document.getElementById('map'),
renderer: 'canvas',
overlays: [overlayPopup],
layers: layersList,
view: new ol.View({
    maxZoom: 18, minZoom: 4, projection: new ol.proj.Projection({
    code: 'EPSG:3857',
    extent: [-20037508.342789, -20037508.342789, 20037508.342789, 20037508.342789],
    units: 'm'})
}) });

That is where I tried changing the extents to the extents of my project in QGIS, which seemed to move the initial view a little but not much (I returned those extents to as they were because my guess is those numbers include the entire world, which is ok so long as I don't need to change them to get the map to load on my project in the initial view).

Here are some of the edits I tried, which resulted in a blank map.

    var map = new ol.Map({
controls: ol.control.defaults({attribution:false}).extend([
    expandedAttribution,new measureControl()
]),
target: document.getElementById('map'),
renderer: 'canvas',
overlays: [overlayPopup],
layers: layersList,
view: new ol.View({
    maxZoom: 18, minZoom: 4, projection: new ol.proj.Projection({
    code: 'EPSG:3857',
map.setCenter(new OpenLayers.LonLat(-74.006998, 40.722120)), 10,
    extent: [-20037508.342789, -20037508.342789, 20037508.342789, 20037508.342789],
    units: 'm'})
}) });

I think part of the problem may be that Openlayers has changed some of its terminology.


texts: gdalwarp - correct for a slight spatial misregistration?

texts:
Some Sentinel-2B band images have slight spatial misalignments as I explained on the step forum.
https://forum.step.esa.int/t/sentinel-2-tci-10m-level-2-images-slight-spatial-misregistrations/14235

Furthermore the gml/jp2 band image header goes uncorrected.
I need to extract pixel values from these images as if they had been spatially aligned properly - e.g. by taking into account a 1 pixel vertical shift downwards.

Can I do this using command-line parameters to gdalwarp : https://www.gdal.org/gdalwarp.html - the program I use to extract polygons of pixels?


texts: What's the use of the buffer function in field calculator?

texts:
I'm trying to find out how the buffer command in the filed calculator works.
I tried to create a new field to my point shapefile (integer) with the following:

buffer ($geometry, 100)

but no geometries are added to the canvas and the filed in the table in empty.
The $geometry attribute should refer to a specific field of the shapefile or not?

Someone has successfully used this function?


texts: QGIS 3 won't work on my Windows 10 computer anymore

texts:
After making a rather lenghty SQL query, QGIS 3 started not working properly and after a while it finally ceased to work completely. Trying to start up it back again won't do anything and it will always redirect me to the crash log report.
For some reason this problem is only affecting any instance of QGIS 3 installed in my Windows 10 computer as I tried running QGIS 2.18 which was previously installed on my computer and I haven't came across any problems so far. Running SAGA or GRASS alone also work perfectly.
I've tried erasing any trace of QGIS 3 within my system and trying to reinstall it back again, but the problem persists.
This is one of the crash logs the crash report sent me:
Stack Trace

QThread::currentThread :
QObjectPrivate::~QObjectPrivate :
_CallSettingFrame handlers.asm:50
__FrameHandler3::FrameUnwindToState frame.cpp:1159
__FrameHandler3::FrameUnwindToEmptyState risctrnsctrl.cpp:241
__InternalCxxFrameHandler frame.cpp:321
__CxxFrameHandler3 risctrnsctrl.cpp:267
RtlpExecuteHandlerForUnwind :
RtlUnwindEx :
_C_specific_handler :
RtlpExecuteHandlerForException :
RtlDispatchException :
RtlRaiseException :
RaiseException :
_CxxThrowException throw.cpp:75
qBadAlloc :
QString::reallocData :
QSettings::qt_metacast :
QSettings::endGroup :
QSettings::sync :
QThreadPool::maxThreadCount :
QUrlQuery::setQueryItems :
QXmlStreamStringRef::~QXmlStreamStringRef :
QSettings::QSettings :
QgsSettings::QgsSettings :
QgsPluginRegistry::restoreSessionPlugins :
QgisApp::QgisApp :
main :
BaseThreadInitThunk :
RtlUserThreadStart :


texts:
I'm developing a solution to serve a pan-tropical (and perhaps ultimately global) raster dataset to a cartoDB application. I'd prefer to stick to GDAL/open source tools if possible, but would be glad to accept other suggestions as well. The final output should be in http://domain.com/tiles/{z}/{x}/{y} format

I've built the basic OSM tile server on Ubuntu 14.04, implementing the postGIS >> mapnik >> renderd >> mod_tile tool chain with OSM data.

I'm currently customizing it to process .tif data, but my question is-- is there a better approach given that we're dealing with .tif data? The above is optimized to handle OSM input; is there something else I should be using for rasters? And is there a reason all the blog posts on tile generation are from 2013 or earlier? Should I be using vector tiles instead?


texts: Multiple answers for an attribute in QGIS

texts:
I am configuring the attribute form in a QGIS project to use it in QFIELD, I have a field called "LOCATION" where I need to input the place where a tree has a sickness (can be in the trunk, leaves, roots, braches or twig).
I set up a value map with all this possible values, the question is, how can I combine all of them without having to add all possible combinations into the value map?
I was thinking it could be possible to make it editable in a way that after you input a value from the value map you can still type in manually, but it doesnt work. Also some kind of checkbox with all these values to choose would be nice.
The outout format i was thinking it would be just all values separated by comas. Example:
LOCATION = trunk, roots, brunches


texts: Why are the icons in openlayers edit tool bar missing?

texts:
I have an application using Geoserver 2.2, Openlayers 2.12 and Postgresql. I have been following the tutorial at:  http://www.gistutor.com/openlayers/22-advanced-openlayers-tutorials/47-openlayers-wfs-t-using-a-geoserver-hosted-postgis-layer.html.

I have managed to display my vector data and I can see the save icon but not the modify, pan editing and the others. If i hover over their supposed locations, i can read their (pan, modify etc) titles and even pick and use them but i can't see the icons.

I tried to download the openlayers stable release and placed the theme/default/img folder in the root directory of the Geoserver but i still get no response. am using windows 7.
I think the problem is with the path but i can't still figure it out.

I also tried to change the images to gif and jpeg but they still can't load.


texts: I have two shapefiles with some of the parcels (surveying) overlapping each other.  How do you get those overlapping parcels and prevent them from doubling up?  Is there a way you can have it dissolved or union etc... so that second parcel with the same "tax id number" from being included in the merge?


texts: Generate pattern inside polygon in Leaflet

texts:
I'm new to Leaflet and mapping. I'm currently looking for solutions to how I can autogenerate a line pattern inside a polygon drawn on the map with a pre defined distance between each line.

Does anyone have any suggestion to how I can achieve this?


texts: Writing SQL WHERE Clause using Multiple Choice values in FME?

texts:
I published a parameter which is of Choice Multiple type. I wish to
use values from this parameter for WHERE Clause when Reading Inputed
Data in the begin of Workbench.

So there is no problem if I have type of published parameter just
Choice, then I use following:

WHERE Clause Tile = '$(Name)'

but I do not know how to write SQL WHERE Clause statement if I have
Multiple Choice and want to read all data for Name of Tiles selected
in Multiple Choice, because SQL form above does not work ...

Example:

In multiple choice i select the names: Anna, John, Tim.
Then, I want that data for tiles that have those names are read.


texts: Is it possible to add summary statistics to a legend in QGIS print composer?

texts:
I would like to be able to add summary statistics to a legend in print composer e.g.

GB Agricultural Land Classification:
- Grade 1 (Area: xx hectares, xx% of total)
- Grade 2 (Area: xx hectares, xx% of total)
...

The vector layer (GB Agricultural Land Classification) is styled by category (Classification - Grade 1, Grade 2 etc).

Is this possible?


texts: Obtaining vegetation data for a coordinate location (lat/long) using ArcMap

texts:
I have location coordinates (lat/long) for a province in Canada and I want to map the type of vegetation to that location. I am trying to extract the vegetation data from land cover files which are in .tif format.

To do this should it be a point query, or does creating a buffer around the coordinate and then taking some weighted average value makes more sense?


texts: Importing GDB WKB geometry error (OGR Geometry Factory err:6)

texts:
I have prepared all my layers in File GDB (converted the projection to EPSG3857) and uploaded the entire File GDB to PostGIS.

Somehow, it returned errors described below

Yet, it is uploaded to my PostGIS server and seemed to work just fine.

The error messages though, are those just an eyesore? Or a critical issue?


texts:
I spent the whole day trying to install GDAL binding for Python. And now it seems to work. At least this import statement now does not produce any error messages:

from osgeo import gdal

now I want to use this creature, which is supposed for creating tiles.I tried two different ways:

import gdal_retile

and

from osgeo import gdal_retile

But they do not work.

How can I use this functionality?


texts: Auto-incrementing attribute fields using ArcPy cursors and conditional statements?

texts:
I am trying to create unique IDs for an attribute named "FACILITYID" (str) and have got a ton of help from related posts and online python tutorials. I selected the rows I wished to create unique ID's for, and used this autoIncrement() function in the field calculator which works fine.

def autoIncrement(pStart=1,pInterval=1):
    global rec
    if (rec == 0):
        rec = pStart
        return rec
    else:
        rec = rec + pInterval
        return rec

However, I would like to keep track of the variable rec so that I can make new unique IDs later. After concluding that I must move out of the field calculator, I have attempted to use cursors instead.

To update the maximum ID value, I created a data table in the environment with a single attribute "FACILITYID" that can be updated when the code runs. I have had success with this except for one small detail: I only want to update rows that do not already have a FacilityID. Specifically, I only want the code to work for rows in with the FACILITYID is equal to one space character.

FACILITYID == ' '

I have tried countless if statements, while loops, the getValue() function, and more, but can't seem to get it right. Can anyone help me make this code work for my target attribute fields?

#import modules
import os, sys
import arcpy

#Set variables
rec = 0
workspace = "Database Connections\FE.sde"
arcpy.env.workspace = workspace
feature_name = "FE.SManhole" #change feature class if required
field = "FACILITYID"
table_name = "FE.SWR_FacilityID"

tRows = arcpy.SearchCursor(table_name)
tRow = tRows.next()

#Pull in the current highest value from Facility_ID table
maxID = int(tRow.getValue(field))

# define automatic increment function
def autoIncrement(pStart=maxID,pInterval=1):
    global rec
    if (rec == 0):
        rec = pStart
        return rec
    else:
        rec = rec + pInterval
        return rec

#loop through table to update FacilityID if needed.
cursor = arcpy.UpdateCursor(feature_name) #There is no guarantee of order here
for row in cursor:
#####location where I have been trying if statements.####
    row.setValue(field, str(autoIncrement()))
    cursor.updateRow(row)

#save new rec value to FE.SWR_FacilityID table
maxID = rec
tRows = arcpy.UpdateCursor(table_name)
tRow = tRows.next()
tRow.setValue(field, str(maxID))
tRows.updateRow(tRow)

del tRows, tRow, cursor, row
print 'Highest "FACILITYID" Value: ' + str(maxID)

I am very new to python and I got a majority of this code from previous posts including Auto-incrementing in Python Script with ArcPy cursor?


texts: Why are Garmin provided maps inaccurate?

texts:
Collected tracks on my Garmin GPS (Oregon 450) or in BaseCamp seem to be not in the correct position. There's only a little detail with major highways but it looks as though the route is shifted south, or the highway is shifted north. When viewing the route in Google Earth or Strava, the route follows the road almost perfectly.

Are the Garmin provided maps just really inaccurate or is there something else going on?

Here's a couple pictures so you can see what's going on.  In the Strava image, I'm travelling east-west when I'm north of the highway (which is correct), but on BaseCamp, it shows me travelling east-west when I'm south of the highway.  This has happened with all the tracks I've recorded so far on this route.


texts: With the latest update, is anyone else experiencing issues with copying features from one layer to another, or copying features within attribute tables. I am noticing that the keyboard shortcuts for copy and pasting features seem to be causing issues for me.


texts: How to replace http for https before DOM Construct on load KMZ file with ArcGIS API JS

texts:
I have an application and it is loaded with https, this application loads a map with a KMZ file, and the content of that file has external reference to a PNG image that is loaded with http. This causes the page's https to break.

Look at this example of how this is done.

var kmlUrl = "https://s3.amazonaws.com/geodera-bucket/filesds/5c2e82457cd9e62b0a453ce3-1578081436616-bushfireAlert.kmz";
var kml = new KMLLayer(kmlUrl);
map.addLayer(kml);
kml.on("load", function() {
  domStyle.set("loading", "display", "none");
});

https://codepen.io/krekto/pen/NWPYmRb

I want to replace these http references with https before loading external content.

My application uses AngularJS, JQuery and Dojo, the solution can be with any of these technologies.


texts: Is it possible to add to the map a layer containing results from a sql query

texts:
I'm developping an ArcMap Extension and I need to join an existing feature class with data from a sql query located on another DBMS.

I already set up a Workspace with an OleDbWorkspaceFactory. I'm even able to add to the document any table contained in this database. I'm now trying to add a sql query to the document the same way I do it with tables. The query can be complex and include a couple of table with inner/left joins.

I looked in the ESRI documentation and I can't find the good object to do this.

Thanks


texts: Changing cache size for geocoding via ArcPy?

texts:
I recently posted this query: Arcpy script stalls, CPU usage goes to 0, no error message returned and got as a suggestion to increase the cache size in Arcpy for the geocoding of addresses that I am performing.

The link here describes how to do this for the GUI version of ArcMap.

But, the support page for the GeocodeAddresses_geocoding doesn't mention how to manage the cache via Arcpy, nor can I find other references to it.

How can this be accomplished in Arcpy?


texts: gdal_calc raster calculator syntax for logical operators and other functions

texts:
In documentation for gdal_calc it is stated Command line raster calculator with numpy syntax. Later on there are few examples where in one of them:

gdal_calc.py -A input.tif --outfile=result.tif --calc="A*(A>0)" --NoDataValue=0  - means set values of zero and below to null

Unfortunately there is no example on logical operators like:

--calc="A*(A>0 and A>B)"- means keep A if A bigger zero and bigger B and set the rest to null

Based on Numpy/Scipy logic functions I would expect to write logical operators as:

--calc="A*logical_and(A>0,A>B)"

I tried this and it seems to work but I would like to be assured that is correct.

In the similar way if you want minimum of A and B:

--calc="A*(A<=B)+B*(A>B)"

You can just write:

--calc="minimum(A,B)"

My issue is I can't find any cookbook to make sure I get this right. Is there some good cookbook with advanced examples of what is and is not possible with gdal_calc?


texts:
I want to get the NDVI of MODIS from Google Earth Engine in a country (Ireland), here is my code. It works, but it doesn't show on the map.

var dataset = ee.FeatureCollection('USDOS/LSIB/2013')
            .filterBounds(ee.Geometry.Point(-7.681, 53.43));

 var image = ee.Image().float().paint(dataset, 'iso_num'); ///import the vector data of the country globally

 var visParams = {
 palette: ['f5ff64', 'b5ffb4', 'beeaff', 'ffc0e8', '8e8dff', 'adadad'],
 min: 0.0,
 max: 894.0,
 opacity: 0.8,
 };

 Map.addLayer(image, visParams, 'USDOS/LSIB/2013'); ///show the data on the map

 var image1 = ee.ImageCollection('MODIS/MOD13A1')

          .filter(ee.Filter.date('2018-01-01', '2019-12-31'))

          .filterBounds(dataset);     ////import the MODIS ndvi data, filter the bounds with the dataset above

  var ndvi = image1.select('NDVI');

  var ndviVis = {
  palette: [
  'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
  '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
  '012E01', '011D01', '011301'
   ],
   };
  Map.addLayer(ndvi, ndviVis, 'NDVI');

it runs, but only shows the boundary of Ireland.
How to get the MODIS NDVI separately?


texts: Blank screen when using tiled wmts in Geoserver and Openlayers

texts:
I am currently following the book GeoServer beginners guide (ISBN: 978-1-78829-737-0). I am on chapter 8 and have been using the following code sample from the book:

<!DOCTYPE html>

    <title>WMTS</title>
    <link rel="stylesheet" href="https://openlayers.org/en/v4.2.0/css/ol.css" type="text/css">

When checking the tile layers in geoserver I am able to preview them, however when viewing the output of this code in a web browser, I just get a blank screen. Other resources suggest putting the name of my workspace (packt) in the url:

so it would be:

url: 'http://localhost:8080/geoserver/gwc/Packt/service/wmts',

however this does not work either.
what should the URL be or is the problem somewhere else?


texts: Why iso cluster unsupervised classification detect less clases than desired?

texts:
I have this RGB jpg georeferenced map image:

When I do unsupervised classification with 5 classes

It only gives 4 classes

The original image was generated from CS6 and is georeferenced.


texts: Removing overlapping points using ArcGIS for Desktop?

texts: I am using ArcGIS for Desktop and I have a layer that visually appears to have 50 points. However, the attribute table shows 450 rows, indicating that many points are overlapping. I would like to show only 50 unique points (and have the attribute table reflect that). It is not important which of the overlapping points are saved as I need to do a geospatial analysis and the actual attributes in the table are not important (it's only important that there are 50 of them). How do I create a file with only 50 unique points?


texts: Creating contiguous pixel based cartogram using ArcGIS for Desktop?

texts:
There are several types of cartograms out in the world, but I feel the best one are the continuous cartograms made with separate pixels. I recently made one on paper and then transcribed it to excel by filling in each cell with different colors. This was time consuming and disheartening.  Is there a straightforward way of doing a similar cartogram in ArcGIS?

Examples


texts: Distance from edge of each polygon to nearest point, such that points within each polygon aren't included

texts:
Using ArcGIS 10.5, I have a cluster of polygons that represent districts (together this cluster of polygons creates a nation). I also have thousands of points that represent offices.

I would like to calculate the distance from the edge of each polygon (market) to the nearest point (office) outside of that polygon (market).

My goal is to calculate the distance of the nearest point (office) from each polygon (district) - such that a value is recorded next to each polygon (office) FID. I used the Near tool in toolkit, but it calculates a distance of 0 if there are points (offices) within the polygons (districts). I only want to calculate the distance from the nearest point (office) OUTSIDE of the district (polygon), but not INSIDE of it.

Could you help me with this?

To clarify, every polygon needs a value that is larger than zero. That value should represent the distance from the edge of that polygon to the nearest point OUTSIDE (and never inside) of that polygon.


texts: How to filter OSM data for addr keys with Osmosis

texts:
I'm trying to import the address data of a specific country from an OSM file into a MySQL DB.

I tried to insert the data with osmosis and also to filter them with osmfilter before running OSM but without success.

After a long time of google search and frustration, I decided to post here my problem.

Can someone tell me how I can use the osmosis to import only data with following keys ->

'addr:city', 'addr:housenumber', 'addr:postcode', 'addr:street', 'addr:country', 'addr:housename', 'address',  'addr:place', 'addr:suburb', 'addr:province', 'addr:state', 'addr:number'


texts: Selecting distinct coordinate pair in PostGIS

texts:
I have a PostGIS table containing four columns:

x, y, k_x, k_y

Each row represents a relationship between two points.
Many of the relationships are duplicates albeit in the other direction, like so:

x  |  y  |  k_x  |  k_y
-----------------------
1  |  2  |   3   |   4
-----------------------
3  |  4  |   1   |   2

I'd like to create a table with a unique ID per relation, regardless of direction. So in the example above both those rows would share the same ID.

Selecting by Distinct gives me the unique directed relations, which wont suffice.

I tried this query to at least identify unique rows, but it's not completing within a reasonable time frame so I think I'm missing something there:

SELECT a.*
FROM
    pairs as a
INNER JOIN
    pairs as b
ON
    true
    AND ((a.x != b.k_x) AND (a.y != b.k_y))


texts:
I read about the new version of the digital elevation model at 30 meters (SRTM 30 meters), but I can not find from were I can download it.

some sources:

http://glcfapp.glcf.umd.edu:8080/esdi/index.jsp

http://asterweb.jpl.nasa.gov/data.asp

http://www2.jpl.nasa.gov/srtm/


texts: Merging Two Identical Feature Classes in ArcPy?

texts:
I used the arcGIS tool feature class to feature class to successfull merge two identical feature classes [one with partial data and the other resting on a database connection]. However when I attempt to do this with arcPY I get a ERROR 000258: layer already exists, failed to execute .

The call looks like this

arcpy.FeatureClassToFeatureClass_conversion("roads", C:\fullPathTo.sde, "roads")

Other enviromental varibles have been set accordingly, but I am looking for what the GUI tool does in this situation, an automatic merge, taking the left [local] data and meging it into the .sde.


texts: How to ensure different colors for surrounding polygons?

texts:
Currently in the client I have an array of colors that is looped through for a given set of geometry features.
Features are grouped together on the map  and each group has a color assigned to it.  It is not taking into account adjacent group polygons such that two adjacent groups may have the same color if the array of groups is longer than the color array.

Is there a simple way to assign colors such that no two adjacent groups have the same color?


texts: Programmatically created and added vectorlayer can't be edited

texts:
I'm trying to create a vector-layer (point, shapefile) and open it in my project. It all works fine but once I try to edit that layer later it doesn't work. addFeatures gives me back false and when I try to edit it via the Qgis Tools the adding Points is greyed out.
The odd part is, that I can open the layer I created via the Layer-Menu and it works fine. Also when I type in the code to open it in the console the layer is editable. But somehow in my plugin it doesn't work the same way.
This is the part of the code that should be relevant:

filename = self.dlg.lineEdit.text() # I get the directory from a lineEdit
dirname = os.path.dirname(filename)
path = dirname + '/name.shp'

project = QgsProject.instance()

fields=QgsFields()
fields.append(QgsField('ID', QVariant.Int))

writer = QgsVectorFileWriter( path, 'utf-8', fields, QgsWkbTypes.Point, project.crs(), "ESRI Shapefile")

iface.addVectorLayer(path,"", "ogr")


texts: Why do I get huge memory and CPU spikes setting up the geocode tool

texts:
I am running the geocode tool against a database table which contains values for address, city, state, zip. I am using an address locator that I downloaded locally

http://gis.ats.ucla.edu/

Whenever I bring up the window for the tool and then select the address locator and then subsequently select the fields I get a huge CPU and Memory spike in ArcMap. Almost 1GB of RAM and 100% CPU. This last about 10 seconds for each field and the UI is completely unresponsive during this time. The process eventually works but it looks broken. Using Geocode Address from the right click menu doesn't give these problems, only the tool itself. It's not a huge problem but it makes demonstrations look very shaky.

Thanks for any help.


texts: BIOMOD2: error message: Error in .local(sp, env, ...) : object 'mask.tmp' not found

texts:
I am trying to model a plant species (composed by 298 occurrences) with BIOMOD2 using the consensus landcover as explanatory variables. However, R shows me an error message during the formatting of input data that I do not understand at all:

Error in .local(sp, env, ...) : object 'mask.tmp' not found

I have a feeling that the problem is related to the training area or something similar, because If I run the analysis in other study areas (with the same explanatory variables, but changing the occurrences and the training area), the code works perfectly!

Could someone help me in fixing this error and understanding where it is coming from?

The following code is an example of my used code:

#Loading libraries
library(biomod2)
library(raster)
library(rgdal)

#Loading occurrences
sic.occ.EUR<- as.data.frame(readOGR(dsn = "C:\\Users\\Invader\\Desktop\\minba1", layer = "occ_JAP"))
sic.occ.EUR <- sic.occ.EUR[,c(4,5)]
sic.occ.EUR$Sycios <- 1
names(sic.occ.EUR) <- c("X_WGS84","Y_WGS84","Sycios")

#Loading explanatory variables

Tbio <- raster("C:\\Users\\Invader\\Desktop\\minba1\\ consensus_full_class_4_crop.tif")
Pbio <- raster("C:\\Users\\Invader\\Desktop\\minba1\\ consensus_full_class_7_crop.tif")
myExpl.EUR <- stack(Tbio,Pbio)
myRespName.EUR <- "Sycios"
myExpl.proj.EUR <- stack(Tbio,Pbio)
final.EUR.myResp <- as.numeric(sic.occ.EUR[,myRespName.EUR])
myRespXY.EUR <- sic.occ.EUR[,c("X_WGS84","Y_WGS84")]

#Formating input data
myBioData.EUR <- BIOMOD_FormatingData(resp.var=final.EUR.myResp,
                                       expl.var=myExpl.EUR,
                                       resp.xy=myRespXY.EUR,
                                       resp.name=myRespName.EUR,
                                       PA.nb.rep = 5,
                                       PA.nb.absences = 10000,
                                       PA.strategy ='random',
                                       PA.dist.min = NULL,
                                       PA.dist.max = NULL,
                                       PA.sre.quant = 0.025,
                                       PA.table = NULL,
                                       na.rm = TRUE)


texts: QGIS - automatically export multiple maps varying only the symbology source column

texts:
I have a polygon shapefile where each row represents a sampling grid cell, and each column represents a species, which may be either present or absent in each cell. I need to export one map for each species. The map has no legend, just the grid and a basemap. My solution so far is go to symbology, select a species (column), go to the composer, export image, repeat.

This seems like it should be a fairly common workflow, so I was wondering if there might be some automated way to do it that I'm missing?


texts: Iterating over set of feature classes and spatial join to create only one output?

texts:
I have a feature dataset with 50-some odd feature classes (polygon) that all have fields that I need to spatially join to one feature class (point).

So in the end, I only want one object with all of the necessary fields from all of the feature classes, so the final object will look somewhat like this:

Point Feature Class:

Column 1: Point_FC_ID
Column 2: Layer1_value
Column 3: Layer2_value

...and so forth.

I don't want to have to create 50 some odd outputs for 50 some odd polygon feature classes.  Ideally I would like to just have my base set of inputs and only one output.

How can I do this?

So far... all I can think of is loading all of the feature classes ito the in_memory workspace and then joining all of the fields I want to one output object. I don't know if you can use the same input object as the output object in a spatial join. I would like to use Spatial Join because in side-by-side tests with the other Overlay GP tools, it performs the fastest.

I am using ArcGIS 10.6.1 and arcpy 2.7.


texts: How can I save a layer from WMS in QGIS?

texts:
I´m using a layer from this WMS (http://gis.invemar.org.co/arcgis/rest/services/SIGMA/MANGLARES_COLOMBIA/MapServer). However, I can´t save the layer in my computer, I can´t make geoprocesses or any other action.
How can I download this layer on my computer? Either in .shp or .tiff format (I only manage to download it in .lyr but I don't have ArcGIS and in .kmz but it doesn't open in QGIS).


texts: Replicating how ArcMap assigns colours in colour ramp in QGIS?

texts:
I am new to QGIS and want to replicate the way ArcMap displays values using an HSV colour ramp.

ArcMap seems to assign colours repeating the colour ramp many times throughout the entire list of values

Whereas QGIS assigns colours from beginning to end using the entire colour ramp

The result ends up in very different looking maps between the 2 programs.

Any ideas on how to replicate the ArcMap results in QGIS?

I really want to move away from using ArcMap!


texts: Points to path algorithm - is it possible to keep the attributes in the new file (lines) from the source file (points)

texts:
In my dataset of points, I need first to group points by ID, sort them by time and for each group to draw a line by time order. And for that processing algorithm points to path works perfectly. But my problem is that I lose all the attributes of points after running the algorithm. It would mean to me that at least ID of each group of points became ID for the corresponding line. I have some attributes that are for each group the same and it is really important not to lose them after running the algorithm.

In order to better understand my problem, my dataset looks like as:

ID      TIME        kv           X       Y
1        00:00      a2         x coord   y coord
1        00:05      a2         x coord   y coord
1        01:00      a2         x coord   y coord
1        03:00      a2         x coord   y coord
2        00:00      a6         x coord   y coord
2        17:00      a6         x coord   y coord
3        05:00      ab         x coord   y coord

And my output should be the set of lines with the kv column, so for each user, I should have a line between points by time order:

ID    kv    GEOMETRY
1     a2         geom
2     a6         geom
3     ab         geom

And the real output in after finishing algortithm looks as this:

ID                  GEOMeTRY
new-random-id1       line-geom
new-random-id2       line-geom
new-random-id3       line-geom

Are there any alternative approaches?


texts:
I want to get points at the intersections of the road network. I have a single shapefile. It contains 100k features and some of them continue each other.

As can be seen from 2 pictures above, the selected feature doesn't have any intersections between its neighbors, while it has them from QGIS's point of view. I've tried the solution from Find intersection point in a road network of line shapefile:

Dissolve without attribution the original road network shapefile, which contains 100k features to get one feature containing all roads. It worked perfectly. As a result, I have a new layer with one feature.
Multipart the new layer to Singleparts. After this step I get a new layer, but with the same problem: the one road may have several segments. That's why the Line Intersection doesn't work well.

As you can see, QGIS made an intersection point at the selected road segment, while it must not be done. Any suggestions?


texts: ee.Image(Number) gets interpreted as float

texts:
I never encountered this error before, albeit working with bandmath on GEE for a while.

This is some Function converting DN to radiance. An ee.Image is created from a float value for example LMAXlambda

function radiance(Im){
var LMAXlamda = ee.Image(Im.get('RADIANCE_MAXIMUM_BAND_6'))
var LMINlamda = ee.Image(Im.get('RADIANCE_MINIMUM_BAND_6'))
var QCALMAX_BAND61 = ee.Image(255.0)
var QCALMIN_BAND61 = ee.Image(1.0)
print(LMAXlamda,'1')
return LMAXlamda.subtract(LMINlamda).multiply(Im).subtract(QCALMIN_BAND61)
            .divide(QCALMAX_BAND61.subtract(QCALMIN_BAND61)).add(LMINlamda)

}

following error is Obtained

print(LMAXlambda) --> 15.303 ; but no ee.Image

i hope my Confusion is relateable


texts: How to organise complex styling rules in QGIS?

texts:
For example, my "roads" vector layer has the following attributes:

road type: primary, secondary, local, track, or path
surface: sealed, or unsealed
segment type: road, bridge, tunnel, floodway, culvert, and few more types
relevancy: 1..9 (related to scale / zoom level)

Which yields almost 1000 unique styles but only 10% is present in data (for example, there are no unsealed primary roads)

In general styling rules are simple:

any path should be dashed black line (0.25mm pen width)
any unsealed road should be dashed
primary and secondary roads should be red (0.6mm pen width)
local roads and tracks should be orange (0.4mm pen width)
tracks should be dashed
each segment type has specific rendering:

road: no change
tunnel: keep colour and pen width but change line to dotted and draw portals
and so on ..

So far, I've tried to use "Single Symbol" based on segment type and then use "Expression" for colour, pen width & style. I couldn't figure out how deal with relevancy so I switched to "Rule-based".

"Rule-based" can deal with relevancy, but number of rules is high. I can't figure out how use "Refine selected rule" > "Add categories to rule" to change pen style only and keep the colour.

In both cases I have to copy attribute expressions across the rules (not sure if there is a concept of "global" expression) and it's tedious to keep then in sync.


texts:
We are going to let our users add points, polylines and polygons to a map. Are there any international standards for map symbology?

For example for emergency symbols?


texts: Error: No data value -9.99999999E8 cannot be applied to band #0 of type Short, when importing raster in Google Earth Engine

texts:
I have a question, what does mean the message in the screenshot I uploaded.

Error: No data value -9.99999999E8 cannot be applied to band #0 of type Short<0,255>

Should a raster dataset have some characteristics before importing it in Google Earth Engine, in the Assets tab?
My raster is an 8bit raster dataset, which contains values 0, 1 and 10 for some land cover classes.


texts:
I am using extras/clusterlayer to render clusters for a uniquevalue renderer.. all that works well and good except that the label colors are all white/black/whatever you set them to..

I though I'd just loop through the cluster layer object and change the label colors on a per graphic basis based on the color it got from the renderer.. This does not work. All label change color even if you specifically change one graphic[i]. The loop is working and it seems to swap the colors properly but in the end the last color applied to the label has been propagated through all the graphics..

the little test loop below should alternate white and black labels for each cluster.. It doesn't, the best i can get it to do is change all the labels from black to white and back again.

        clusterLayer.setRenderer(UniRenderer);
        var graphics = clusterLayer.graphics;
        console.log(graphics.length);
        for(var i=0; i<graphics.length;i++)
        {
            if(i%2==0)
            {
                console.log("black");
                //have also tried clusterLayer.graphics[i]._graphicsLayer.clusterLabelColor
                clusterLayer.graphics[i]._layer._clusterLabelColor= "#000";
            }
            else
            {
                console.log("white");
                clusterLayer.graphics[i]._layer._clusterLabelColor= "#FFF";
            }
            //console.log(label);
            //label = "#000";
        }
        console.log(clusterLayer);
        var symbolList = clusterLayer.renderer._symbols;


texts: Does Unicode Common Locale Data Repository (CLDR) library have rules for longitude/latitude coordinates?

texts:
Does the Unicode Common Locale Data Repository (CLDR) library contain rules for formating longitude/latitude coordinates?

I am looking for information how coordinates should be formatted in each language/country.


texts: Interpolate dxf polyline vertex elevation from DTM?

texts:
I am working on several autocad polylines (with several vertices...) that have the exact N,E coords but now I would like to add the elevation data for every polyline's vertices from a dtm.

I thought it was an easy task but I can't reach my goal so far.

I imported dxf polyline in QGIS (I tried also in GRASS).

There it seems that v.drape doesn't work.

In fact, I exported the vector as dxf (with flag -z) but I have all vertices at a same elevation (find from centroid in dtm) and the result is "flat" polylines at different elevation which is no good.

I would like to have each of the polyline's vertices with a different elevation.

Is there a way to do an easy thing like this?

I want to transform a 2d dxf in a 3d dxf by draping it on a dtm.


texts: QGIS Select double feature by attributes

texts:
There is this point shapefile in QGIS (2.12). There is one column with the street name and one with the housenumber.

Is it possible to select all features that have the same street name and housenumber (not one by one, but automatically all that are alike)? If yes, what tool do I have to use or what is the expression if I "select by expression"?


texts:
I have a geotiff and I checked it's metadata using gdalinfo.exe:

Driver: GTiff/GeoTIFF
Files: C:/image_test.tif
Size is 3956, 34151
Coordinate System is `'
GCP Projection =
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
GCP[  0]: Id=1, Info=
          (1,1) -> (119.851499142389,-30.7134418800722,390.999576434493)
GCP[  1]: Id=2, Info=
          (3956,1) -> (119.352774909547,-30.6076899722145,390.999578144401)
GCP[  2]: Id=3, Info=
          (1,34151) -> (118.788734486763,-34.4555300913518,390.99952176027)
GCP[  3]: Id=4, Info=
          (3956,34151) -> (118.265259471947,-34.3453414924152,390.999523188919)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2019:04:07 13:55:02
  TIFFTAG_RESOLUTIONUNIT=3 (pixels/cm)
  TIFFTAG_XRESOLUTION=0.00079999998
  TIFFTAG_YRESOLUTION=0.00079999998
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,34151.0)
Upper Right ( 3956.0,    0.0)
Lower Right ( 3956.0,34151.0)
Center      ( 1978.0,17075.5)
Band 1 Block=3956x1 Type=UInt16, ColorInterp=Gray

I am trying to create an ImageMosaic in geoserver but the corner coordinates are just the pixel ranges and not the actual lat/longs and so the WMS layer is abit crazy.

I've looked at other geotiff's that work and they include within the coordinate system, e.g.:

Origin = (-70.119788596117274,18.999683535544683)
Pixel Size = (0.000256210858208,-0.000256210858208)

Is there a way I can set this and then also apply this to the corner coordinates? Maybe using gdal_transform.exe?


texts: Lambert projection with one or two standard parallels and proj4

texts:
I'm sure there is a gap in my understanding, but I haven't been able to identify where it is, nor how to fill it. My issue is this:

I'm transforming lon lat coordinates into a lambert conformal conic projection.
The lambert projection uses a single standard parallel 46.244 and the central meridian is 17.0. And it's on a spherical ellipsoid, but that's not particularly relevant here. I am using packages that use the proj4 library.
I test the transformation using linux tool cs2cs as follows:

echo "15.220483 48.437298" | cs2cs +proj=lonlat +to +proj=lcc
+lat_1=46.244 +lat_0=46.244 +lon_0=17. +a=6367470 +b=6367470
+ellps=sphere +datum=WGS84

-131346.80      245318.16
I also understand that the one-standard parallel form is the same as the two-standard parallel form if the two standard parallels are the same. It is more convenient for me to give all parallels in my scripts. So I also test with:

echo "15.220483 48.437298" | cs2cs +proj=lonlat +to +proj=lcc
+lat_1=46.244 +lat_2=46.244 +lat_0=46.244 +lon_0=17.
+a=6367470 +b=6367470 +ellps=sphere +datum=WGS84

-131290.66      245281.34

yet I get different results. Any ideas as to why?

And at very least, which should I use?


texts:
Is it possible with the GDAL library to read an image (e.g. GeoTiff) from a memory stream, like stdin instead of a defined location on the local disk?

It would be useful to read data in streaming environments. In my case I would like to use the Hadoop Streaming API, which reads the datastream from HDFS, so passing a file URI to the GDALOpen() function is not possible.


texts: Join Tables related by a composite primary key in qgis

texts:
I have two tables, one has the geometry and other has attributes. The table with geometry has unique number primary key, and 4 columns admin1,admin2,admin3,featureid.

Table two is a csv file has the same columns mentioned above and their combination forms a composite primary key.

How can i join the two tables?


texts: How to convert shapefile to AutoCAD DWG using a 1:1 scale in ArcGIS Desktop?

texts:
I have a potential job that has a requirement of delivering a boundary in AutoCAD. These are the requirements: Digital copy of wetland delineation in AutoCAD shall be provided, with the survey at
1:1 in model space. Survey shall be tied to North American Datum 1983 State Plane
Illinois East Federal Information Processing Standard 1201 (feet).
I do not have the AutoCAD program and am wondering if this deliverable is possible using the convert to CAD tool in ArcGIS Desktop. The reason I am skeptical is I recently tried doing something similar for an engineer for a current job. I am relatively certain my projection matched theirs (same as the one requested above) but they said the scale was 10:1 when they got it, not 1:1, so obviously something was off. They were able to fix it in their CAD program.
Does anyone have a solution or suggestion?
Or is it a projection issue?
Or maybe something that can only be adjusted through AutoCAD?


texts:
I was able to find the right coordinate system online NAD 1983 Alaska

But now I don't know what to do.  I can't find the coordinate system on QGIS.  When I choose the one I think is correct it warps everything.

I have a lot of shp files that were made in the coordinate system NAD 1983 EPSG 102631 which does not show when I use filter and try to set the CRS. When I load the shp files into my project and the CRS says custom everything seems fine. But when I set a CRS for my project it goes crazy.  I want to use these shp files with other data and it seems I need to convert or at least get something sorted out before I can use this data with other data.

It seems I need to custom set the proj strings... I've been reading the manual and I haven't found how to do this yet.


texts: OpenLayers - MousePosition show Y coordinates as wrong for ImageStatic with Pixels Projection

texts:
I have a map with a pixel projection and load a static image in a layer (My image height is 4298).

When I view the coordinates in mouseposition I see that the X is showing correct, but not the Y.

For example: If I move the cursor at the top left corner of my image it shows: 0, 4298 - Which is wrong as it is suppose to show 0, 0. When I move my cursor at the bottom left corner of my image, it shows 0, 0 - Which is also wrong at it is supposed to show 0, 4298. I presume that map click events will also handle the Y wrongly.

How can I correct this?

var imgExtent = [0, 0, 6065, 4298];
var imgProjection = new ol.proj.Projection({
    code: 'my-image',
    units: 'pixels',
    extent: imgExtent
});
var pageSource = new ol.source.ImageStatic({
    projection: imgProjection,
    url: '/url/to/my/image.png',
    imageExtent: imgExtent
});
var pageLayer = new ol.layer.Image({
    source: pageSource
});
var imgView = new ol.View({
    center: ol.extent.getCenter(imgExtent),
    minZoom: 0,
    maxZoom: 10,
    rotation: 0,
    projection: imgProjection,
    extent: imgExtent
});
var imgMap = new ol.Map({
    controls: ol.control.defaults({
        attributionOptions: {
            collapsible: true
        }
    }).extend([
        new ol.control.FullScreen(),
        new ol.control.ZoomSlider(),
        new ol.control.MousePosition({
            coordinateFormat: ol.coordinate.createStringXY(0),
            projection: imgProjection,
            undefinedHTML: '&nbsp;'
        })
    ]),
    layers: [
        pageLayer
    ],
    renderer: 'canvas',
    target: 'panelEditMainMap-body',
    view: imgView
});


texts:
I downloaded data from Natural Earth Data and added it to TileMill (git revision 8044976 on Mac)

I downloaded raster data, borders, cities, airports and timezones. However, the timezones show only a few odd ones. What am I doing wrong? Or is this a bug in the data, TileMill or Mapnik?

[edit] I realised it may have something to do with the polygons. Those odd timezones are "closed" while most other timezones are more ore less straight lines from top to bottom. Maybe mapnik gets confused about what is the inside and outside of the polygon?


texts:
I have a list of addresses and a map including train stations, schools etc., and I would like to calculate the shortest distance from each address to the nearest train station, school etc.

I'm aware that the google maps API can calculate the shortest distance to any specific point, but my problem is complicated by the fact that I also need to find the nearest train station etc. Is there an easy way to do this in QGIS?

I use QGIS
Windows 7 64bit

Thanks!


texts: How to insert a Polygon into MySQL database Table?

texts:
I have shapefiles which I converted to .csv format. I want to insert those into a MySQL Database table. I used ST_PolyFromText according to the MySQL reference manual.

## create table
cur.execute("CREATE TABLE `db1`.`table1` ("
  "`Id` INT NOT NULL AUTO_INCREMENT,"
  "`geometry` POLYGON NOT NULL,"
  "PRIMARY KEY (`Id`));")

sql_statement = "INSERT INTO table1(geometry) VALUES (ST_PolyFromText(%s,0))"

When I try to insert polygon values into the Database. I get the error    3037, 'Invalid GIS data provided to function st_polyfromtext.'

How should I use MySQL Spatial Functions to insert my Polygons ??

My Polygon geometry column looks like this :


texts:
I have a question regarding DEM. I have downloaded SRTM DEM 30 m and I did a mosaic for the data using Mosaic tools in Arcmap to a Mosaic to a new raster. However, I got a DEM raster with values of 0-255 . I need to have a DEM values with meters not with the values of 0-255. I did convert the raster to 32-float but I still get the same value 0-255. and the cell size is 0.0002777.

Any idea how to convert DEM raster with values of 0-255 to values represent elevations?


texts:
I'm trying to make DEM from points in QGIS.

I've read all the other questions regarding the same issue and tried them all. But, none of them worked.

For example, at the end of inverse distance weighted interpolation, it looks like this.

Edit: it worked after adjusting the size of pixels.


texts: How to convert Multipolygon to GeometryCollection?

texts: I get this error (Script Error, File too large: larger than 512kb) in Google Earth Engine and my layer is a Multypolygon, just with 86 polygons. I have another layer with 300 polygons without any problem, and I noticed that it's a GeometryCollection (Type), but I don't Know how it changes. Do you know how can I do that?


texts: mapbox gl - multiple styles and filters per layer

texts:
The intention here is to determine the proper convention for adding multiple styling rules for a single layer in 'layers' list. My impression was such that multiple layer entries would mean a call to draw the whole layer multiple times which should be avoided. However, in cases such as the need to style both fill and stroke (lines) it appears the official recommendation is to add one layer with Fill type, and another layer with Line type. In this instance, will the engine draw the whole layer twice? Does that only happen when the type is repeated? Or rather, is it contingent on ID uniqueness?

This is what I initially imagined handling for an "else rule" and one with a filter would look like:

      {
         "source-layer" : "myLayer",
         "type" : "fill",
         "paint" : {
            "fill-color" : "blue"
         },
         "id" : "myLayer",
         "source" : "test"
      },
      {
         "source-layer" : "myLayer",
         "type": "fill",
         "paint" : {
            "fill-color" : "red"
         },
         "id" : "myLayer",
         "source" : "test",
         "filter" : [ "==", [ "get", "ATTRIBUTE" ], "colorMeRed" ]
      },

Will this draw myLayer, then draw it again, overlapping?