texts: Successfully extracting values from a GRIB2 file...is it possible?

texts:
Let me say first that I'm starting to work with GRIB2 files, I've read the docs, and more or less know how they are structured. My goal is to be able to parse some GRIB2 files from NOAA, more precisely, those from the WaveWatch III model. So please forgive my ignorance if some of what I say sounds too naive :D

I'm not attached to any technology, although I feel comfortable with Java and Python. To be more precise, you can find the GRIB2 file I'm using as a sample at this link.

I'm also using the ToolsUI from Unidata for a more "user friendly" inspection of the GRIB2 file. As you can see, there are some variables  in the provided GRIB2, but I'm only interested in 2 of them:

Wind_direction_from_which_blowing_surface
Direction_of_swell_waves_ordered_sequence_of_data

So far, I've tried this code, which is based on this question also, and I'm able to open the GRIB2 file:

Grib2Input grib2Input;
String pathfile = "glw.grl.WDIR.grb2";

//Create RandomAccessFile
ucar.unidata.io.RandomAccessFile gribfile = null;
    try {
    gribfile = new ucar.unidata.io.RandomAccessFile(pathfile, "rw");
}
catch (IOException e) {
    e.printStackTrace();
}

//Create grib2Input file
grib2Input = new Grib2Input(gribfile);
Grib2Data gd = new Grib2Data(gribfile);
try {
    grib2Input.scan(false,  false);
    List<Grib2Record> records = grib2Input.getRecords();
    float[] data = null;
    for (int i = 0; i < records.size(); i++){
        //First get record information
    System.out.println("record " + i +  "-> Param. Categ.: " +   records.get(i).getPDS().getPdsVars().getParameterCategory());
    System.out.println("record " + i +  "-> Param. Number: " + records.get(i).getPDS().getPdsVars().getParameterNumber());
    try {
    data = gd.getData(records.get(i).getGdsOffset(),records.get(i).getPdsOffset(), records.get(i).getId().getRefTime());
    }
catch (Exception e) {
    e.printStackTrace();
}
//Second print some data
for (int j = 0; j < 50000; j+=5000){
    System.out.println("data[" + j + "] -> " + data[j]);
}

System.out.println("--");
}

}
            catch (IOException e) {
              e.printStackTrace();
            }

Now the output, for the first 2 records, is this:

record 0-> Param. Categ.: 2
record 0-> Param. Number: 0
data[0] -> NaN
data[5000] -> 64.34
data[10000] -> NaN
data[15000] -> NaN
data[20000] -> 38.13
data[25000] -> NaN
data[30000] -> NaN
data[35000] -> NaN
data[40000] -> 66.55
data[45000] -> 107.06
--
record 1-> Param. Categ.: 0
record 1-> Param. Number: 7
data[0] -> NaN
data[5000] -> NaN
data[10000] -> NaN
data[15000] -> NaN
data[20000] -> NaN
data[25000] -> NaN
data[30000] -> NaN
data[35000] -> NaN
data[40000] -> NaN
data[45000] -> 328.5

Well, here are my question:

1.- Is there a better way to do this with NetCDF-Java library? It looks to me that for bigger size GRIBS this method is going to take ages. Besides, I should remove all NaN values from the data array for each record

2.- Should I go for more GRIB2 specific libraries like pygrib or pynio? If so...any working examples on parsing GRIB2 files with those?

3.- I also tried using wgrib2 command line tool...like...generating a CSV file with that tool and then parse the CSV file. Thing is that the generated CSV file is not very comprehensive either

Thanks!!!
Alejandro

EDIT: You can find an answer to this (using Python, tough) here. I must also say that I was able to do this too using Netcdf4Java, although with GRIB-API from ECMWF is way simpler.


texts:
I am trying to read from a local instance of the OpenStreetMap database all roads intersecting a boundary coming from another database.

I actually get a result with the following piece of code:

#postgis_engine is the SQLAlchemy postgis engine
#boundary_geom is the polygon representing the boundary
roads = geopandas.GeoDataFrame.from_postgis(text('SELECT name, way FROM planet_osm_line WHERE highway IS NOT NULL AND ST_Intersects(way, ST_GeomFromText(:b, 3857)) = true'),
            postgis_engine, 'way', {'init': 'epsg:3857'}, params={'b': boundary_geom.wkt})

As you can see I convert the boundary to its wkt representation and use it as a parameter for the SQL query.

I am wondering if, in a such way, I can break some text parameter length limitation of the underlying database access technology (i.e. SqlAlchemy).

Is there another way to pass a geometry parameter to the GeoDataFrame from_postgis method?


texts: Gdal_calc.py problems from command line?

texts:
I hope someone can help me with this issue with gdal_calc.py. When I call gdal_calc.py from the command line, I get the following error (after displaying the error message I will describe my current environment):

C:\Users\xx> gdal_calc.py
RuntimeError: module compiled against API version 9 but this version of numpy is
7
Traceback (most recent call last):
File "C:\Anaconda2\Scripts\gdal_calc.py", line 46, in <module>
from osgeo import gdalnumeric
File "C:\Python27\ArcGIS10.3\lib\site-packages\osgeo\gdalnumeric.py", line 1,
in <module>
from gdal_array import *
File "C:\Python27\ArcGIS10.3\lib\site-packages\osgeo\gdal_array.py", line 26,
in <module>
_gdal_array = swig_import_helper()
File "C:\Python27\ArcGIS10.3\lib\site-packages\osgeo\gdal_array.py", line 22,
in swig_import_helper
_mod = imp.load_module('_gdal_array', fp, pathname, description)
ImportError: numpy.core.multiarray failed to import

I am running Windows 7, installed Anaconda and the command line calls Anaconda's python version, and I have installed ArcGIS (my employer requires to script in ArcGIS from time to time).

When I call gdal from the command line, it runs without problems. Also, whenever I run gdal_polygonize.py there is no issue. So far, I have only run into problems when running gdal_calc.py as it seems to access the osgeo installation in ArcGIS and not the osgeo installation in Anaconda.

The osgeo path in Anancondo is:
C:\Anaconda2\Lib\site-packages\osgeo

The osgeo path in ArcGIS is:
 C:\Python27\ArcGIS10.3\Lib\site-packages\osgeo

Also, to confirm the Python installation accessed:

Python 2.7.11 |Anaconda 4.0.0 (32-bit)| (default, Mar  4 2016, 15:18:41) [MSC    v.
1500 32 bit (Intel)] on win32
Type "help", "copyright", "credits" or "license" for more information.
Anaconda is brought to you by Continuum Analytics.


texts: Replacing rgb value (0,0,0) with (nodata,nodata,nodata) in GDAL?

texts:
I have a set of rasters in which I want to replace any rgb value of (0,0,0) with (nodata,nodata,nodata)

I've tried using

gdal_translate -a_nodata "0,0,0" in.tif out.tif

But this seems to convert all 0s to nodata, for example (0,10,10) becomes (nodata,10,10).  I want these to stay the same and only convert (0,0,0) to (nodata,nodata,nodata).

I'm looking for any solutions using gdal tools like gdal_translate or gdalwarp that I can add into a larger shell script.


texts: ERROR 000358 When trying to select by attribute using "IS NULL" query

texts:
I am trying to select a layer by attribute (field is null) using a stand alone python script.

    #Select layer by attribute to find any null values
    arcpy.AddFieldDelimiters(mayIntLyr, "DEVICE_ID")
    select = '"DEVICE_ID" IS NULL'
    arcpy.SelectLayerByAttribute_management(mayIntLyr, "NEW_SELECTION", select)

It keeps returning error 000358, indicating the expression is invalid.
Not sure what is wrong with my expression, and I tried adding field delimiters as well and it still doesn't work.

I am using ArcGIS 10.4, and have attempted to incorporate an answer from a similar thread to my code (shown below) to allow for null values, blank values, and spaces, but it still doesn't work.

     select = "DEVICE_ID IS NULL and DEVICE_ID != '' and DEVICE_ID != ' '"


texts: How to subtract lines with thickness from polygon

texts:
I have a polygon and some lines crossing it. I want to subtract the lines with a certain "thickness" from the polygon, so it splits up into various polygons which do not touch each other anymore. You can see the two situations in this picture:

Is there a way to do this?


texts: Javascript App works fine...and then it doesn't

texts:
This question probably won't fit into the Stack Exchange guidelines because of a number of factors, so bare with me.

I have built a GIS web app at http://matthewiannowlin.com/planning/thesis. I used the JavaScript API. At a certain point, a problem started showing up where about 4 out of 5 times, the Esri stuff won't load. So the map and all the widgets. If I refresh enough times I'll finally get the whole page to load. Then it works perfectly.

When I look at console log, I see this error on the times Esri widgets don't load:

[Error] TypeError: null is not an object (evaluating 'e.renderer.toLowerCase')
ia (init.js, line 33)
fa (init.js, line 33)
g (init.js, line 35)
(anonymous function) (init.js, line 40)

I don't even know what to ask, I'm afraid. Besides...what the hell? I feel like this became a problem as the code become more complex and/or as I called more widgets.

Here are some random thoughts: does this have something to do with the timing of how things are load? Does this somehow relate to dojo/ready or dojo/DomReady? Or parsing.


texts:
I was told to implement a web service which, given 2 inputs (street name, cross street name) the output is the crossing geographical coordinates.

What is the best way i should go?

Is there a standard for Geographical Web Services that queries a database?

Are there same similar examples? (I haven't found anyone)

We use PHP, and a Postgres database with Postgis.

I have developed a solution whith nusoap and xml but I really don't know if
this is a good solution.


texts: Combine Elevation Raster and MultilineString to create MultilineStringZ in PostGIS

texts:
I want to show how steep a given stream is, with PostGIS and QGIS.

My thought process for solving this problem of steepness is as follows:
1). Given Raster and Multilinestring data, combine them into a MultilineStringZ column
2). Perform some elementary calculus on the new MultilineStringZ column to find steeper sections. Colorize them in QGIS if possible to examine the results.

Elevation Data:

Stream Data as MultilineString:

I need a better understanding on how a MultilineString and a raster layer are combined to create a 3D stream, or alternatives to completing the task of encoding "steepness" into line strings.


texts:
I was wondering if you could help me with my problem. I want to execute a Map Algebra expression with the QGIS Raster Calculator. It would be easier to use the GDAL Raster Calculator, but since the rasters cells of my two rasters are not congruent, the GDAL tool ist not working. The reason for that is that I create a raster from a 3d polygon with TIN Interpolation and firstly the cell size is 0.999 instead of 1 and secondly I have no chance to snap my new raster to the old one. Maybe someone has a solution for that?

Here is my Python code and I cannot find the mistake :/.

ras_lyr1 =  QgsRasterLayer(input1)
ras_lyr2 =  QgsRasterLayer(input2)

entries = []
entry = QgsRasterCalculatorEntry()
entry.ref = 'dgm@1'
entry.raster = ras_lyr1
entry.bandNumber = 1
entries.append( entry )

entry2 = QgsRasterCalculatorEntry()
entry2.ref = 'flaeche3d@1'
entry2.raster = ras_lyr2
entry2.bandNumber = 1
entries.append( entry2 )

crs = QgsCoordinateReferenceSystem('EPSG:31287')

calc = QgsRasterCalculator('(dgm@1/flaeche3d@1)<0',out11,'GTiff',ras_lyr1.extent(),crs,
ras_lyr1.width(),ras_lyr1.height(),entries)
calc.processCalculation()

What I am doing wrong here?
Thanks a lot for your help!
Niki


texts: Edit or update ArcGIS Online Web Map data using the REST API

texts: Web Map data is described on the ArcGIS REST API site under the Common data types section.  I am curious whether I can request the current JSON configuration as shown in the examples (using /data in the url), modify it, and POST it back?  As an example, I could request the web map data, add or modify the bookmarks list (by editing the JSON web map data, see example here), and submit the update to modify the web map item in ArcGIS Online.  Does a REST operation or alternative exist for this?


texts:
I m using ArcGis 10.4.1 and Network Analyst extension. When i load locations in order to be used as Stops i lose information. ie in the following image is shown the attribute table of the location to be used

When i open Stops attribute table the following information is presented without any field from locations attribute table

How can i get the information of my locations?


texts: Creating buffer around particular raster cells using ArcGIS for Desktop?

texts: I would like to draw 100m buffers around particular raster cells (artificial surfaces) in my LandCover raster layer in ArcMap10.2. The Buffer tool only works with polygon and the Expand tool always crashes the whole ArcMap, because its too much data, I think. That´s probably because the expand tool wants to create buffers around every single 100x100m raster cell, if I understand it right. However many of the cells, displaying artificial surfaces, are next to each other, so I only want a Buffer around these cell groups and not around every single cell. Does anyone knows how I can do this?


texts: gdalwarp "ERROR 1: Couldn't determine X spacing"

texts:
Can someone tell what is causing this error in gdal warp "ERROR 1: Couldn't determine X spacing"?

I have a set of .xyz files in a folder that I want to convert to .tif files. The variable file_list is a list of filenames.

# warp coordinates of each tile
for f in file_list, desc='warp':
    bashCommand = "gdalwarp -t_srs EPSG:4647 -overwrite " + f + ' ' + f[:-4] + '.tif'
    process = subprocess.Popen(bashCommand.split(), stdout=subprocess.PIPE)
    output, error = process.communicate()

For most of the files this command works fine, but for one it gives me the error.


texts: Efficient way to find nearest feature between huge postgres tables

texts:
I have two tables (g1 and g2), both with around two million features. I need to know for each feature in g1, which is the closest feature from g2, with a condition. Since they are both so huge, I also restricted the distance to 150 meters. I'm using pgadmin, and I'm using this query:

SELECT DISTINCT ON (g1.gid) g1.gid As gref_gid, g2.gid as gnn_gid, g2.code_mun, g1.codigo_mun, g2.text, g2.via INTO igcvia1
FROM u_nomen_dom As g1, u_nomen_via As g2
WHERE g1.codigo_mun = g2.codigo_mun AND ST_DWithin(g1.geom, g2.geom, 150)
ORDER BY g1.gid, ST_Distance(g1.geom,g2.geom)

The thing is that it is taking more than 24 hours to finish the query, and I don't know how much longer it's gonna take. I am wondering, is this the most efficient query? Technically, this is measuring all distances and then just show the closest one. Is there a way it can just find the closest one and then stop measuring?
It would be great if you guys could help me find a more efficient way of doing this. I also tried QGIS but since these tables are so big, it just shuts down.


texts:
The task is to scale a polygon in QGIS.
The procedure is performed tracing the polygon and then apply a scaling factor ( the factor used is " Distance on the ground" / " Distance on map " ) .

Both ArcGIS and AutoCAD there is the "scale" tool.

Is there any way to do Qgis or plugin?


texts:
I have two line layers . One is the stream network (dissolved) and other is the point layer converted to line layer by converting each point into short line segments. I am using "Split by lines" tool available by processing toolbox which splits the input line layer at the points where other line layer intersects .The resulting split layer has more splits and at locations other than the intersecting line layer points.

Is there any alternative for the tool or any update to the used approach to achieve the same.


texts: Creating new point geometry MapInfo File (TAB) in PyQGIS

texts:
I am trying to create empty tab file with different approach's.
1st approach
fn = r'E:\data\spilce.tab'
shpfields = QgsFields()
id = QgsField('id', QVariant.Int)
shpfields.append(id)
gtp = QgsWkbTypes.Point
cr = QgsCoordinateReferenceSystem('EPSG:4326')
QgsVectorFileWriter(fn, 'utf-8', shpfields, gtp, cr, 'MapInfo File')

2nd approach
fn = r'E:\data\spilce.tab'
cr = QgsCoordinateReferenceSystem('EPSG:4326')
ipointlayer = QgsVectorLayer("Point?crs='EPSG:4326'", "temporary_points", "memory")
ipointlayer.dataProvider().addAttributes([QgsField("Id", QVariant.String)])
ipointlayer.updateFields()
QgsVectorFileWriter.writeAsVectorFormat(ipointlayer, fn, 'utf-8', cr, 'MapInfo File')

Above codes creating non-geometry tab files.
How to create Point Geometry TAB file?


texts: How to center a MapBox GL JS map at the point coordinates from a point edited with MapBox-GL-Draw plugin

texts:
I'm using MapBox GL JS and the mapbox-gl-draw plugin in a simple MapBox map .... here you're the code ...

    <meta charset='utf-8' />
    <title></title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

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

    <link rel='stylesheet' href='https://api.mapbox.com/mapbox-gl-js/plugins/mapbox-gl-draw/v0.18.0/mapbox-gl-draw.css' type='text/css'/>

When the user click / edit a point on the map I'd like to set the map center at the lat / lon coordinates of the edited point.

Suggestions / examples?


texts: SAGA has a feature that determines flats. Does anybody know how they define a flat? If I start with a raster, each pixel is flat - so it has to be more sophisticated than that. I would like to be adjust the thresholds to consider it a flat if two pixels are within 15m vertical of each other. Would that be possible? How?


texts:
I am running a shell script for batch processing in GRASS 7 (calling it with sh [...sh] in the command line window), which creates a series of rasters through classifications, raster calculator and regression analysis.
What I would like to do now is extend that script, so that it calculates the surface roughness of the resulting raster. r.roughness only works so far in GRASS 6, however. I have GRASS 6 installed as well.

Is it possible to call GRASS 6 and its extensions with the GRASS 7 command line without starting GRASS 6 or do I have to write a script that calls both programs from outside of either GRASS environments?


texts:
I have a task where I need find the density of building footprints. I have a vector layer consisting of building footprints in polygons.

To find the density I was thinking about drawing a set of squares and then make QGIS calculate how many polygons are within each square and then color them according to that number, for example from black to white, where black would be the highest density, and white the lowest. If a polygon occurs in two or more squares it should be included in the one where most of the polygon is present.

I hope some of you have any knowledge regarding the task, and perhaps a step by step idea on how to do it, as I have very limited technical skills in QGIS.


texts: Performing Weighted Overlay using vector data with no inherent values in ArcGIS Desktop?

texts:
I have been asked to do a weighted suitability analysis using all vector based polygon feature layers (8 different layers). Each layer has been ranked from 1 to 8 in order of importance.

My problem is that I need to convert each layer into rasters in order to do a weighted overlay, but many of the layers only contain boundary polygons that don't contain specific values in which to create a raster with the ability to be analyzed in a weighted overlay.

I get once I have rasters, I can reclassify them to match to do the overlay, but what happens when one (or two or three) of the those layers contain only polygons that represent boundaries and not a range in values (ie. one layer containing 5 equally important zip-code boundaries, one containing 3 equally important zone features, one containing 44 equally important census tracts, etc.)?

It was also requested that I make it look cool like a heat map.

Can this be done with a weighted Kerenel Density (then how do you weight kernel density...create a "population" field with ranked values 1 thru 8?)?

I am using ArcGIS Desktop 10.5 Advanced with most all extensions.


texts:
I would like to write only specific fields to an excel export. My current script is written as this:

fc = r"C:/Users/maps/table.xls"
xlsLocation = r"C:/Users/maps"
xlsOutputName = "SanFrancisco"
xlsPath = xlsLocation + '\\' + xlsOutputName + '.xls'
keepFields = ["name", "UID", "notes", "records", "surveyname"]

# I wrote the following, but it returned error 999999: Error executing function. I assume this maybe because the listfield function can't read xls?

fields = arcpy.ListFields(xlsPath)
dropFields = [x.name for x in fields if x.name not in keepFields]
print "Deleting extraneous fields in xls"
arcpy.DeleteField_management (xlsPath, dropFields)

Not certain about how to proceed. Our projects requires xls, so no other format can be considered. The excel file contains over 50 fields. I am using python 2.7


texts: Is it possible to distinguish between points with the same coordinates?

texts:
I have several points with the same coordinates in one layer and of course they are located on top of each other.

Is there a way to separate these points automatically (not manually)? I don't want to label them differently nor can I change the coordinates (although there's actually no other way...). But is there a trick or a tool I can use? I'm working with QGIS 2.0.1.


texts:
This is a follow on from previous question. You will see I am a new PostGIS user.

I have looked at the validity of the topology in my data and find around 1700 topology issues in the landcover (polygon) dataset. These are all "Ring Self_intersection".

Using some useful information in this presentation by Paul Ramsay(page 20), I have attempted to buffer my data accordingly:

  UPDATE schema.data
  SET the_geom=buffer(the_geom,0.0)
  WHERE isvalid(the_geom)=false and isvalid(buffer(the_geom, 0.0))=true;

But I get the response:

  NOTICE:  Ring Self-intersection at or near point (edited)
  ERROR:  new row for relation "sgm_buffer" violates check constraint
  "enforce_geotype_the_geom"

  ********** Error **********

  ERROR: new row for relation "sgm_buffer" violates check constraint
  "enforce_geotype_the_geom" SQL state: 23514

I have also tried just:

  UPDATE csgn.sgm_buffer SET the_geom=buffer(the_geom,0.0)

and get:

  ERROR:  new row for relation "sgm_buffer" violates check constraint
  "enforce_geotype_the_geom"

  ********** Error **********

  ERROR: new row for relation "sgm_buffer" violates check constraint
  "enforce_geotype_the_geom" SQL state: 23514

So do I assume that this buffer technique will not work on my issues? Or am I doing something wrong?

Mr Ramsay then goes on in the same document to suggest an additional technique to deal with "figure of 8" polygons. I don't entirely understand how to write this out in full in code (or what it is doing and whether it will help me):

  SELECT
          ST_BuildArea(
              ST_Union(
                  ‘LINESTRING EMPTY’,
                  ST_ExteriorRing(
                      ‘POLYGON((...))’

So.... can anyone help? This is proving a real nightmare with my tight timescales and limited skills in PostGIS.


texts: How to properly extract point data from multi-raster grib file in R?


texts: Getting all Null when spatial joining lines with polygons (ArcGIS Desktop 10.5.1)

texts:
I'm trying to append the attributes of the lines (0 or 1) to all the polygons that they intersect (take the mean if multiple lines exist), as shown below.
The polygons, lines and the DataFrame are all in the same Geographic coordinate system.

Then I ran a spatial join (see below) with the polygon as target feature, and the line graph as join feature.

The resulting attribute table, however, gives all Null values for all appended attributes ("Test_1" is the field where I put 1 and 0 for the lines).

I saw there's a bug related to the spatial join tool which should be fixed in version 10.5, and the coordinate systems I used for both feature classes are the same.

It worked when I use the "Join and Relate" - "Join" - "Join based on spatial location", however, the spatial join tool in the toolbox still couldn't solve the problem. Would love to know why or the difference between these two tools.


texts: R raster::terrain function curvature function missing (version 2.6-7)

texts:
I am trying to apply a predictive model based on terrain to a geographic area. I am trying to do this using R and the package raster. The model was developed using GRASS, but I have a limited understanding of GRASS and thought it would make sense to learn how to do it using R. But the package raster currently lacks the functionality I'm looking for.

The model contains different properties of the terrain. Aspect and slope is no problems to calculate using the terrain() function. However, also tangential and profile curvature is included in the model, but those options are not available. I found a guide at R-bloggers (2012) giving examples of how to calculate curvature using “plan.curvature” and “prof.curvature”, but it seems like this function has been removed from the newer versions of raster. A Stackexchange post points out some principal problems of calculating curvatures, and I imagine that this might be the reason behind the removal of those functions.

However, even if there might be limitations of the method, I still need to apply it to the model I have. Is there any way of calculating curvature using the current version of the  Raster package, or do I have to rely on an older version?


texts: Using the new OS Maps API (XYZ Version) and I keep getting 404

texts:
My request URL is this -
https://api.os.uk/maps/raster/v1/zxy/Outdoor_3857/5/47/-765.png?key=apiKey
I'm passing a valid API key, but all I'm getting is a 404. Any ideas on where the get further information as to why it's a 404?
OpenLayers 6 in an Angular 9 application.
if (this.osMap == null) {
  console.log('centreLng:' + data.centreLongitude + ' centreLat:' + data.centreLatitude);
  this.osMap = new Map({
    target: 'map',
    layers: [new Tile({
        source: new XYZ({
          url: serviceUrl + '/Outdoor_3857/{z}/{x}/{y}.png?key=' + apiKey,
          projection: 'EPSG:27700',
          tileGrid: tilegrid
        })
      })

    ],
    view: new View({
      projection: 'EPSG:3857',
      center: fromLonLat([data.centreLongitude, data.centreLatitude]),
      zoom: 12
    })
  });
}


texts:
PlanetScope basic scene images are 3.7m gsd (nadir) and orthorectified scenes are resampled to 3m per pixel.

How do you improve spatial resolution?

Does this processing affect the characteristics of the image?


texts: Calling a function by its name in Earth Engine code-editor

texts:
I have a module where I process some raster images according to one of different functions/ statistical models.
I am trying to find a way to allow the user to select which function to use. I have tried storing the functions as variables or properties but this doesn't work.
So my question here is how can I call a function using its name (string format)?
Any other workaround to dynamically choose which function to launch?

Here is the link to a simplified example.

Here is the code for the example:

//first possible function to use
var sumFun = function(a)
{
  return ee.Number(a).add(3)
}

// second possible function to use
var subtFun = function(a)
{
  return ee.Number(a).subtract(2)
}

//Storing functions as properties (DOESN'T WORK)
var modDic = ee.Dictionary({'sumFun':sumFun, 'subtFun': subtFun})
print(modDic)
var fun=modDic.get('sumFun')//calling function by its name
print(fun(4)) //doesn't work either atm

EDIT: corrected cast as dictionary


texts: Optimize Earth Engine script script profiler

texts:
I created a workflow where I try to apply an algorithm over a selection of transects (line geometry) and images in a ImageCollection. I noticed that processing time goes up very quickly when the amount of transects is increasing and I'd like to explore the possibility to imrpove the script (and processing time). I started looking into the profiler and I am thinking where to start digging and how to approach this issue.
The workflow roughly follows this procedure:

Create an imagecollection (SpectralUnmixing + Billateral filtering on images): for now 1 year of landsat observations (+/- 20 images)
iterate over all transects (+/- 250 with a length of 30km): break up  the transects in point, reduceRegion for each point and extract pixel values at transect.
Build Arrays with for each pixel the extract values of Spectal unmixing output
Do some array calculations such as first/2nd order derivatives for peaks & valleys on the transect
return the required output.

https://code.earthengine.google.com/8edf8f719c2f8ae01b1d8629bcdf9dde
So I was expecting that most of the computation time would be in the reduceRegion (step2) and billateral filtering (step 1) and not so much in the Array computations (step 4). So when I look at the profiler I see that most computation time is indeed plumbing, loading assets and applying image algorithms to it (related to billateral filtering function I assume). Yet, most of the memory goes to the algorithm array group, which I assume are the steps that fall under number 4.
Does the peak memory below in the image indeed indicate that my array computations can be optimized or is it mainly the image computations that could be improved? And where in the alogrithms is there a possibility to improve computation time spent on overhead/plumbing?
The screen shots are from tests applied to 1 image and 1 transect.


texts:
I have a GeoTIFF and I am trying to overlay a scale (miles or km) using GDAL.

Is there a function in GDAL that will do this?

I believe ArcPy has this functionality, but I am guessing GDAL does not - but wanted to check before I move on to my next question.

If there is no function for doing this - What would the best way be to get the scale?

My thoughts are, since I am using gdal.translate with projWin = [-2053514.581, 2702991.983, -2035515.049, 2691426.327], I can get the scale because these values are meter values.  Then I just need to create an overlay image based on that.


texts: resetStyle in geoJSON pointToLayer circleMarkers

texts:
I am creating a geoJSON layer, with pointToLayer returning circleMarkers. In that layer's oneachfeature, I am triggering a style change on hover, and this works fine. However, I cannot reset the style on mouseout. It all works with lines and polygons, but not with points created as circleMarkers via pointToLayers.

How can I reset the style of such circleMarkers on mouseout?


texts:
I have installed OSGeo4W following this: Installing QGIS with FileGDB?

Now I have QGIS Desktop 3.8.3 and gdal-filegdb 2.4.1-1 installed on windows 10.

I open directory that include the following file extensions: gdbindexes, gdbtable, gdbtablx, spx, atx, freelist , selecting ESRI FileGDB as described here too Installing File Geodatabase (*.gdb) support in QGIS?

Error "Invalid data source" is displayed.

How to open these files in QGIS?


texts: Function to remove borders and clouds from SENTINEL2

texts:
I've been trying to test SENTINEL, and found some functions to clean the clouds and edges of the images, but I do not understand exactly how this function works, and I'm very curious to understand better.

Which are, for example:

var cloudMask1 = function(image){
  var cfmask_conf = ee.Image(image.select(['B1']));
  return image.mask(cfmask_conf.lt(1600));
};

//Função que aplica uma máscara que remove os pixels com nuvem
//Essa função recebe cada imagem landsat e retorna a imagem sem nuvem
var cloudMask2 = function(image){
  var cfmask_conf = ee.Image(image.select(['B1']));
  return image.mask(cfmask_conf.lt(2700));
};

Could you explain me? Why the 1600 and 2700 values?


texts: Using ogr2ogr to save KML file including styling of lines/other features

texts:
I am trying to create an automatic procedure for using ogr2ogr to read a PostGIS table and output a KML file containing line features.
However, I need to define the styling of these lines including both width and color. Potentially I could save the styling information as attributes if that could help.

Is this possible with ogr2ogr?


texts: ArcMap Raster Calculator with Python code

texts:
I'm trying use raster calculator that can export a raster to do different average values for for multiple layers. For example, I want to do averages for only the values and layers that are greater than 0. I'm trying to input Python code but it keeps failing. Here is a video where I try and explain my project if you need more context. https://www.youtube.com/watch?v=bt6ecPGjjDA&t=24s

I've also tried putting single quotes on the line but there is still a syntax error.

>>> import arcpy
... from arcpy.sa import *
... f1 = arcpy.Raster("S:/npenme2/Faunalyze\Feature_KB3_Right_Radius_14497/value")
... f2 = arcpy.Raster("S:/npenme2/Faunalyze\Feature_KB1_Right_Radius_14471_RC11/value")
... f3 = arcpy.Raster("S:/npenme2/Faunalyze\Feature_KB4_Right_Radius_14510/value")
... if value > 1:
...     outraster = (f1 + f2 + f3)/3
... else:
...     outraster = (f1 + f2 +f3)/3 - 1
...     f123 = f1 + f2 + f3
... outraster = Con(f123 > 0, f123 / 3, f123 / 3 - 1)
... outraster.save(S:/npenme2/Faunalyze/rasterCalc)
...
Parsing error SyntaxError: invalid syntax (line 12)
>>>

I've also tried SQL on the actual raster calculator tool with this query. But it too does not work:

Con("Feature_KB3_Right_Radius_14497" + "Feature_KB4_Right_Radius_14510" + "Feature_KB1_Right_Radius_14471_RC11">0,(("Feature_KB3_Right_Radius_14497" + "Feature_KB4_Right_Radius_14510" + "Feature_KB1_Right_Radius_14471_RC11")/3), ("Feature_KB3_Right_Radius_14497" + "Feature_KB4_Right_Radius_14510" + "Feature_KB1_Right_Radius_14471_RC11")/3-1))


texts: Concatenate attributes ignoring NULL/missing values with FME

texts:
I've got a table with different attributes (named "source:geommetry", "source:name", "source:power"...) I'd like to concatenate to one attribute named "source_all" using FME Desktop. A blank space is supposed to be the delimiter.

The issue is that sometimes attributes are NULL which leads to an output with two (or more) blank spaces somewhere in the resulting string.

I used the StringConcatenator-Transformer (but wondering if it's possible to implement an if-statement.

Another try was using the ListBuilder and the ListConcatenator resulting in a completely emty string even if only one attribute is NULL. (Although "Drop Empty and NULL Attributes" is set to YES)

Any suggestions how to build a pretty string ignoring NULL attributes?


texts: Using ModelBuilder (and Iterator?) to rename multiple shapefiles files

texts:
I have thousands of shapefiles that I would like to rename, all are currently named "cYYYDDD_Union_Clip_Clip" and I would like to ultimately rename each to "YYYYDDDD_abc")? I've tried using Iterator in ModelBuilder with in-line variable substitution but as I understand both I can only substitute the existing shapefile name (e.g. the best I can do with Iterator and in-line substitution is maintain the input name as the output name).
Any one know what I can do to rename multiple shapefiles in an automated fashion?


texts:
When I try to generate chart by region, I get following error:
Error generating chart: String: Unable to convert object to string.
How can I revise my code to remove this error?
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_026047_20140216')
    .select(['B[1-7]']);
Map.addLayer(image, {bands: ['B5', 'B4', 'B3'], min: 0, max: 0.5});

// Define and display a FeatureCollection of 7 known locations.
var points = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(83.978710, 283202722)),
  ee.Feature(ee.Geometry.Point(84.008199, 28.218142)),
  ee.Feature(ee.Geometry.Point(84.2250, 28.371391)),
  ee.Feature(ee.Geometry.Point(83.947601, 28.215001)),
  ee.Feature(ee.Geometry.Point(83.895411, 28.588669)),
  ee.Feature(ee.Geometry.Point(84.042354, 28.523533)),
  ee.Feature(ee.Geometry.Point(84.085761,28.170022))
]);
Map.addLayer(points);

// Define a list of Landsat 8 wavelengths for X-axis labels.
var wavelengths = [0.44, 0.48, 0.56, 0.65, 0.86, 1.61, 2.2];

// Create the chart and set options.
var spectraChart = ui.Chart.image.regions(
    image, points, ee.Reducer.mean(), 30, 'label', wavelengths)
        .setChartType('ScatterChart')
        .setOptions(options);
var classNames = ee.List([ 'built_up', 'grassland', 'forest', 'water','snow','rock/barren','cultivation']);

// Define customization options.
var options = {
  title: 'Landsat 8 TOA spectra at three points near  City',
  hAxis: {title: 'Wavelength (micrometers)'},
  vAxis: {title: 'Reflectance'},
  lineWidth: 1,
  pointSize: 4,
  series: {
    0: {color: '00FF00'},// built_up
    1: {color: '00FF00'}, // grassland
    2: {color: '228B22'},//forest
    3: {color: '0000FF'}, // water
    4: {color: 'ADD8E6'}, // rock/barren
    5: {color: '808080'},//snow
    6: {color: 'FFFF00'}// cultivation
}};

// Define a list of Landsat 8 wavelengths for X-axis labels.
var wavelengths = [0.44, 0.48, 0.56, 0.65, 0.86, 1.61, 2.2];

// Create the chart and set options.
var spectraChart = ui.Chart.image.regions(
    image, points, ee.Reducer.mean(), 30, classNames, wavelengths)
        .setChartType('ScatterChart')
        .setOptions(options);

// Display the chart.
print(spectraChart);

code link: https://code.earthengine.google.com/0227789fa9d4b62979aa111d2d1e99ab


texts: How to create an exact tracking route from GPS data with animated icon?

texts:
Probably this is not purely GIS related question, its related to developing issues. I am just asking if there is any kind of suggestions from Gurus. I am building an tracking application where GPS data is receiving after every 15 seconds. I am using Leaflet, PostGIS and OSM solution. The GPS longitude and latitude are not on the road every time. I have illustrated 3 cases of tracking route here.

Case One:

Consider the following Figure One. Tracking route starts from Blue Marker (Marker 1) and ends at Red Marker (Marker 2). The actual tracking route should be the Blue line following the Yellow Major Highway Roads of Map. The green dots are the GPS positions which is receiving after every 15 seconds simultanously. If I simply draw a polyline according to the GPS points the tracking route will be the Green Line. As it is not the real route of vehicle there are so many disadvantages during routing and distance analysis.

Case Two:

Now Consider the following Figure Two for case two. Here I have optimize the GPS positions and calculate a nearest locations of GPS points on road. The Red dots are the false positions of GPS data. If we see deeply the second, fourth and fifth points are not in right position. Considering these three points, they were calculated and plotted from nearest road but not on the right position. They should be on Blue Line. In most road junction this can be happened. If I draw a tracking polyline it should be same as Red Polyline of Figure two. This returns better result than case one but not solving the issue accurately.

Case Three:

Now consider the following Figure Three where I have implemented both false position of GPS data and Routing. Its return more better result than case one or two but not solving the issue exactly. If we see deeply the routing result not showing accurate tracking path due to those three wrong false points and the result returns same as Gray Polyline of Figure three. On second point its going to a local road and again returns to the major highway. For the fourth and fifth position its following another major highway avoiding the accurate vehicle path.

What will be the best solution to get the accurate tracking path. What will be the best optimization in this case. I want to build up a solution similar to the following example of instacart site where vehicle icon is following the route and shows with an animation that show vehicle icon moving forward on road.

https://www.mapbox.com/blog/tracking-orders-live-instacart/


texts: Converting annotations to feature labels in ArcMap 10

texts:
The Situation:
I am organizing the labels for street names. In order to get the positioning correct, I converted the feature labels to annotations. This has allowed me to hand-craft the positioning of the labels, which was the desired result.
However, (to my knowledge) these annotations are now unable to scale with the feature.
For example: My 1:1 Print zoom looks excellent, but zooming in on a feature, the annotation remains at Print scale and appears giant at the new zoomed scale.

The Question:
Can I convert these annotations back to feature labels and maintain their new hand-crafted positions?

Details:
I am using a "Basic" license. This means I do not have access to the Feature-Linked Annotation.
It is not a problem if the Street Name labels/annotations remain separate from the Streets Feature. The main issue is keeping the scale of the label-to-street ratio when I zoom in or out.

Summary:
I want to place each street label by hand and have those labels scale accordingly with the Streets feature layer. (Add and Remove labels with free range of motion)


texts:
I have a beautiful grid/polygon, each square having its own "ID", and several observations for each square. But, I have three categories of observation, each with its own row, matching the same "ID". Here "gid" and "id_carreau".
I would like to have the three different "NUMCOUNTS" on the same row.
My final goal is to create a pie chart for square of the grid, proportional to these categories.
Here is my question, how do I merge these rows? see selected.


texts: Querying GeoJSON polygons from Django models in Leaflet?

texts:
I'm developing an application in GeoDjango and using Leaflet to show my data.I have loaded my data successfully into the map and added a sidebar where I have my textboxes to input values and a search button. I want to filter data depending on the inserted parameters in the textboxes.
I found an example at http://ipasic.com/article/query-and-filter-leaflet-map-layers/ which uses point data.
Was wondering how to handle polygon data.

Textboxes section

LR No:

            <input type="text" id="lrno" value="" size="30" />

Parcel ID:

            <input type="text" id="idno" value="" size="30" />

Parcel No:

              <input type="textbox" id="parcel_no" value="" size="30">

          <button id="searchbutton" role="button" class="btn btn-primary">Search</button>
          <button id="resetbutton" role="button" class="btn btn-primary">Reset</button>

Results

Map.js

    var parcel = L.geoJson();
    var dataurl = '/data/';
    $.getJSON(dataurl, function (data) {
    parcel.addData(data).setStyle(featureStyle);
    parcel.eachLayer(function (layer) {

layer.on('click',function(e){
    var layer = e.target;
    this.setStyle({
        weight: 4,
        opacity: 1,
        color: 'orange',
        dashArray: '3'
    });
    if (!L.Browser.ie && !L.Browser.opera) {
        this.bringToFront();
    }
    info.update(layer.feature.properties);
    map.fitBounds(e.target.getBounds());
});
layer.on('mouseout',function(e){
    parcel.setStyle(featureStyle);
});

});});

I can't really figure out what function I should use for this..
Any help ?


texts: Writing custom Prepare Scripts for Open Data Cube?

texts:
How do I write custom prepare scripts for the Open Data Cube format ?

The docs say the following:

  We expect that many new Data Cube instances will require custom prepare scripts to be written. It is generally a straightforward task of mapping metadata from one form to another and writing out a YAML document.

That's it.

Also the official training center is not helpful as it does not say anything about it.

Can someone explain this a bit more in detail or does someone have better sources or tutorial(s) other than analyzing the  non-commented sample prepare scripts ?


texts: Topology - To be, or not to be, that is the question

texts:
Topology, I mean "real" topology, is not a new concept and it's available since a long time in GRASS gis. I won't speculate the reasons myself but, after a long period of spaghetti model dissemination, recently some attention is leaning towards Topology implementation. One tool of major importance is Postgis2.0 and one great example of a kind of topology in action is openstreetmap. Qgis developers is also generating some discussion about the subject.

My experience with Grass was great, there is nothing more pleasant than two areas sharing the same boundary and "edit time" topology check. This presentation also describes great advantages of using real topology.

Now the question:
Considering web mapping and desktop gis and with a nice argumentation based on personal experience, industry requirements, tools available, etc. Do you think that real topology is to replace spaghetti model and be widely adopted in a near future?


texts: Adding QGIS coordinate, scale, magnifier and rotation toolbar into custom PyQGIS canvas

texts:
I was wondering if it is possible to add a toolbar with coordinate, scale magnifier, and rotation such as the one in QGIS Desktop to a custom PyQGIS canvas on a standalone application.

below is a simple code for map canvas in Main Window
from PyQt5 import QtCore, QtGui, QtWidgets

class Ui_MainWindow(object):
    def setupUi(self, MainWindow):
        MainWindow.setObjectName("MainWindow")
        MainWindow.resize(800, 600)
        self.centralwidget = QtWidgets.QWidget(MainWindow)
        self.centralwidget.setObjectName("centralwidget")
        self.gridLayout = QtWidgets.QGridLayout(self.centralwidget)
        self.gridLayout.setObjectName("gridLayout")
        self.widget = QgsMapCanvas(self.centralwidget)
        self.widget.setObjectName("widget")
        self.gridLayout.addWidget(self.widget, 0, 0, 1, 1)
        MainWindow.setCentralWidget(self.centralwidget)
        self.menubar = QtWidgets.QMenuBar(MainWindow)
        self.menubar.setGeometry(QtCore.QRect(0, 0, 800, 26))
        self.menubar.setObjectName("menubar")
        MainWindow.setMenuBar(self.menubar)
        self.statusbar = QtWidgets.QStatusBar(MainWindow)
        self.statusbar.setObjectName("statusbar")
        MainWindow.setStatusBar(self.statusbar)

        self.retranslateUi(MainWindow)
        QtCore.QMetaObject.connectSlotsByName(MainWindow)

    def retranslateUi(self, MainWindow):
        _translate = QtCore.QCoreApplication.translate
        MainWindow.setWindowTitle(_translate("MainWindow", "MainWindow"))
from qgis.gui import QgsMapCanvas

if __name__ == "__main__":
    import sys
    app = QtWidgets.QApplication(sys.argv)
    MainWindow = QtWidgets.QMainWindow()
    ui = Ui_MainWindow()
    ui.setupUi(MainWindow)
    MainWindow.show()
    sys.exit(app.exec_())

I have no idea on how to execute this. I have searched google for PyQGIS Canvas Toolbar but it does not yield the results I want.

The end result would be something similar to this (ignore the layer on the canvas)
Wondering if anyone could point me in the right direction. Thanks!


texts:
I do have a problem with QGIS.
I Have to create a map which shows different areas in a city. It is important, that there are labels with names of the streets, so that my clients can work with the map. As Background Map I have to use MapQuest.
My problem now is, that the displayed Background in the Composer is different from the background which is exported (Composer: only motorways are labeled, exported: every street is labeled but so small, that nobody can read it.
I’m able to label the streets from an external source (OSM) but then there are both labels displayed, the small ones and mine, and this looks terrible.
So I think I have two options, first I can export the map without labels, so I can use my extra ones, second I have to reach, that the labels are exported bigger. (both I do not know if and how it is possible)
Does somebody have an idea how to solve my problem?


texts: QGIS 2.10 for OS X dos not start on OS X 10.9.5

texts:
I have installed QGIS 2.10.1-1 from Kyngchaos.com following the provided instructions and installation requirements.
The install process of GDAL Complete 1.11, Matplotlib and QGIS (in this order) was ok, apparently.
Then, when I launch the QGIS software, after the splash screen the program crashes without any message during the loading of the GUI.

What is missing?


texts: Raster Calculator - filtering out values below 0.01

texts:
I am trying to use the raster calculator in QGIS to filter out depths lower than 0.01cm (in geotif which represent a flood map, respectively the depths), while keeping the original values of the raster. When I use the syntax I always get output with 0s and 1s, or I get the raster ranging from 0 up to max, where 0s represents the values below 0.01,and I am trying to avoid that. I am trying to get raster with values ranging from 0.01 up to max.

("HQ100_2D_PD@1" >= 0.01) * "HQ100_2D_PD@1"

I know that vectorizing the raster to points, then filtering out the values in attribute table will definitely work but I am looking for more elegant and faster approach using Raster Calculator.

Any ideas?


texts: Workaround for Intersection with points is not yet supported' GRASS Error

texts:
I am dealing with a line layer which I am unable to import properly into GRASS.
After some work solving dangles and invalid geometry issues I managed to import it, nevertheless when I try to use v.clean.break the software is giving this message:

"￼Intersection with points is not yet supported"

I have no idea how to deal with it.

Looking for help.

Link to file: line layer.


texts:
I have to create a map which is near at the MapQuest style, but for some reason I can't use a Lyer plugin, so I'd like to build the "Background" with files from geofabrik.

Are there qml styles (or other types of files which can be importet to qgis) on the internet, that I can download?


texts:
Is it possible to nest a QComboBox inside of a Qmenu within QGIS?

Below is my current code and an associated screenshot with the unexpected results.

(orange = current incorrect placement, green = ideal placement)

# CSV attribute editor toolbox layer population(s)
self.projCombo = QComboBox(self.iface.mainWindow())
self.projCombo.addItem('example slot')
self.projCombo.addItem('another example slot')
self.projCombo.addItem('the last example slot')
self.ShapeLayerList = QWidgetAction(self.projCombo)

# build the sub-menus for CSV attribute editor
self.popupMenu2A = QMenu( self.iface.mainWindow() )
self.popupMenu2A.setTitle("Shape Data")
self.popupMenu2A.addAction(self.ShapeLayerList)
self.toolButton2A = QToolButton()
self.toolButton2A.setMenu( self.popupMenu2A )
self.toolButton2A.setDefaultAction( self.editAttribStress )
self.toolButton2A.setPopupMode( QToolButton.InstantPopup )

self.popupMenu2B = QMenu( self.iface.mainWindow() )
self.toolButton2B = QToolButton()
self.toolButton2B.setMenu( self.popupMenu2B )
self.toolButton2B.setDefaultAction( self.editAttribProp)
self.toolButton2B.setPopupMode( QToolButton.InstantPopup )
self.popupMenu2B.setTitle("Grid Data")
self.popupMenu2B.addAction( self.editAttribProp )

# build the main CSV attribute editor dropdown toolbutton
self.popupMenu2 = QMenu( self.iface.mainWindow() )
self.toolButton2 = QToolButton()
self.toolButton2.setMenu( self.popupMenu2 )
self.toolButton2.setDefaultAction( self.editAttributes )
self.toolButton2.setPopupMode( QToolButton.InstantPopup )
self.popupMenu2.addMenu(self.popupMenu2A)
self.popupMenu2.addMenu(self.popupMenu2B)


texts: Code for SelectionChange Event on Comboboxes in VB.Net

texts:
I am programming 3 combo boxes so that each one feeds off the information from the previous combobox.
For example, in the first box they would select a year, then the second combobox would provide specific events for that year and then the third combobox would provide specific maps that are available based on the the first two answers.
I know the process is written in the Sub Change() location but I am stuck on how to relate them to each other. Can anyone possibly provide me some sample code for this procedure. Or maybe a website where I can go to learn how to write this kind of code. Any help would be great!!


texts: Can valuemap widget be set to allow new entries

texts: I am testing qfield for use in tree surveying. I wish to have a dropdown list of tree species. But I also need to be able to capture values not in the list. Can the valuemap widget be configured to allow entries not in the predefined list?


texts: How to make a column as a layer of label?

texts:
I have a table with many columns, some of them are "lon/lat" and "Name". After I draw objects on a map.

How can I make a layer by the content of the column "Name"?
How can I make a legend automatically? for example:

red ---2G

Blue---3G

Green--2G&3G


texts: Geoserver Geowebcache clear layer with parameters tile cache with Rest API

texts:
I'm using this command to clear tile cache of a certain layer:

curl -v -u admin:geoserver -H "Content-type: text/xml" -d "<truncateLayer><layerName>my-layer</layerName></truncateLayer>"  "https://xxx.xxx.com/geoserver/gwc/rest/masstruncate"

This is working fine. But layer 'my-layer' is a parametrized layer (Viewparams). Is there a possibility to truncate only tiles with certain parameters?


texts: Within the UK, do any commercial web services offer building level geocoding to a coordinate?  If so what level of accuracy is available?


texts: Using Time Slider and Time Format in ArcMap?

texts:
I have an .xls file that contains Lat/Longs as well as a time stamp in this format:

hh:mm:ss AM/PM

When I create a feature class shapefile out of this .xls file, ArcMap reads it in and displays it fine, but the TIME field just shows 12:00:00 AM.  What am I doing wrong here?  I don't understand why ESRI makes it so difficult to do really simple things, like read in times correctly!

Can anyone here show me the enlightened path to setting up my .xls properly or some other tricks that will enable me to enable the Time Slider feature?


texts: How to set/save the extent of a map in QGIS?

texts: I'm creating a bunch of maps with TimeManager in QGIS in order to create an animation. In order to overcome some of the limitations of TimeManager for dealing with raster files, I'm creating animations 30 frames at a time, which I can export and later combine the frames to make a longer animation. The problem with this approach is that it becomes important that the extent and scale of the map is identical for each of the 30 frame projects. So far, the only way I know to set the extent of a map in QGIS is to select the hand tool and drag it around, setting the scale numerically at the bottom of the map window. It would be far preferable if I could set the extent numerically as well, so that I could ensure that it is the same each time. Is there a way to do this?


texts: How to create label with variable size in ArcMap?

texts:
I'd like to create a circle label in ArcMap in the middle of each Polygon feature in my map document to shows the value of a particular attribute associated with that feature (like number of available stations). I'd like to represent the number of stations by the size of the label, so that the larger number of station in a particular feature is represented by the larger circle in that feature. Is there any possibility in ArcMap to do so?
any help would be much appreciated.


texts: Geographically Weighted Regression Steps

texts:
I am looking for either an answer or a resource that will outline the steps
I need to take to perform a GWR, including screening variables and the interpretation.

I have access to ArcMap with all extensions, R, GWR4 and good old excel and pencil and paper. I have my data in a .shp file and in an excel file with the XY coordinates for the polygon centroids calculated.

I have been trying to piece this process together by reading research papers, text books etc, but aside from seeing a multitude of the same matrix equations and one video tutorial using ESRI which I am not that confident in, I am really lost on where and how to begin the actual analysis.

Note:  If possible, I would like to avoid using the GWR algorithm in ESRI environment. I do have ArcMap for mapping if needed.


texts: Omitting certain style string labels in QGIS 3 rule-based labeling

texts:
I’m labeling airport taxiways. The taxiway labels are all strings but they tend to be one of three varieties and I would only like to show one of those. I’ve been trying to mess around with Rule-Based Labeling for this layer and haven’t had any success yet.

The strings can be one of the following format:

A, B, C...AA,BB,CC… A1, B1 etc
A.B, A.B.C, A1.B2.A, etc <— the use of the “.” In these makes them different from the first variety
$UNK<— this is when the feature has no known value but the field must still be filled in

Below I've attached an example of what the labels look like in QGIS.

I would like to only show the labels of the first variety.
I think the best way would be to create a rule that doesn’t print anything that equals $UNK or has “.” in it. I’m just unsure of how to write that as a rule or even if rule-based labeling is the right way to achieve my goal.


texts:
I have updated our ArcMap extension to work with ArcMap 10 - mainly just change reference from ESRI.ArcGIS.ADF to ESRI.ArcGIS.ADF.Local and set value of 'Specific Version' property for all references to false, and also use the new ESRIRegAsm.exe to register it.
So I am wondering now if I can build a single extension that will work with both 9.3.1 and 10.0 or do we need separate binaries for each product?


texts: Leaflet Layer Control appended to panel but duplicated

texts:
So I've used jquery append to attempt to move my leaflet group layer controls on to a panel.  This was successful, however the original still exists on the right side of the map object.  I've tried to remove this using jquery and css with no luck.  I think my problem is I have set an id to the control but the id is only applied to the second instance of it, so the original control remains...Any ideas?

$("<input type='button' name='layercontrol' value='EMEA' id='buttonstyle' /><br />").click(function() {
$("#emealayer").toggle("slow");
}).appendTo(".panel");

var EMEA = {
"EMEA":object
};
var layerControlEMEA = L.control.groupedLayers({},EMEA, {collapsed:false} );

layerControlEMEA.addTo(edit.map);
$('.panel').append(layerControlEMEA.onAdd(edit.map));

layerControlEMEA._container.setAttribute('id','emealayer');


texts: Searching for fields in geodatabase and adding them if they do not exist with ArcPy?

texts:
I have a list of files that need to be present in my featureclasses - myFields[].
I need to add these fields if they are not present in my FeatureClasses.

import arcpy

myField="Field01"
myFields = [field01, field02, field03]

env.workspace = "D:/test/data.gdb"
fcs = arcpy.ListFeatureClasses()

for f in fcs:
    fieldList = arcpy.ListFields(f)
    for field in fieldList:
         if field.name == myField:
             arcpy.AddField_management(fc, myField, "TEXT")

This works when I only have one field, but I have a list of fields! How do I adapt this when I have many fields and many different formats? Do I have to write a separate function for each field, or is there a cool python way of doing this?


texts: Java GeoTools: Creating Shapefile - FeatureCollections deprecated

texts:
I'm trying to create & populate Shapefile in Java using GeoTools.

Basically following this link (quite similar to the official documentation & other tutorials specifically for the point I'm asking below).

For below code I'm getting first warning that FeatureCollections is deprecated and error that add() cannot be resolved.

FeatureCollection<SimpleFeatureType, SimpleFeature> collection = FeatureCollections.newCollection();
...
collection.add(feature);

How can I fix that or any other way to create Shapefile using Java?


texts: RuntimeError: ERROR 010093: Output raster format UNKNOWN is unsupported

texts:
I have written a stand-alone Python script to perform an ExtractByMask operation on numerous raster files in ArcMap, that have a .bil extension format.  However, when I try to save my output, I get the following runtime error:

Below is what my code looks like currently:

    #Import necessary system modules
    import arcpy,os
    from arcpy import env
    from arcpy.sa import *

    #Location where raster (.bil) files are located.
    arcpy.env.workspace=r"D:\Mike_Daily_PRISM_PPT_from_2007_through2016\
    2007_Test"
    arcpy.CheckOutExtension('Spatial')

    #Location where output files should go after "for" loop has completed.
    outws=r"D:\Mike_Daily_PRISM_PPT_from_2007_through2016\2007_Test\Output"

    #Location of the shapefile used to mask the PRISM raster(.bil) files
    mask=r"D:\Mike_Daily_PRISM_PPT_from_2007_through2016\2007_Test\
    MissouriBound.shp"

    #Only search for files that have .bil extension.
    rasterlist=arcpy.ListDatasets("*.bil","Raster")

    #Print the list of raster files that have a .bil extension
    print rasterlist

    #for i in rasterlist, perform an ExtractByMask function on each raster.
    for i in rasterlist:
    outExtractByMask=ExtractByMask(i, mask)
    outname=os.path.join(outws, str(i))
    outExtractByMask.save(outname)

I believe the error has something to do with the format that it's trying to save the output format raster with.  [Note:  Actually, I would prefer to save the output in a TIFF format (.tif).]  Can anyone confirm whether that is the issue with my code, and how I can then modify it to save the output with a .tif extension?


texts: QGIS Gratuated Symbology gives differents results on every classification with Natural Breaks (Jenks)

texts:
When styling a polygon layer using the Graduated styling and the Natural Breaks (Jenks) mode, QGIS is classifying the data with different values.

So every time I click on "Classify" the values are different.

Why?

EDIT 1:

Using QGIS 3.4.4 on Windows 10


texts: Downloading satellite images in GeoTIFF format

texts:
I need to extract at a particular zoom level (17) and resolution (1280x1280) satellite images in GeoTiff (.tif) format that cover my polygons and points in my GeoJson file.
To obtain static images in .png format with Google Maps and Mapbox APIs there is no problem, but how can I get them in .tif format concretely?
The Mapbox and Google Maps APIs do not allow exporting in this format.
For example to export Point data as static images with Mapbox Static Images API, this simple script works very well (using geopandas).
gdf = geopandas.read_file(path_file)
for _, row in gdf.iterrows():
    image = open(name + ".png", 'wb') # image name is automatically generated
    image.write(requests.get('https://api.mapbox.com/styles/v1/mapbox/satellite-v9/static/'+ str(row['longitude']) +','+str(row['latitude'])+','+str(zoom)+',0/'+str(width)+'x'+str(height)+'@2x?access_token='+Mapbox_access_token).content)
    image.close()

Here is the image returned for latitude = 23.093865104007651 and longitude = 39.216685022463729 with the zoom and size quoted above. We can see the point located (like a wheel) in the center of the image.

I would just like to get the same result but have the image in GeoTiff format. Then I can process it, with rasterio for example. As it was advised to me in a previous question I tried with a WMS or WFS as with sentinelsat Python API, but I did not succeed in doing it exactly, especially with my need to have an image with a specific zoom.


texts: Intersect tool unable to find feature class in Python script

texts:
I keep receiving this error:

from this piece of my script:

# Process: Clip
arcpy.Clip_analysis(Full_Grid, NAR_U9, os.path.join(Output_Folder,
"Cruise_Grid_U9"), "")

Cruise_Grid_U9 = os.path.join(Output_Folder, "Cruise_Grid_U9")

# Process: Intersect
arcpy.Intersect_analysis(Cruise_Grid_U9, os.path.join(Output_Folder,
"Plots_Before_Dissolve_DELETE.shp"), "NO_FID", "", "POINT")

I can't understand why it's unable to find the feature used for the intersect. In the error message there is a '#' at the end of the file path. Is this just part of the error message or is it actually trying to look for that file path?


texts: how to split polyline where start/endpoints of the linestring touches or intersects with the nearby lines in qgis using python or postgis

texts:
I have a line layer and i want to find if the start/endpoints of the  linestring  touches or intersects with the nearby lines in qgis using python.

I have tried the following code but it is not working.

def linestouches(l_layer):
    ln_spIndex = QgsSpatialIndex()
    for ln in l_layer.getFeatures():
        ln_spIndex.insertFeature(ln)
        #print(ln['id'])
    for ln in l_layer.getFeatures():
        list1 = []
        print "****ln_id = ",ln['id']
        ln_geom = ln.geometry().asPolyline()
        ln_id = ln.id()
        print('The line id is',ln_id)
        # extract the segments seg_start, seg_end of the line
        X,Y = ln_geom[0],ln_geom[-1]
        print X,Y
        list1.append([X,Y])
        list1 = [val for sublist in sect_list for val in sublist]
        for i in list1:
            print( 'the i is',i)
            nearestgIds = ln_spIndex.nearestNeighbor(i,2)
            print('the nearest ids are',nearestgIds)
            for nid in nearestgIds:
                if nid != ln_id:
                    sel_feat = nid
                else:
                    pass
            print('the selected feat id is',sel_feat)
            l_layer.select(sel_feat)
            vertex_geom = QgsGeometry.fromPoint(i)
            sel_feat = l_layer.selectedFeatures()[0]
            sel_feat_geom = sel_feat.geometry()
            print('vertex_geom,sel_feat_geom',vertex_geom,sel_feat_geom)
            if sel_feat_geom.st_ intersects(vertex_geom):
                print('intersects')
            else:
                print('not intersects')

Please tell me what is going wrong here.


texts: ArcGIS Pro 2.4 .NET SDK: calling templatized class instance methods from Python?

texts:
I'm trying to do some learning with the ArcGIS Pro 2.4 .NET SDK in the Pro application's Python window:

import clr
import ArcGIS.Desktop.Core as ADC

for i in ADC.Project.Current.Items:
    print(i)

This works. Does anyone know the Python syntax for calling methods like Project.GetItems<T>(), if it's even possible?


texts: Trying to add explanatory variable to run logistic regression in R, to explain landcover change?

texts:
I'm trying to use landcover data from the NGS to run a logistic regression, using the same methodology from a study, Analyzing land cover change with logistic regression in R.

The problem is that I don't know how I can incorporate explanatory variables to account for the change. The data used in the study had variables such as proximity to roads and settlements, whereas the NGS data is simply descriptive of landcover type.

Do I need to manipulate the raster files in ArcMap in some ways before importing the csv into R, i.e. by somehow computing distance from each pixels to a road?

Or should I look for another set of data?


texts: Prevent deletes in a versioned table at the database level

texts:
I have a table in an Oracle 18c/10.7.1 geodatabase. The table is versioned with the option to move edits to base.
As far as I can tell, it's not possible to prevent ArcGIS Desktop users from deleting records, while still allowing inserts and updates, if the table is registered as versioned.
Question:
Is there a way to prevent deletes in a versioned table, while still allowing inserts and updates, by implementing something directly in the DB?


texts: Prediction mapping using temporal point data in ArcGIS Desktop?

texts: I have a temporal point data set of PM2.5 from 2002 to 2017.Based on this data sets I would like to generate some prediction (forecast) maps for 2020, 2025 and 2030. What would the be the main tools to derive these maps using GIS or ArcGIS?


texts: How to load Ordnance Survey OpenSpace in QGIS?

texts:
The version of gvSIG that I have on my Android tablet will load OS OpenSpace layers and I wonder whether or not it would be possible for QGIS desktop to do the same? And if so, how? I have my own OpenSpace API key.

This appears to be the (long) line in gvSIG's 'layers.txt' file that does the job of fetching the layers:

  101|Ordnance Survey (UK)
  OpenSpace;6[>],http://openspace.ordnancesurvey.co.uk/osmapapi/ts?FORMAT=image/png&KEY=#KEY#&URL=#URL#,image/png,11,0,200,0.0,0.0,0.0,0.0,1000000.0,1000000.0,EPSG:27700,2500.0:1000.0:500.0:200.0:100.0:50.0:25.0:10.0:5.0:4.0:2.0:1.0,
  ,1.1.1,88013ECD7A8F379FE0405F0ACA607F60:8824326E25F94E17E0405F0AC86047BC,http://www.prodevelop.es/maps/map.htm|http://www.gvsigmini.org/maps/map.htm


texts: Filter features within bounding box on provider side

texts:
I have a large polygon layer and want to filter only the features within a given bounding box (a other layer), BOFORE loading it into QGIS (on provider side). I tryed it with a SQL virtual layer but that consumes a lot of time:
select a.geometry, b.geometry from PolygonLayer a, BoundingBoxLayer b where ST_within(a.geometry,ST_envelope(b.geometry))

Is there a other faster and better way?
Currently I am using an Oracle DB as provider, but will be migrated to PostGIS in future.


texts: Replace field values with arcpy update cursor

texts:
I am a bit new for the arcpy. I have a database table so I can't use field calculator. so I would like to use Python script for solving the issue.

import arcpy
... # Your input feature class
... table = r'D:\werk_map.gdb\test_gw_bsl'
... # Start an update cursor and change values from 0 to 100 in a field called "your_field"
... with arcpy.da.UpdateCursor(table, "Thickness") as cursor:
...     for row in cursor:
...         if row[0] == "0,1":
...             row[0] = "0,5"
...         cursor.updateRow(row)

I just want to replace all “0,1” values in the thickness field with “0,5”.
The above written arcpy script is running without error but not changing the values in the table.
I think may be it has do something with the comma in the field value.


texts:
I'm trying to load some DRG data from the USGS into QGIS. The source data uses UTM coordinate and I can't seem to get it to overlay with WGS 84 vector data.  I (think) I've followed the recommended procedure to transform the coordinates but obviously I'm doing something wrong.

I'm loading the .jpg file and assume the associated .jgw file is loading with it.  QGIS reports translating the CRS to WGS 84 (project default) when I load it.  When I move the mouse over the image the coordinates look like utm values in the range from the world file and of course the data does not overlay with my vector data.

The data did not come with a .crs file.

I'm a new QGIS user.

World file contents:

1.464843750000
0.000000000000
0.000000000000
-1.464843750000
412500.732421875000
4190999.267578125000


texts: What GPS geographic coordinate system is this?

texts:
I received a couple of gps data on animal movements from a research paper but was unable to recognize the format and map it.

Any idea what format these GPS coordinates are?

The GPS locations should be in Singapore.

N 01’24 487, E 102 46.950
N 01’ 24 422, E 103 46.831
N 01’ 24 600, E 103 46.805
N 01’ 24 558, E 103 46.821


texts: Create layer based on input parameters and zoom in OpenLayers 5

texts:
I have a vector layer which i statically load all the features from the JSON file. What I want to do is to create layers dynamically based on the user input.

The user will choose district code, for example, and when button is clicked it will create the layer and zoom to the specified extent (call the ZoomToParcel function) . I haven't found any example in the docs.

var getDistrict, getTownVal, getQuarterVal; //this is user input (suppose they are dynamic)
var vector = new ol.layer.Vector({
    source: new ol.source.Vector({
    url: 'getjson.json', //loaded on client
    format: new ol.format.GeoJSON()
  }),style:style_simple
 //,extent:ol.proj.transformExtent([...], 'EPSG:4326', 'EPSG:3857') //This should be dynamic
});

var map = new ol.Map({
  target: 'map',
  layers: [roadlayer,vector],
    controls: ol.control.defaults().extend([new ol.control.ScaleLine()]),
    view: new ol.View({
    center: ol.proj.fromLonLat([33.33386, 35.14710]),
    zoom: 14,
    maxZoom: 19
  })
});

function ZoomToParcel(){ //the extent here should be dynamic
  let extent = vector.getExtent() || vector.getSource().getTileGrid().getExtent();
  map.getView().fit(extent,{duration: 1590},map.getSize());
}

HTML:

<input type="text"> District

<input type="text"> Town

<input type="text"> Quarter

   <button onclick="ZoomToParcel()" class="btn btn-primary">SEARCH</button>

I log the (e) inside the select event, see the dist_code inside the values_ array:

UPDATE:
Mike's answer helped with a bit modification, i make the if statement directly in the style. The filtering works fine, but by default the map should not show any layer (only baselayer).

var style_simple = function (feature, resolution) {

  let featureParcel = feature.get("parcel_nbr")+"";
  let featureDist = feature.get("dist_code")+"";
  if(featureDist == 1) {  //show districts with code 1
  return new ol.style.Style({
  fill: new ol.style.Fill({
    color: 'khaki'
  }),
   stroke: new ol.style.Stroke({
      color: 'blue',
      width: 2
    }),
    text: new ol.style.Text({
      text: featureParcel,
      font: '13px Calibri,sans-serif',
    fill: new ol.style.Fill({
    color: 'blue'
    })
    })
  });
}
}

As Mike suggested, Instead of re-loading the geojson for each change it is easier to filter the features on that property, in the style function so only the one we are interested in are displayed, and when setting the extent.

The user input is dynamic getDistrict (think of it as a drop down with 5-6 districts)

The point is to be able to update the map when input changes. I am a little confused where i should put the change event and how i can update the style accordingly


texts: Add columns with information OSM by osm2po?

texts:
I imported the europa.osm network within my spatial database using the osm2po tool. I use this with pgrouting.

How can I know the information provided by the tag address (city, country, etc) and not just the information that osm2po added by default?

I think I should change the osm2po.config file. But I have no idea how I should do this.


texts:
I created a shape in OpenJump and exported it to shape(ESRI) file  and then imported it to SQL2008 via shape2sql ,

the problem is that it shows bigger than what I want,

In OpenJump it shows 1meter for a small shape,but when I map it to real world map,it shows 100km.

I noticed that my problem is setting SRID or projection mode,
as I know there are some Projection Mode,
what is the best Projection Mode for drawing a building or factory.

Area about 20m^2 to max 100.000m^2 ,
Width about 4Meter to 1000Meter,
Length about 4Meter to 1000Meter.

It is for World Wide users , but every user just have a map of his Home/Factory.

What SRID should I have to set ?


texts:
I am developing an application using the dynamic layer.

However, I have trouble rotating the map (Not the same as rotating data frame by default).

See more details: https://geonet.esri.com/videos/1617

This is an example of the directory:

C:\Program Files (x86)\ArcGIS\DeveloperKit10.1\Samples\ArcObjectsNet\MyDynamicDisplayApp\CSharp

Link to download this example:

https://docs.google.com/file/d/0B0_zDbofmbfZOVFlSDhkQXVWM3c


texts: Programatically assigning NoData value in raster output?

texts:
I have written a tool that operates on raster data.  When the tool finishes it outputs a raster into a gdb, but the output does not have a NoData value listed under properties (input does, output does not).  I've tried myriad tricks to get a NoData value to export, no luck.

ESRI recommends exporting to a tiff to set a NoData value if a raster does not have one, but I would have to do that hundreds of times.  I'd like an arcpy command that will get the tool to output a raster with a NoData value.

I need to output a raster with a NoData value listed in the properties.

Can anybody help me?

Here is the end of my .pyt file.

output_raster = arcpy.NumPyArrayToRaster(tempras0, \

arcpy.Point(desc_ras.extent.XMin,desc_ras.extent.YMin), \
desc_ras.meanCellWidth,desc_ras.meanCellHeight)

outSetNull=SetNull(output_raster, output_raster, "VALUE < -10000")

outSetNull.save(parameters[2].valueAsText)

arcpy.DefineProjection_management(parameters[2].valueAsText,desc_ras.spatialReference)

arcpy.SetRasterProperties_management(parameters[2].valueAsText, nodata="1 -10.5")

return


texts: Adding x,y coordinates to Attribute table in QGIS

texts: If I use the 'Add coordinates to point' it adds the coordinates to in the file, however splits the data into the amount of objects in the original file. If I try to add separate columns in the attribute file calling them 'Coordinate X' or 'Y', then using the $x or $y expression, the output is just NULL. I need to coordinates but all in one file.


texts: Is there a difference between projecting shapefiles then rasterizing & rasterizing shapefiles then projecting?

texts:
I currently have a project where I rasterize shapefiles into simple presence/absence (1 & 0) rasters and then add each one to a running total to create a total richness raster at the end. The shapefiles are in WGS84 projection, but I need the final raster to be in Behrmann Equal Area projection.

What I want to know is if there is any difference in the final raster if I project all the shapefiles to Behrmann and then rasterize and merge each of them compared to rasterizing and merging, then projecting the final raster to Behrmann at the end?


texts: TypeError: unsupported operand type(s) for *: 'float' and 'geoprocessing Layer object'

texts:
I am getting a rather unclear type error, and after searching online I cannot find any solutions that help. I am trying to create my own python toolbox for wifi router signal loss, and I am trying to create a raster based on a few inputs that will populate from an equation. The error that I am getting lies on line 107, which is the equation line

fspl= (4 * math.pi * distmod * dist_float * hertz)/(2.99792458 * (10**8))
output_rast = fspl

I cannot figure out what is causing this error, initially it did not like me trying to use a raster instead of a float, so I changed the distance raster into a float, and now it does not seem to like the float being part of the equation. Here is the rest of my code incase the error stems from an earlier mistake. I cannot find anything that relates to my specific error, and I have not been able to find any way to to work around this issue.  Any suggestions?

import arcpy
import math

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox Label Property"
        self.alias = "Toolbox Alias Property"
    # List of tool classes associated with this toolbox
    self.tools = [FSPL, WAP_Buffer]

class FSPL(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Free Space Path Loss"
        self.description = "This python script tool will create Free Space Path Loss to determine the dB range output from the WAPs."
        self.canRunInBackground = False

def getParameterInfo(self):
    """Define parameter definitions"""

    param0 = arcpy.Parameter(
    displayName="Wireless Access Points",
    name="wireless_pts",
    datatype="GPFeatureLayer",
    parameterType="Required",
    direction="Input")

    param1 = arcpy.Parameter(
    displayName="Network Frequency Type",
    name="network_freq",
    datatype="String",
    parameterType="Required",
    direction="Input")
    param1.filter.type="ValueList"
    param1.filter.list = ["2.4 GHz", "5 GHz"]

    param2 = arcpy.Parameter(
    displayName="Distance Raster Float",
    name="dist_float",
    datatype="GPLayer",
    parameterType="Required",
    direction="Input")

    param3 = arcpy.Parameter(
    displayName="Distance Raster Units",
    name="units",
    datatype="String",
    parameterType="Required",
    direction="Input")
    param3.filter.type="ValueList"
    param3.filter.list = ["Feet", "Meters"]

    param4 = arcpy.Parameter(
    displayName="Output Raster",
    name="output_rast",
    datatype="GPRasterLayer",
    parameterType="Required",
    direction="Output")

    return [param0, param1, param2, param3, param4]

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""

    return

def updateMessages(self, parameters):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

def execute(self, parameters, messages):
    """The source code of the tool."""

    #Get inputs
    wireless_pts = parameters[0].value
    network_freq = parameters[1].value
    dist_float = parameters[2].value
    units = parameters[3].value
    output_rast = parameters[4].value

    #Create expression

    if network_freq == "2.4 GHz":
        hertz=2400000000
    else:
        hertz=5000000000

    if units == "Feet":
        distmod=0.3048
    else:
        distmod=1

    #equation

    fspl= (4 * math.pi * distmod * dist_float * hertz)/(2.99792458 * (10**8))
    output_rast = fspl

    return


texts:
is there a way to load MBTiles file using GeoTools?

I need to work with huge map and so far I have tried PNG + world file which was unusable (slow speed due to large files - 150MB). So I am thinking about using some sort of database and loading it by tiles would be the best way. I can easily create MBTiles (using MOBAC).

I have found MBTiles class in documentation http://docs.geotools.org/latest/javadocs/org/geotools/mbtiles/mosaic/package-summary.html but there is only a little explanation.


texts: Getting error message "An error has occurred while executing Python code: ZeroDivisionError: float division by zero"

texts:
I am getting the following error while executing Python code:
An error has occurred while executing Python code:

ZeroDivisionError: float division by zero
Traceback (most recent call last):
  File "C:/Users/goodwin/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py", line 120, in request
    self.setup_map()
  File "C:/Users/goodwin/.qgis2/python/plugins\openlayers_plugin\openlayers_layer.py", line 170, in setup_map
    qgisRes = rendererContext.extent().width() / self.targetSize.width()
ZeroDivisionError: float division by zero

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

Python Path:
C:/PROGRA~1/QGISES~1/apps/qgis/./python/plugins\processing
C:/PROGRA~1/QGISES~1/apps/qgis/./python
C:/Users/goodwin/.qgis2/python
C:/Users/goodwin/.qgis2/python/plugins
C:/PROGRA~1/QGISES~1/apps/qgis/./python/plugins
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
C:\PROGRA~1\QGISES~1\bin\python27.zip
C:\PROGRA~1\QGISES~1\apps\Python27\DLLs
C:\PROGRA~1\QGISES~1\apps\Python27\lib
C:\PROGRA~1\QGISES~1\apps\Python27\lib\plat-win
C:\PROGRA~1\QGISES~1\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGISES~1\bin
C:\PROGRA~1\QGISES~1\apps\Python27
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGISES~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/goodwin/.qgis2//python
C:\PROGRA~1\QGISES~1\apps\qgis\python\plugins\fTools\tools
C:/Users/goodwin

I am getting this error when moving from the map screen to the print composer.  I am using QGIS version 2.14.3-Essen on Windows 10. After I get this error, the software starts running very slowly with frequent 'not responding' status.


texts: Importing Python modules into Manifold GIS?

texts:
I am using Manifold 8, with the IronPython 2.6 distribution.

I've gotten Manifold to at least be able to run IronPython scripts (courtesy of http://www.georeference.org/Forum/t100925.11) by deleting the Scripts > Reference to IronMath.dll and changing the IronPython.dll reference to full file path.

However, despite being able to load up IronPython itself and import, say, urllib2 just fine, I can't do the same in Manifold. Instead I get:

  Exception has been thrown by the target of an invocation. No module
  named urllib2 [sic]

At the same time, I can do import math and then junk like math.sqrt(), and if I don't import math then math.sqrt() throws an error. This and the fact that Manifold knows the error suggests that the software should be able to import modules, but I don't know why it can't get urllib2 or cmd or whatever.

I've tried adding a reference to IronPython.Modules.dll in references, as well as a bunch of other DLLs that are in the IronPython 2.6 folder.

Given how little stuff there is about Manifold on here I'm not hopeful about finding a positive answer, but who knows. Maybe someone out there will have an idea or relevant experience that might assist despite not knowing Manifold GIS specifically.

Thoughts?


texts: Measure length of each segment for a polygon in Postgis

texts:
I have a polygon table in PostGIS and need to calculate the length of each line segment.

How can I achieve this? I have done it in QGIS but would prefer to do it in the database. I have found a solution in QGIS to do this.


texts: QGIS plugin window qt creator formatting is ineffective

texts:
I am developing a QGIS plugin. I have assigned the QGIS main window as the parent window of custom plugin using:

QtGui.QDialog.__init__(self, iface.mainWindow())

But after doing so the problem is that all the custom formatting in my plugin window which I did using qt creator is gone. I guess it is inheriting these properties from its parent window. Any ideas on how to restrain it from doing this?


texts:
A friend of mine is working on a small project to map receding North Pole summer sea ice minimums by decadal averages on a wall-paper map.

Unfortunately, it seems as though the map projection for his wall-paper map is non-standard (or closely related to the Robinson projection, but not exactly).

The map in question can be found here:
http://www.pinboardking.com/product/big-red-grey-worldmap/

Due to the fact that the sea ice min/max vectors in the North are very sensitive to projection warping/aberration, I'm trying to find the exact EPSG to be as accurate as possible.

How can the projection be identified or approximated?

Edit: Uploading a picture of said map that I had on hand for better reference.


texts: QgsDistanceArea.convertMeasurement not enough arguments error

texts:
I have setup a custom processing/python script to determine the shortest path on a network based on the code for Dijkstra's algorithm here,

http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/network_analysis.html#id4

Similar to the Road Graph plugin, I want to obtain the distance of the path as well. I've setup my code as below but I keep getting errors for the convertMeasurement function.

QgsDistanceArea.convertMeasurement(float, QGis.UnitType, QGis.UnitType, bool)
-> (float, QGis.UnitType): not enough arguments See log for more details

Not sure what I'm doing wrong. I've been researching and googling this problem but can't seem to find a solution. I've included some commented out lines of code as well to show what other options I've tried.

distance_measure = QgsDistanceArea()
distance_measure.setSourceCrs(crs)
distance_measure.setEllipsoidalMode(True)
#distance_measure.setEllipsoid('GDA94')
dist_unit = QGis.UnitType()
dist_unit_deg = QGis.UnitType(2)
#dist_unit
print dist_unit
print dist_unit_deg
#print type(dist_unit)
#QgsDistanceArea.convertLengthMeasurement(distance_measure, dist_unit)

dist_val = distance_measure.measureLength(rb.asGeometry())
print dist_val
QgsDistanceArea.convertMeasurement(distance_measure, dist_unit_deg, dist_unit, False)
print distance_measure.lengthUnits()
#print QgsDistanceArea.formatDistance(dist_val,QGis.Meters,5)


texts: Dissolving feature/field by interval using ArcGIS Desktop?

texts:
I have a huge dataset with flowpaths containing a catchment area for each line. The attribute table has ~ 780.000 rows, and thereby 780.000 features in the mapwindow.

Is it possible by using a field to merge the data in intervals, so you get a multi-part feature instead?

For instance merge all the features so those with a catchment area of 1-10.000 become one, from 10.001-20.000, 20.001-30000 and so on? The interval doesn't have to be 10.000 but can be another value.

I was messing around with the dissolve feature, but other than doing a manual population of a field and then dissolving by that attribute, I couldn't find another way to do it.

Basically I want what the symbology window does, but I want to do it to the dataset.

I have ArcGIS Desktop 10.5.1 with the Advanced License.

Edit:

import pandas as pd
import arcpy

fc = r'C:\Users\JGJ\Desktop\DatabaseVejle.gdb\StromningsVeje' #Change to match your data
bins = [-1,10000,50000,100000,200000,500000,1000000,5000000,10000000,100000000,1000000000] #Add/remove bins

df = pd.DataFrame.from_records(data=arcpy.da.SearchCursor(fc,['OID@','flow']),columns=['OID','flow'])
df['Interval'] = pd.cut(df['flow'], bins=bins).astype(str)

intervaldictionary = dict(zip(df.OID,df.Interval))

with arcpy.da.UpdateCursor(fc,['OID@','Intervalfield']) as cursor:
    for row in cursor:
        if row[0] in intervaldictionary:
            row[1] = intervaldictionary[row[0]]
        else:
            row[1] = 'Unknown interval'
        cursor.updateRow(row)

Ran into a snack, when processing a large dataset:

Runtime error
Traceback (most recent call last):
  File "<string>", line 7, in <module>
  File "C:\Python27\ArcGIS10.5\lib\site-packages\pandas\core\frame.py", line 950, in from_records
    values += data
MemoryError

Any suggestions on how to get passed this?


texts: Trouble Copying Data into File Geodatabase from Personal Geodatabase

texts:
I've been unable to successfully copy and paste personal geodatabase feature classes to a new file geodatabase in ArcCatalog 10.1. As the screenshot shows, the Data Transfer dialog box contains all the feature classes to be copied:

After clicking OK, the following error is generated:.

I've experimented with copying/pasting only one file at a time, and the same error continues to appear. Has anybody encountered this error when trying to paste feature classes from a personal geodatabase into a file geodatabase? If so, what steps did you take to troubleshoot it?

Thanks!


texts: Setup ArcGIS 10.6.1 SDK runtime in Visual studio 2017

texts:
I have ArcGIS 10.6.1 and I wish to setup ArcGIS Add-in in Visual Studio 2017 community edition.

I installed Visual Studio followed by ArcGIS Desktop. Then I installed ArcGIS 10.6.1 runtime SDK.

However, ArcGIS link is not visible in Visual Studio.


texts: Converting between coordinate systems but retain visual shape?

texts:
We have some data that is stored in GeoJSON format (using EPSG:4326), but now we want to display the data in the mapbox-gl-js, which only supports EPSG:3857 for display (support for other projections is not yet implemented). The data is arbitrary, non-earth geographic data.

How can I transform the coordinates from the old system, so that visually map in both systems look identical?

I do not care about actual geographical accuracy, but rather the visual accuracy.

As I understand, I should be working with Secant Projection but to me it just doesn't seem to work.


texts: How to find a point on map with GPS coordinates?

texts:
I need to import data from my smarphone using GPS.

The area is north of Italy so I believe I have to choose SR: WGS84 for the new project, haven't I?

Data are collected in the following csv table:

id - x - y
120 - 45,0770533 - 7,6610054
124 - 46,115925 - 12,098411

I use QGIS 2.6.1


texts: How can I view a .kml or .kmz file on a mobile device offline?

texts:
I need to view a KML/KMZ file on a mobile device in a remote area where there will not be service.  I have downloaded applications that have maps of the area; these maps will be available offline. However, I don't think they have a capability to import the KML/KMZ file.

What can I do to make these files accessible without pulling data from the network?

Keep in mind I also need GPS capability (all the mapping apps I've downloaded have it, so this shouldn't be a problem).

If that's not possible, can I save the text of the KML/KMZ and at least have that on my mobile device (IPhone 5)?


texts:
I have a sanitary sewer line feature class where each line has a unique identifier.

I have a table that lists the date a line was cleaned. It may have been cleaned multiple times so I can have multiple dates for each unique identifier.

I also have another table that lists the date a line was inspected. It may have been inspected multiple times so I also have multiple dates for each unique identifier.

I need to select each cleaned date and find an inspection date that lies within a 30 day range from when it was cleaned.

For Example,

From Table_Clean:

Line A is cleaned 1/1/2014
Line A is cleaned 4/1/2014
Line A is cleaned 8/1/2014

From Table_Inspected:

Line A is inspected 1/2/2014
Line A is inspected 4/2/2014
Line A is inspected 11/2/2014

I need to list that the first two were cleaned and inspected within 30 days. The third was cleaned but not inspected within 30 days.

Is this many to many relationship possible? And does anyone have any suggestions on how to symbolize this on a map? Obviously the same exact line feature may be represented twice being cleaned and inspected on different dates.


texts: Mapping datapoints over several floors of a building

texts:
I'm new to web mapping and am planning and investigating ways on how to best plot and visualize user flow through a building across several floors.

We are planning on tracking users location based on proximity to iBeacons. Each data point will track both time and location (based on proximity to the iBeacon). We're considering building a torque visualization to demonstrate how users flow through and use the building.

How would you recommend best visualizing this data and representing the multiple (4) floors to see user flow around the building (hopefully on one map).


texts: GDAL Rasterize with "where" option in Python

texts:
I've got a shapefile with set of points with a column "test_class" of String type in the attribute table.

I need to rasterize these points with setting pixel's value for each class and typically I use:
gdal_rasterize -burn 1 -where "test_class='classA'" points.shp raster.tif

and it works.
At now I need to automate the procedure in Python.
A simple test works fine:
import gdal
raster = gdal.Open("raster.tif", gdal.GA_Update)
shp = "points.shp"
OPTIONS = gdal.RasterizeOptions(burnValues=[1])
gdal.Rasterize(raster, shp, options=OPTIONS)

But when I try tu use "where" option:
OPTIONS = gdal.RasterizeOptions(burnValues=[1], where='"test_class"="classA"')

I'm getting an error:

ERROR 1: "classA" not recognised as an available field.

I was trying to use "where" with many combinations of characters ', ", without, and it still doesn't work.
I am open to any other solutions. I was tested rasterio library but it is very slow.
I use Python 3.7 and GDAL 3.0.2.


texts: How to create one shapefile from multiple rasters?

texts:
Ok, I have about ten maps that connect along their edges to each other linking rivers and properties.  When I link them all in Inkscape I end up with a large file that won't save (and a smaller one with 5 of the properties saves but is too big and slow in QGIS).

So I am thinking that I will load them, about three linked together at a time, into QGIS and create my shapefiles.  Is it possible in QGIS to then merge/link the shapefiles along their edges so that the rivers and properties line up and create one large shapefile?   additionally, is it possible to then take this large shapefile and georeference it to a base map or can that only be done with rasters?

thanks


texts: Counting the number of shared point within polygons

texts:
I have a set of points on a map, each with a unique location id, with one or more points within a number of polygons with unique polygon ids. I am trying to associate each location id with the number of shared points in its respective polygon.

To do so I used the intersect tool on the points feature and polygon feature, intersecting by polygon id. I then used the dissolve tool with the polygon id as the dissolve field and the location id COUNT in the statistics field. Then unchecked "create multipart features".

This leaves me with a table of points with the polygon id each point is in, and the number of points in each polygon. But removes the location id associated with each point. How can I keep the location id associated with each point?

Here is my table before being dissolved:

Here is my dissolved table, I want to simply keep the location id for each point, but getting the count removes that field.

I did something similar using spatial joins but need to see the differences in output using the intersect and dissolve tools.


texts:
I meet a problem while I try to calculate the length of a polyline in QGIS in the attribute table.

I created a new column with the function "$length" inside (type: real, 10 digits and 3 digits of precision). My projection is WGS84.

My problem is that the result is completely wrong!

For example, I used the measurement tool on the desktop and measured a street: I got 525m. On my attribute table I have 0.007...

I tried to change my projection but I still have the same problem...

Does anyone have an idea to help me ?