texts: QGIS: how to set transparency for OSM tiles layer?

texts:
I'm using QGIS 2.8.3-Wien in an Ubuntu14. My problem is the following:

I'm using the OpenStreetMap tiles to draw maps of some specific regions. The problem is that, for some urban/rural regions, these tiles offer very little information at non-urban regions.

By using the OCM Landscape layer, I can get more landscape information for rural areas, although the infos for urban areas is poorer.

Finally, I though I could merge both layers by setting some level of transparency on it.

However, QGIS does not offer a properties dialog for this kind of layer.

Is there any way of setting transparency for a tiles layer?


texts:
So I want to remove some markers that are created on a button click. If I use map.removeLayer() it removes the layer and also I set the array.length = 0, which resets the containing layer. However, if I run the function again, it shows the same amount of created markers and array elements (fine so far), BUT the _leaflet_id count goes up, as if the previous markers haven't been deleted.

For example,

I just figured this issue out because also the markers are overlapping at the same position, so they appear bold and the browsers slows down.

Any ideas how to solve this? Basically, if I clear the elements, I want to reset the variable and the visualization to zero, i.e the initial state of the map before I create any markers.

SOLUTION

I had to use markers.clearLayers() in order to remove all the layers inside the marker featureGroup. Before I only did map.removeLayer(markers), which does not delete the featureGroup content, but only removes the markers layer


texts: App for collecting points with hierarchical legend encoding

texts:
In the frame of a field campaign for land cover mapping, we try to develop an App which allows for the collection of points on a base map and fill a set of attributes. So far, our best solution is to use ArcGIS Collector because we often need to work offline. Using the domain, the attribute table is filled with pre-encoded values.

The problem is that we have a hierarchical legend with three levels. There are less than ten items at the first legend, then around 25 at the second level and more than 50 at the third. The drop down list is therefore very long at the third level, and if everything is possible we could have inconsistent labels between two levels. (e.g. grassland then cereals instead of grassland then graminaceous plant).

So, the question is : how can we filter the lists proposed to the field collector based on his previous choices ? In other words, how can we adapt the domain to the values of another field.

For instance, after selecting "forest" in the first drop down list, a new drop down list should be available with only "mixed forest", "broadleaved forest" and "needleleaved forest".

EDIT: there was already an ArcGIS Idea on this topic, so the question probably doesn't have a built in solution. Any workaround with ArcGIS or similar capabilities to ArcGIS collector with another app would be worth the bounty (as long as it works). To sum up, here are the main features :

1) provide a base map (high resolution image)

2) works offline (GPS but no Internet)

3) point capture tool with colour symbol

4) when a choice has been made based on a drop down list for the first field, a second drop down list becomes available conditionnally to the first choice, then a third drop down list becomes available conditionnally to the second choice

5) If possible not specific to a single mobile device

I can already do 1,2,3, and 5 with arcgiscollector. So my question is: how to get  conditional drop down list (point 4).


texts:
I'm nota GIS guy so I need some help.  I need to reproject from one SRID to another.  I'm trying to take global lat/longs (4326), reproject them to a local SRID (2336).  This is basically the line of code I'm using.  I keep getting back NaN.  I've tried reversing the lat/long.  Should I be using a different Z value?  What am I doing wrong?

DotSpatial.Projections.Reproject.ReprojectPoints(
    new[] { 29.4660863, -81.25784159999999 },
    new[] { 0d },
    ProjectionInfo.FromEpsgCode(4326),
    ProjectionInfo.FromEpsgCode(2336),
    0,
    1);


texts: How to extract zip code from reverse geocoding in PostGIS

texts:
I have a bunch of lat/lng coordinates that I can reverse geocode, after having installed PostGIS and the relevant Tiger files.

I throw in a query like this:

select reverse_geocode(ST_Point(lng, lat))
from plus_viewable
limit 5

And so far, get results like the following:

({0101000020AD10000041F2CEA10C2055C08AE76C01A1F54040},"{""(1343,,Franklin,Dr,,,Marietta,GA,30067,)""}","{""Cobb Pkwy S""}")

I just need to be able to extract the zip code.

How would I be able to do that?


texts: How do I make a wind vector map using wind data (u and v in one nc file) and overlap this wind vector map over sea surface temperature map in QGIS?


texts:
GIS student hoping for a bit of wisdom on the Path Distance tool.  I am doing a ski slope analysis, trying to find runs using cost grids (that incorporate slope, aspect, avalanche risk, etc.)  While I am familiar with the cost distance tool for finding the path of lowest cost, I wanted to see if there is a way to find paths that are: 1) of lowest "cost" in terms of my cost grid, 2) always sloping downhill (or at least remaining flat).  I have played with the path distance tool for a while now, and the documentation is limited.  Is this the right tool to use?  I've been setting the vertical factor parameters with my DEM to have a zero factor of 1, a low-cut angle of -90 degrees, and a high-cut angle of 0 degrees.  This produces a grid containing only one cell, with value 0.

Any thoughts or advice?


texts: Conditional rules for all columns of an attribute table in QGIS

texts:
is it possible in the attribute table of a vector, to use a condition for the whole table in a single method?

I have computed a correlation matrix between all the fields and it could be very useful to highlight all the values, for example, > 0.8.

I tried to use the Full row option and to add @value > 0.8 but it did not work..

Else I had to add a rule for each column and it could take a while..


texts: GeoServer, select feature, highlight it, and export it to PNG

texts:
Like any GIS Web application, you need to provide a facility to print a report with a part of map, and specially the position of a given ParcelID on the map with different zooms...

I have GeoServer, simply, I want a script that accepts a filter (something like ParcelID='ABC123' and returns me a PNG URL of this Parcel highlighted on the map with all the other parcels and layers...

I know the WFS with GetFeature and Filter return some row data about this parcel, and I know the WMS can return a map with CQL_FILTER but only for the needed parcel without any other parcels, so it's not enough for me.

Would you please help!


texts: Inserting search box using dataset column using CartoDB?

texts:
I have got some problems to add a search box that can runs searching points from a column of my dataset. I've read this article:

Linking search box to dataset column in Carto?

but I'm not able to resolve the problem beacuse I don't know HTML language, so I don't know what I have to insert, ore what I have to change!

Can anyone help me?


texts: Set arcpy command line parameter same as Tool parameter

texts:
I have a python script that works correctly when used in a toolbox but not from the command line. How do I duplicate the toolbox magic that happens when specifying a parameter Data Type?

There is GetParameterAsText() but apparently not GetParameterAs('FeatureLayer'), or GetParameterAs('Workspace'), ... in the GetParameter docs.

The parameters as handled in the script:

mxd = arcpy.GetParameterAsText(0)
clip_layer = arcpy.GetParameterAsText(1)
out_gdb = arcpy.GetParameterAsText(2)

In the toolbox I have the 2nd with Data Type set to Feature Layer.

When run from the command line this emits:

  arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
  ERROR 000840: The value is not a Feature Layer.

...while from the Tool it is:

  Executing: clipAllLayers Q:\Current_AOI.mxd Q:\current_aoi.gdb\Converted_Graphics Q:\clipped.gdb
  Start Time: Wed Feb 25 10:33:30 2015
  Running script clipAllLayers...
  layer_one
  layer_two
  ...snip...
  Completed script clipAllLayers...
  Succeeded at Wed Feb 25 10:33:55 2015 (Elapsed Time: 24.81 seconds)


texts:
I would like to create my own custom toolbar for leaflet.draw plugin but I don't know how to approach this solution.

To be more clear, I am trying to create a side-panel with buttons for each action, and with each onclick event I call the draw action (replace the default toolbar with a custom one).

Is that even possible?


texts: No result when performing Attribute Transfer Mapping for a line layer to new line layer in ArcGIS

texts:
I am performing Attribute Transfer Mapping in ArcGIS, but not getting results as desired.

My project is to survey the roads, bridges, culverts and habitation with some attributes by using a GPS device. I deployed my survey team and every week end they send the data in shapefile format, then my in-house team take the survey data reference and extracts the road center line on base images. Now, here the problem arises, when I perform the Transfer Attribute Mapping with matching fields, nothing is working and result is zero.

However, when I use the interactive transfer tool, it's working...

Can any one suggest me what to do in this situation.

See the below image...


texts: Getting mean values of reclassified raster classes inside a fishnet grid

texts: I have a raster DEM and a reclassified raster DEM  with different elevation ranges. I have a fishnet grid box (polygon). What I need to do is calculate the mean values of data lying inside those ranges for each grid. Say, I have set 3 elevation ranges 0-500, 500-1000 and 1000-8879m. Now for each grid, I need the mean of values of raster DEM lying inside those ranges. For eg. 2, 43, 24 lying in first elevation range should give their average value for a grid cell. Similarly, two other values are to be found accordingly for the other ranges. Hence I will have three values for a grid. It is possible that some grids may have only one or two values. The raster is 90m and grid is 25km resoultion.


texts: I cannot install QGIS 2.4 on my Windows 7, 64 bit machine

texts:
I am running Windows 7, 64 bit and cannot install QGIS 2.4 after several days of trying. I have also tried downloading 2.2 with no avail. I keep getting this error:

The text of the error is:

  NSIS Error

  Installer integrity check has failed. Common causes include incomplete download and damaged media. Contact the installer's author to obtain a new copy.

  More information at:
  http://nsis.sf.net/NSIS_Error


texts:
Ciao, sto creando una pagina in php con una mappa ... in OpenLayers Examples ho trovato due interessanti esempi:

esempio che mostra come usare diversi tipi di poupup (http://openlayers.org/dev/examples/popupMatrix.html)
esempio per inserire icone personalizzate (http://openlayers.org/dev/examples/markers.html)

Il mio problema è che non riesco ad integrare i due esempi...

dal primo:

echo"
        ll = new OpenLayers.LonLat($lng,$lat).transform(
            new OpenLayers.Projection('EPSG:4326'), new OpenLayers.Projection('EPSG:900913')
        );
        popupClass = OpenLayers.Popup.Anchored;
        popupContentHTML = '$T';
        addMarker(ll, popupClass, popupContentHTML, false, true);

dal secondo:

echo"
        ll = new OpenLayers.LonLat($lng,$lat).transform(
            new OpenLayers.Projection('EPSG:4326'), new OpenLayers.Projection('EPSG:900913')
        );
        var size = new OpenLayers.Size(50,50);
        var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
        var icon = new OpenLayers.Icon('{$row['enclosure']}',size,offset);
        markers.addMarker(
            new OpenLayers.Marker(
                (ll
                ),icon
            )
        );
";

Ho notato che nel primo addMarker richiama una funzione, che sembra simile a markers.addMarker

... sinceramente non so dove mettere ", icon" nel primo o nella funzione che richiama?

EDIT - Translation to English (Tradotto in inglese)

Hello, I am creating a php page with a map ... in OpenLayers Examples I found two interesting examples:

• Example showing how to use different types of poupup (http://openlayers.org/dev/examples/popupMatrix.html)

• example to insert custom icons (http://openlayers.org/dev/examples/markers.html)

My problem is that I can not integrate the two examples ...

the first:

echo"
            ll = new OpenLayers.LonLat($lng,$lat).transform(
                new OpenLayers.Projection('EPSG:4326'), new OpenLayers.Projection('EPSG:900913')
            );
            popupClass = OpenLayers.Popup.Anchored;
            popupContentHTML = '$T';
            addMarker(ll, popupClass, popupContentHTML, false, true);

The Second:

echo"
        ll = new OpenLayers.LonLat($lng,$lat).transform(
            new OpenLayers.Projection('EPSG:4326'), new OpenLayers.Projection('EPSG:900913')
        );
        var size = new OpenLayers.Size(50,50);
        var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
        var icon = new OpenLayers.Icon('{$row['enclosure']}',size,offset);
        markers.addMarker(
            new OpenLayers.Marker(
                (ll
                ),icon
            )
        );
";

I noticed that in the first addMarker calls a function, which looks similar to markers.addMarker

... I honestly do not know where to put "icon" in the first or in the function calling?


texts: Is there a way to save a text box's formatting and text so that you can insert it into a map with only a couple of clicks? For example, I have to insert a box that reads "CONFIDENTIAL NOT FOR PUBLIC RELEASE" on lots of maps, and it would save some time and annoyance if I didn't have to go through the process of changing borders, backgrounds, text size and alignment, and all that jazz every time.


texts:
Currently I am looking list of reputed (or top 50) GIS Journals and Publications across the world to present/submit GIS paper.

Found few links :

GIS Journals and Other Publications
GIS Online Publications

Any additions  will be great !!!


texts: Showing offline map downloaded from openmaptiles.com with OpenLayers?

texts:
I've downloaded .mbtiles file from https://openmaptiles.com/downloads and want to show the map in localhost. I've used following code but it did not show anything.

<!DOCTYPE html>

<title>Local Map</title>
<link rel="stylesheet" href="css/ol.css" type="text/css">

Go to map

Following scenarios were tested and non of them were successful

putting .mbtiles file and referring to it as url: 'pathtofolder/file.mbtiles'
putting .mbtiles file and referring to its parent folder as url: 'pathtoparentfolder'
extracting .mbtiles with mbutil and obtaining many many numbered folders with metadata.json file and referring to it as url: 'pathtoparentfolder'
extracting .mbtiles with mbutil and obtaining many many numbered folders with metadata.json file and referring to it as url: 'pathtoparentfolder/{z}/{x}/{y}.png'

What is the problem and how can I debug OpenLayers? In fact I hesitate between using ol.source.OSM and ol.source.XYZ.

It's noteworthy to mention I have many .png files extracted in that folders but they can not be displayed.


texts: There is not much information in the form of documentation on http://geonode.org/ development. I need to do some custom "App" development for geonode so that my users can have more functionality within the maps ie: search functions, geodesy transformations, etc.  My question is what is the recommend development path for this is it frontened (JS) or backend (django)..or something else?


texts: Does India Have a CORS (Continuously Operating Reference Stations) Network?

texts: Does India have a CORS (Continuously Operating Reference Stations) network? If yes, can anybody provide more details about it?


texts: Understanding the TimeWholeYear configuration in ArcGIS Area Solar Radiation

texts:
I am a bit confused with the interpretation of my results and it would be great to have some clarification!

I used the TimeWholeYear configuration (without choosing the {each-interval} option, so the result should be a single output for the whole year) for the 'Area Solar Radiation' tool in ArcGIS10. I set the analysis to run with intervals of 14 days, and 0.5 hours.

With that I understand that every 14 days, and every 0.5 hours of those days, the tool will measure the radiation that falls on each cell of the raster.

Now that I have the results, I am a little bit confused as to what they mean. Basically, I have 2 questions regarding how to interpret the results:

Is the result value in each cell an average or the sum of the radiation received by that cell in the whole year?
If I go to layer properties > statistics in the output raster, I can see the minimum, mean and maximum value of the global radiation. How do I interpret that?

a) it refers to the minimum, mean, and maximum radiation measured during the year (i.e. measured on one occasion every 14 days and 0.5 hours)?...

...or b) it is the min, mean, and max values amongst all the cells of the raster?

I'm sorry if these are a very basic questions!!

Thanks,
JJ


texts:
I would like to know how to disable the login window [PICTURE A] when I start ArcMap ?

In the code I am starting ArcMap, but then I get prompted with username and password dialog.

How can I disable this dialog ?

Keep in mind that I have the save username and password in the ArcCatalog connections properties so that is not that issue here.

var mxDocument = new MxDocumentClass();  // This statment is where the login screen pops up and arcmap starts
arcMapApp = mxDocument.Parent;
if (arcMapApp == null)
  return
arcMapApp.Visible = false;

PICTURE A


texts: Calculating field based on content from another table with no join in ModelBuilder?

texts:
I'm trying to work out how if (at all) it is possible to do a field calculation where there is no join or relationship between two attribute tables in ArcGIS 10.3.1. I've searched all over the place for a solution.

In my example, my model already selects some features based on an attribute, resulting in a featureclass with approximately 2000 features, with population figures assigned to them. I already have a table with the significant loadings that each of the population figures has towards a final score.

I basically want to be able to do the following calculation.

[Featureclass Field Value]*[Table Field Significant loading]

I'm trying to avoid "hard coding" the significant loadings into the field calculation so that I can just update the table, rather than having to change several parts of the model, when/where necessary.

I can't join or relate the tables within ArcMap as they have no common field.

I suspect that this might be able to be done with Python but it's been a while since I used it and I'm getting errors every time I try so I think I may be missing something fundamental, not seeing the wood for the trees.

My first attempt is this but IDLE gives me an errors after I define the parameters (ie: when I tell it I want it to multiply stuff)

import arcpy
import os

# Defined within model
SLO65Sing = arcpy.GetParameterAsText(0) #Significant Loading field 1
SLO65All = arcpy.GetParameterAsText(0) #Significant Loading field 2
ValO65Sing = arcpy.GetParameterAsText(0) #Value field 1
ValO65All = arcpy.GetParameterAsText(0) #Value field 2
FSO65Sing = arcpy.GetParameterAsText(0) #Factor Score field 1
FSO65All = arcpy.GetParameterAsText(0) #Factor Score field 2

def Calc_O65Sing(SLO65Sing*ValO65Sing)

Calc_O65Sing(!FSO65Sing!)

def Calc_O65All(SLO65All*ValO65All)

Calc_O65All(!FSO65All!)

Note: moved this over from Stackoverflow when I realised I had posted it in the wrong place!


texts: OpenLayers3: How to use zIndex for Layers

texts:
In my project I have the following situation:

layer A: zIndex = 1000;
layer B: zIndex = 2000;
layer C: zIndex = 3000;
layer D: zIndex = 4000;

I add the layers in this order:

layer A, layer C, layer D, layer B

I expect that correspondent to the zIndex the layers are in the following order:

layer A, layer B, layer C, layer D

But they are not! I get the order as they are added.

My questions:

Is the zIndex mechanism implemented - I'm using OL3.15.x and OL3.16.0?
How to set the zIndex to get the wanted result?
Are there any limitations on zIndex?
Should I add the layer before setting the zIndex?

Currently I'm porting a huge project from OpenLayers2 to Openlayers3. In OL2 this worked and with OL3 not.

This is a code snipped for one layer. The other layers are created in a similar way:

var sourceA = new ol.source.Vector({
 loader: fetchData,
 strategy: ol.loadingstrategy.bbox,
});
if (sourceA) {
 var layerA = new ol.layer.Vector({
  source: sourceA,
  visible: false
 });
 if (layerA) {
  layerA.setZIndex(ZINDEX_LAYER_A);
  ...
  map.addLayer(layerA);
  ...
 }
}


texts:
There is a very similar problem to mine, but with a different plugin. The error is the same.
The error is this:
RuntimeError: wrapped C/C++ object of type RasterizeStringAlgorithm has been deleted

Traceback (most recent call last):
  File "C:/Users/belli/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\rasterizestring\rasterize_string_provider.py", line 60, in loadAlgorithms
    self.addAlgorithm( alg )
RuntimeError: wrapped C/C++ object of type RasterizeStringAlgorithm has been deleted

Versão do Python: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
Versão do QGIS: 3.10.11-A Coruña A Coruña, d2171173e4.

I am not an expert in programming language (IT).


texts:
I'm trying to do some habitat suitability modelling on QGIS. I've downloaded several vector maps from Ordinance Survey (OS). So far I've found a Cartographic stylesheet for  "OS VectorMap Local". This seems to work on Qgis. When I downloaded these there were GML style sheets which could be uploaded by going into properties --> style --> load style.(http://hub.qgis.org/projects/quantum-gis/wiki/OS_Styles)

However I also need other elements from a topographic OS MasterMap Layer and I can't find any styles that work for it.

I've found Styled Layer Descriptors (https://www.ordnancesurvey.co.uk/business-and-government/help-and-support/products/styled-layer-descriptors.html). However when I downloaded it it was just text files and ".sld" files. Which I have no idea to use as a style in Qgis.

Any suggestions?

I don't think it should make a difference but I'm using a mac.

Using the comment below I added in a the sld style to the area and go this:

What I was expecting was something like this: (from the VectorMap style)


texts: How to separate raster bands in ModelBuilder?

texts:
I have a simple model that iterates through folders to collect imagery, mosaics them and clips the mosaic to a study area boundary.  I need to produce an NDVI from the clipped 4-band mosaic "sa11_4band.tif" by incorporating the raster calculator into the model.  The problem I am encountering is that the individual bands do not appear in the raster calculator when I incorporate the raster calculator tool.  How should I separate the raster bands from "sa11_4band.tif" in order to perform the necessary calculation in the raster calculator?:

Float("sa11_4band.tif - Band_4"-"sa11_4band.tif - Band_3")/Float("sa11_4band.tif - Band_4" + "sa11_4band.tif - Band_3")


texts:
I have a ~1GB shapefile which is running circles around a ~5MB shapefile when it comes to feature look-up time (using Python, with GDAL):

def check(name, lat, lon, attributeFilter:'SQL Query'=None):
    shp = openShapefile(name)
    layer = shp.layer

    p = ogr.Geometry(ogr.wkbPoint)
    p.SetPoint_2D(0, lon, lat)

    layer.SetSpatialFilter(p)
    if attributeFilter:
        layer.SetAttributeFilter(attributeFilter)

    for feature in layer:
        return shp.getFields(feature)

    return None

>>> timeit("check('urban', 50.469815, 30.409954)", "from locationgenerator import check", number=1000)
20.37787419225606
>>> timeit("check('regions', 50.469815, 30.409954)", "from locationgenerator import check", number=1000)
0.26346330684964414

Basically the code is checking if the shape file contains any feature at the given coordinates.
openShapefile is a helper function that opens & caches the shapefile for the given name.

The 'regions' shapefile is http://www.gadm.org/version2.
The 'urban' shapefile is https://koordinates.com/layer/1285-world-urban-areas-110-million/.

As you can see, the larger and much higher resolution file performs around 100 times better, which is very suprising to me.

Things I have noticed:

it's not the call to getFields (another helper function of mine), replacing it with return True exhibits the same behaviour
if I modify the regions file (e.g. remove some features, merge some features, etc.), the performance drops several orders of magnitude - which leads me to believe this behaviour is due to some arcane pre-calculated optimizations performed on the shapefile.

How can I achieve at least comparable performance for the much simpler shapefile?


texts:
I have raw data in the form of groundwater levels collected regularly over a couple of years. The values are collected in an excel file together with xy coordinates for the water wells (of which there are between 10 and 20) as well as z ground level values.

My goal is to present the water levels (the important data) in a way that is easy to understand and maybe a bit more interesting than simply showing the position of the wells and keep the values in the attribute table. I want the viewer to understand how the water levels vary within in the area.

But since I'm kind of new to GIS and ArcMap I'm not sure what's a good way to do this. Do you have any suggestions?

What I have done so far is turn the wells into vector points based on the xy values in the excel file. The rest of the values can be viewed in the attribute table.

Since the water wells are located at different heights I guess it's pointless to do an interpolation of the water levels.


texts: Calculating amount of landcover lost if sea level rises to 1m and 2m Contour Line using ArcMap

texts:
I'm relatively new to ArcGIS Desktop 10.1 and stuck trying to understand how to calculate the amount of each land cover type lost if sea level was to rise to a 1m and 2m contour threshold.

Essentially I have created a basic map of land cover for an area adjacent to a harbor. I have extracted 0m 1m and 2m Contour lines from a Lidar-derived EDM (Using a =<2m query) and plotted them on the land cover map as seen below.

I then tried using the intersect tool with both the land cover shp and contours shp but a very similar map was produced outlining essentially the same image as before, as seen below, the attributes table has included some more information but i'm not sure what to use in calcualting the amount of area lost if sea level was to rise.


texts: Retrieving average building height in urban areas?

texts:
I am looking for a way to collect data on the height profile of some cities over time. I  don't need to produce a detailed 3D reconstruction of cities, nor do I need to detect and measure individual buildings. I would simply need a coarse proxy of how "vertical" a city is - for instance, average building height, or maximum height.

I have in mind a few possible options but I'm not sure how involved they are and what kind of software can do this:

Look at the difference between a DSM and a DTM - I am still looking into possible sources (ASTER?)
LiDAR data
Stereo imagery (e.g. IKONOS)
Retrieving building height from shadows – I think even an ArcGis extension does this, but this is not automated and would require a building by building process.

I was told that software such as Socet Gxp, Leica LPS/XPro, and Pixel Factory can generate height data from satellite imagery automatically. Has anyone heard of those or know of any open source alternatives?


texts: Create regular grid in QGIS from coordinates in excel file

texts:
I have received an excel file that contains data on grid cells of 500X500m. The geographical information looks like this: 13843145006371500. The first four digits is the code for the municipality and the others are the xy-coordinates for the south west corner of the grid cell. Can I use this information to create a grid in QGIS?

I use QGIS 2.10 on windows 10


texts: How to create a new field with a combination based on another field with arcpy?

texts:
I like to create a new txt which has two columns: Id and GRIDCODE.
This is the original Table:

The Id column of the new txt should have just unique values. I’m able to do this with the following code. But the GRIDCODE -column of the new txt should have the possible combinations of one equal Id, for example:
The new Id field 47973 should have the new GRIDCODE field 7, 10
or the new Id field 47990 should have the new GRIDCODE field 4, 7.

To clarify the question I build an example-table with excel, this should be the result I like to write in the new txt. Do I have to work with an update cursor and if, else? What about the comma, how is it possible to use them without getting a new column but the combination in a column?

I would be very thankful to get a tip how I can go on with my code.

import arcpy
from arcpy import env

env.overwriteOutput = True
env.workspace = r"D:\Users\jul\ers\cities_UA\resultfolder"

inputshp = r"D:\Users\jul\ers\cities_UA\resultfolder\at003l_linz_result.shp"
outputfile = r"D:\Users\jul\ers\cities_UA\resultfolder\linz_test.txt"
f = open (outputfile, 'w')
f.write ("ID,Gridcode,\n")
f.close ()

f = open (outputfile, 'a')
ID=[row[0] for row in arcpy.da.SearchCursor(inputshp, ["Id"])]
uniqueID = set(ID)
for ID in uniqueID:
print ID

f.write((str(ID)) + "," + "\n")
f.close() #working

Grindcode = [row[0] for row in arcpy.da.SearchCursor(inputshp, ["GRIDCODE"])]
???

Sorry, but I think I’m not able to write the code, what is wrong.

import arcpy
import os
from arcpy import env

env.overwriteOutput = True
env.workspace = r"D:\Users\julia\erste_aufg\cities_UA\resultfolder"

DataDict = {}
inputshp = r"D:\Users\julia\erste_aufg\cities_UA\resultfolder\innsbruckauswahl.shp"
outputfile = r"D:\Users\julia\erste_aufg\cities_UA\resultfolder\inns_test.txt"
f = open (outputfile, 'w')
f.write ("ID,Gridcode,\n")
f.close ()

#f = open (outputfile, 'a')

ID=[row[0] for row in arcpy.da.SearchCursor(inputshp, ["Id"])]
GRIDCODE = [row[0] for row in arcpy.da.SearchCursor(inputshp, ["GRIDCODE"])]
uniqueID = set(ID)
GRIDCODE = set(GRIDCODE)
for ID in uniqueID:
    print ID
    if ID in DataDict:

        DataDict[ID].append(GRIDCODE)
    else:
        DataDict[ID] = [GRIDCODE]

    print DataDict

It is looking like this at the moment:


texts: General strategies for debugging model builder

texts:
When coding a standard program in a text editor there are many tools for debugging problems (syntax errors, debuggers, logging information to the console, unit tests, etc.)

What are some ways to do these things in model builder? Is it possible?  Should I just use the Python API instead?

Thanks!


texts: Merge DEM TIFF grids and get average value

texts:
I have two different DEM grids (basically single-band GeoTIFFs) and I want to merge them.

The grids overlap and have different values, and I want to create a new DEM grid that takes the average or maybe maximum value in the overlapping regions.
The DEM grids already have a nodata flag in them.

Currently I have this script through OsGeo4W:

set pathin=F:\mergebox_DEM\IN
set pathout=F:\mergebox_DEM\OUT
set filename=Merge_result

gdalbuildvrt -tr 0.3 0.3 -r bilinear index.vrt %pathin%\*.tif
gdal_translate index.vrt %pathout%\%filename%.tif -tr 0.3 0.3 -co COMPRESS=PACKBITS -co BIGTIFF=YES

Currently, my merged result looks like this (pseudocolor scale added for emphasis):

This is not very good - you can clearly see that the files are simply overlapped. I think I see why that happens - the bilinear resampling is applied after the two files are overlapped?

Am I missing something very basic?

I don't want to use a raster calculator, because in the future I might have to merge a large number of files and I would like to be able to do it directory-wise.


texts: QGIS: Exporting/saving extents of WMS raster layer in a batch process

texts:
In QGIS is it possible to run: Export -> Save As as a batch command?
I have WMS-raster layers from which I'd like to extract certain extents (based on quite a number of features from a vector layer) and save them as GeoTIFF fur further statistical analyses. As it is a WMS service, I cannot use "Clip Raster by Extent" (or at least I don't know how). Saving the data for each vector feature extent manually is quite time-consuming, so I'd like to run it as batch function.


texts:
I performed a selectByExpression on one of my layers and then created a new layer.
layer 1 : co,
layer 2 : co_temp
  co.selectByExpression("SECTION=1")
  co_temp = processing.run("native:saveselectedfeatures", {'INPUT': co, 'OUTPUT': 'memory:'})['OUTPUT']
  co.removeSelection()
  renderer = co.renderer().clone()
  co_temp.setRenderer(renderer)
  co_temp.renderer().filter()
  co_temp.triggerRepaint()
  co.removeSelection()

I would like to be able to create a co_temp caption, but only for selected co_temp expressions. There is a lot of useless information in co that I would like not to retrieve for the co_temp caption.
mapt = QgsLayoutItemMap(layout)
mapt.attemptSetSceneRect(QRectF(20, 20, 80, 80))
mapt.setFrameEnabled(True)
mapt.setLayers([co_temp])
mapt.setCrs(co_temp.crs())
layout.addLayoutItem(mapt)
mapt.setExtent(co_temp.extent().buffered(250))

legend = QgsLayoutItemLegend(layout)
legend.setTitle("Legend")
legend.attemptSetSceneRect(QRectF(120, 20, 80, 80))
legend.setFrameEnabled(True)
legend.setBackgroundColor(QColor(200, 200, 200))
legend.setTitle('')
legend.setLegendFilterByMapEnabled(True)

legend.setResizeToContents(False)

layout.addLayoutItem(legend)
legend.setLinkedMap(mapt)

I tried this but it doesn't work.
co :

co legend:

co_temp :

co_temp legend:


texts: Copyrights disappearing at certain scales on GoogleEarth Pro

texts:
I want to use a map (scale 1:10000) of a city (Montreal). I assembled it from several high quality pictures.

We need to cite copyrights from GoogleEarth Pro. I have read online resources :
permissions page

Despite trying this, when zooming in, copyrights displayed before (e.g. Cnes Spot/image, DigitalGlobe) disappear so I cannot cite anything for this map at this scale.

Does it mean that GoogleEarth is the sole mention I should make? Is it a bug?


texts: How can I choose values from stacked raster images according to location of ROI

texts:
Before I choosed ROI on land surface temperature (LST) images and obtained pixel values of LST. Now I need pixel value of another images which are stacked multi-band image. My main purpose is to obtain pixel values of both images "at the same location". I think that I can use ROI file from LST to choose the pixel values of multiband image. Anyone who guide me to accomplish this purpose using ARCGIS, ENVI, ERDAS?

I choosed the ROI in rectangular mood in ENVI. When I opened it, it looks like small polygons (not points) spreading on my study area.

Thanks a lot


texts: When should I run a PCA analysis to classify satelite images?

texts:
I want to classify (using training data) some LANDSAT 8 OLI images.
I've come across some mentions of PCA analyses to avoid band-to-band correlations.
In what cases should one run a PCA? What is the advantage of doing so?


texts:
I've inherited a pretty large / complex (internal) web portal (run under python 2.6.6), which has some mapping capabilities, using Mapnik 2.2.0.  The web service has recently started throwing 500 errors when it tries to create some of the maps that are displayed, with the sole error:
$ Premature end of script headers: app.wsgi

The maps are generated from geojson that is passed as a string (having been created from a python dictionary); I've tracked it down to 1 line of code in a map.py, which adds a layer to a map:
data = mapnik.Ogr(file=json, layer='OGRGeoJSON')

when I comment out that line of code (and the associated m.layers[-1].data = data) then the image is displayed on the portal, but obviously without this layer.  The json is valid (I have checked it on http://geojsonlint.com).

Debugging efforts
Using a separate json file
I wanted to see whether calling the JSON from an external source replicated the problem.  I therefore saved the json string as a separate file (test.json), in the same directory as map.py, and modified the Ogr call to
data = mapnik.Ogr(file='test.json', layer='OGRGeoJSON')

I still get a 500 error, but with a slightly more verbose error message:
datasource = mapnik.Ogr(file='testgeo.json', layer='OGRGeoJSON')
File /usr/lib64/python2.6/site-packages/mapnik/__init__.py, line 536, in Ogr
return CreateDatasource(keywords)
RuntimeError: OGR Plugin: Failed to read GeoJSON data

Checking the geojson
$ ogrinfo testgeo.json
ERROR 4: GeoJSON Driver doesn't support update.
Had to open data source read-only.
INFO: Open of `testgeo.json'
      using driver `GeoJSON' successful.
1: OGRGeoJSON (3D Polygon)

So the json seems valid

Using the python command line
I have also tried running python interactively in the same directory as map.py and test.json:
$ python
Python 2.6.6 (r266:84292, Jul 23 2015, 15:22:56)
[GCC 4.4.7 20120313 (Red Hat 4.4.7-11)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import mapnik
>>> mapnik.Ogr(file="testgeo.json", layer="OGRGeoJSON")
<mapnik.Datasource object at 0x7f6243c27398>

Further info
$ gdal-config --version
1.9.2

$ ogr2ogr --version
GDAL 1.9.2, released 2012/10/08

This returns with no errors (so I assume that this means the geojson is valid, confirming my earlier test of running the geojson through http://geojsonlint.com):
$ ogr2ogr testgeo.shp testgeo.json

OS Version
$ lsb_release -a
LSB Version:    :base-4.0-amd64:base-4.0-noarch:core-4.0-amd64:core-4.0-noarch:graphics-4.0-amd64:graphics-4.0-noarch:printing-4.0-amd64:printing-4.0-noarch
Distributor ID: CentOS
Description:    CentOS release 6.7 (Final)
Release:        6.7
Codename:       Final

These 3 debugging efforts show that the json is valid and that the call works, except when it's run through the webstack. I've checked the permissions and ownership of all relevant files and directories, and they all seem to be as I'd expect.
What else could be causing the issue?!

Update (2 weeks after asking the original question)
I have come into work today and thought I'd have another crack at looking to see if I could tease some error messages out (new week, fresh start).  But, utterly bizarrely, the code does not seem to be failing (for now, at least).  I have changed nothing, and as far as I can gather, there have been no hardware changes to the servers.
This now becomes an investigation of what did cause it, and whether there's anything I can do to prevent it happening again...


texts: Seeking Sentinel L1C images (.SAFE) from 2015-2019 since no longer available on Copernicus Open Access Hub

texts:
I am looking for Sentinel images from before 2019 in L1C quality and .SAFE format, as I want to correct them with sen2cor in SNAP to convert them to L2A quality.
Where can I find them now that they are offline & not available on the Copernicus Open Access Hub?
The Hub took them offline, and they are currently not requestable.
USGS Earthexplorer and Google Earth Engine only provide TIFFs and JPEGs.
So far I have only acquired JSONs from Amazon Web Services (connected to Sentinel Hub), but learning Python might make me able to acquire .SAFEs ..
Details:
I am making a monthly NDVI timeseries from december 2015 - 2020 of my research location in Bolivia. I use Sentinel 2 L2A images for this, and Google Earth Engine. Because there are no L2A images before 2019 I will have to process L1C images myself and upload them to GEE.


texts: Batch loading QgsRasterLayer within a QgsTask

texts:
When I try to load many huge rasterlayer into QGIS it blocks the UI. So I try to implement a QGIS-batch-raster-loading task in order to keep QGIS responsive.
Calling the task itself from within my plugin works so far but I cannot add layers to my project like this:
project = QgsProject.instance()
project.addMapLayers([my_array])

In the QGIS docs this warning is written:

Any background task (regardless of how it is created) must NEVER use
any QObject that lives on the main thread, such as accessing
QgsVectorLayer, QgsProject or perform any GUI based operations like
creating new widgets or interacting with existing widgets. Qt widgets
must only be accessed or modified from the main thread. Data that is
used in a task must be copied before the task is started. Attempting
to use them from background threads will result in crashes.

I wonder if this means that it is just not possible to bulk add layers from within a QgsTask or if I miss something?
The main snippet of my code is:
class RandomIntegerSumTask(QgsTask):
    """This shows how to subclass QgsTask"""

    def __init__(self, description, duration, layer):
        super().__init__(description, QgsTask.CanCancel)
        self.duration = duration
        self.total = 0
        self.iterations = 0
        self.exception = None
        self.layer=layer

    def addRasterLayer(self,layers_to_add, project):
        layerarray=[]
        mycrs = QgsCoordinateReferenceSystem.fromEpsgId(25832)
        for path,name in layers_to_add.items():
            print (path, name)
            rlayer = QgsRasterLayer(path, name)
            if rlayer is None:
                continue
            rlayer.setCrs(mycrs,True)
            layerarray.append(rlayer)

            project.addMapLayers([rlayer])

    def run(self):
        """Here you implement your heavy lifting.
        Should periodically test for isCanceled() to gracefully
        abort.
        This method MUST return True or False.
        Raising exceptions will crash QGIS, so we handle them
        internally and raise them in self.finished
        """
        try:
            project = QgsProject.instance()
            self.layers_to_add={}
            for path in self.layer:
                name = QFileInfo(path).baseName()
                self.layers_to_add[path]=name

            self.addrasterlayer(self.layers_to_add,project)


texts: Accessing additional CAD annotation fields

texts:
I'm trying to read certain values from the annotation layer of a CAD file (MicroStation.dgn).  It seems like for certain entities (Cell) the text isn't stored in the field 'TxtMemo' or 'Text'.  When I use the Import CAD Annotation tool I see additional fields including 'TextString' which has the text value I'm trying to obtain.

Is there anyway to access this field directly using the SearchCursor tool directly on the CAD annotation features?

arcpy.da.SearchCursor(cad_annotations, 'TextString')

results in the error:

RuntimeError: A column was specified that does not exist.

I guess the workaround is to add additional code to do the import to a temporary location and read the data there but given that it can display the 'TextString' value when I add the CAD annotation layer to a map it seems like there should be some way to access it directly without the conversion.  I'm just wondering if I'm missing something.


texts: Native Mapping Application - Leaflet Mapbox

texts:
I am going to build a native mapping application that I would like to offer for Android and IOS devices. I am an accomplished developer of web-based app using javascript, html and css. I am also familiar with node.js. I would like to continue using my existing stack and use Leaflet and Mapbox for the mapping components.

Can I develop my app as I always have using javascript, html and css? What would I then need to do to make my app available in the Apple and Google app stores?

Is my app going to be limited by taking this approach rather than using Mapbox's SDKs?


texts: how can i save a satellite raster image as a layer file in arcmap 10.2

texts: I am trying to save a raster satellite image from esri world imagery as a georeferenced layer file. i still want to be able to zoom in and out as i can when veiwing it through the server.


texts: Extract a polyline from a Geoserver shapefile

texts:
I have a GeoServer which contains a shapefile with some roads on it. Is it possible to programmatically query the GeoServer and extract a specific road from it and display it in a layer? I'm using OpenLayers library locally.

The reason I want to do this is because I want my line to be snapped on to the actual roads on my basemap (Bing maps).


texts: file type needed to add layer containing contour lines

texts:
What file format is usable/recognized by qgis to add contour lines to a map?

I'm working with Linux Mint 15, and qgis 1.7.5


texts:
I have a database with fields including latitude, longitude, city, state, and country.

When mapped, it is clear that a number of the sites are in the wrong location - most commonly because of a sign error:

I was thinking of writing a script to check for such errors, along the lines of

if country = USA and lon > 0 set lon = -1*lon
if country in (select names from europe) and lat < 0 set lat = -1*lat

etc.

This would be followed by additional checks by plotting the sites with their names over google maps, but this would get most of the errors.

What pitfalls to this approach should I anticipate?
Any suggestions for how to do comprehensively way?
Are there existing scripts that do this? It is hard to imagine that I am the first to face this issue.


texts: How to activate the cursor once we deactivated by activating panning tool

texts:
We have activated pan tool on canvas using the following command:

canvasPressEvent(self,event):
   if event.button == 1:
       qgis.utils.iface.actionPan().trigger()
 canvasReleaseEvent(self,event):
   if event.button == 2:
      print "hello"

the problem is that it is not at all going to the canvasrelease event. So i want the cursor to be active to draw a line when right click is pressed after activating pan tool how can i do that using python programmin?


texts: cannot delete row from coverage with updatecursor

texts:
I am able to update without any issue or error but when I try to delete a row I'm getting an error, see below for the code and error

I can't see why it would fail.

import arcpy

def count_row(data):
    i = 0
    for row in data:
        i += 1

    data.reset()

    return i

arcpy.env.workspace = 'C:/Data/cea01apr'
feature = 'point'
fidfilter= 'fid = 1'

featureClassList = arcpy.ListFeatureClasses()

fieldList = arcpy.ListFields(feature)
fieldString = [str(x.name) for x in fieldList]

with arcpy.da.UpdateCursor(feature, fieldString, fidfilter) as cursorupdate:
    print('Number of rows:' + str(count_row(cursorupdate)))
    rowupdate = cursorupdate.next()
    print(rowupdate[12])

    rowupdate[12] +=  'h'
    cursorupdate.updateRow(rowupdate)

with arcpy.da.SearchCursor(feature, fieldString, fidfilter) as cursorupdate:
    print('Number of rows:' + str(count_row(cursorupdate)))
    rowupdate = cursorupdate.next()

    if rowupdate:
        print(rowupdate[12])

with arcpy.da.UpdateCursor(feature, fieldString, fidfilter) as cursorupdate1:
    print('Number of rows:' + str(count_row(cursorupdate1)))
    for row in cursorupdate1:
        cursorupdate1.deleteRow()

raw_input('Press the <ENTER> key to exit')

the output is

C:\Data>test.py
Number of rows:1
 hhhhhhhhhhhhhhhh
Number of rows:1
 hhhhhhhhhhhhhhhhh
Number of rows:1
Traceback (most recent call last):
  File "C:\Data\test.py", line 39, in <module>
    cursorupdate1.deleteRow()
SystemError: error return without exception set

I have tried this

arcpy.MakeTableView_management(feature, 'test')
arcpy.SelectLayerByAttribute_management('test', 'NEW_SELECTION', fidfilter)
arcpy.DeleteRows_management('test')

but got this error

Traceback (most recent call last):
  File "C:\Data\test.py", line 41, in <module>
    arcpy.DeleteRows_management('test')
  File "C:\***\cots\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 15352, in DeleteRows
    raise e
arcgisscripting.ExecuteError: ERROR 999999: Error executing function.
Failed to execute (DeleteRows).


texts: Performance of GeoServer for thousands of spatial data files

texts:
Input:

thousands of large files (from 1mb to 5gb) containing the coordinates of the points (x, y, z) in a specific coordinate system. Points can intersects

It is required to make an online service where on top of the map there will be a colored surface (like a heatmap depending on the Z coordinate) with some interactive functions. More precisely 3 surfaces with different color maps and you can switch between them

Now GeoServer (with GeoWebCache), GDAL utilities, and some web-stack with leaflet+OSM at the front are used.

How it works now:

Preparing files using the GDAL utilities.

'og2ogr' to Shapefile from CSV (with srs changing)
'gdal_grid -zfield Z' from Shapefile to Tiff
3 times 'gdaldem color-relief ... map.txt ... -alpha' for coloring Tiff and 'gdal_translate -of PNG -scale -co worldfile = yes ...' for creating PNG Worldfile.
Output 3 PNG heatmaps in different colors.

Upload PNG to GeoServer. We have 3 groups of layers by color. Adding a new layer to a group of layers and recalculating the boundaries of the group.
On the front by 'L.tileLayer.wms' we show the desired group of layers from GeoServer.

Problems and question:

The main problem is the performance. For the quickness on the client, I run Seed in the GWC of each group, but it works very slowly. (for example, 1 group, 300 layers in a group, zoom up to 17 can go about a day). But after the cache is ready, everything works quickly on the client. But as boundaries grow, productivity gets worse.

This is my first GIS project and I do not have enough knowledge. It seems to me that my decision is not good and not optimal. What optimal approaches and tools can be used to solve such a problem, and what parameters/options can speed it up?


texts: How do I find the percent overlap between two shapefiles in QGIS?

texts:
I have two different shapefiles (first and second). I want to calculate how much overlap there is between the liquefaction zones in the second shapefile and each of the neighborhoods in the first.

I'm using QGIS. My ideal output would be something like:

nhood    percent in liquefaction
Glen Park            .75
Financial District   .60
Tenderloin           .32
....                 ....


texts: Creating a map book with maps and layouts of different sizes in one mxd

texts:
I'm looking for a way to create a map book containing maps of different sizes using Data Driven Pages. Differently sized index features (1, 2, 3 and 4 in the sketch) should run along a line (in this case a railway). The result should be a series of maps with different sized layouts, for example:

page 1 (height/length): 30 cm x 50 cm / 12" x 20 "
page 2: 30 cm x 70 cm / 12" x 28 "
...

If it were only two pages, I would create two mxd files, but for 40 different sized maps it would take too long.

I'm using ArcGIS Desktop Basic, Version 10.6


texts: Leaflet: tool for multiple marker selection and computation of summary

texts: Let's suppose that we've got a leaflet map with some markers and another layer, a raster image let's say. Is there a tool, or a way to create one, that allows one to select with an area, a subset of those markers? And then, how could one calculate a summary of the selected group's corresponding variables? I am looking for a selection tool simmilar to the measuring tool that leaflet already supports, and then a way to calculate the summary and output it to a panel in the map. So far I have looked at the Crosstalk package in R, but I have not seen a way to use the selected values to calculate an aggregate statistic. Also, googling didn't clarify much the situation for me, and there doesn't seem to be a clear answer on the git hub threads I've read.


texts: Create 3D lines at a specific distance and bearing from a 2D polygon using Shapely and Geopandas

texts:
I have a shapefile of the runways(polygons) of different airports:

I want to build a line that represents the 3D trajectory followed by an airplane 9 km after take off and 18 km before landing and 1000 meters above ground.
My questions:

How to calculate the coordinates of the starting points?
How to calculate the coordinates of the target 3D points?
How to use shapely to create the 3D Multilines

I'm trying to solve this as follows:

Creating a line that goes through the polygon exactly in the middle, obtain the coordinates of the vertices of this line and supply them as the coordinates of the known point, making a lineal regression it will achieved it.

Calculate the coordinates of a point 9 kilometers away from the runway and keeping the orientation of the runway, that is, it will be a 3D line projected from the runway, the same to represent the landing.

Use the equation of the direct problem in geodesy, to calculate the coordinates of the points 9km and 18km away from the runway

The output I would like to achieve, is like that:


texts: Google Maps satellite view footprints not congruent to street map view

texts:
I am trying to draw the layout of a building and project them to a web map using Django and Leaflet. I am using a satellite image to get the correct latitude and longitude and then I send the coordinates via an API. The process in short:

Get lat/lon from a satellite image
Send coordinates to webmap with SRID 4326
Display the outlines

My main problem with almost all my coordinates is that the Google Maps satellite layer is not congruent to the street view. For example let's say I pinpoint lat: 50.11843503887959 and long 8.63194562490686 to Google Maps. I get the following image in satellite view:

Note that the point is exactly on the edge. But if I change the layer to the normal map view I get this result:

Note that the point is about 1cm to the right within the building. This happens to almost all buildings I am trying to pinpoint. Sometimes it's to the left and sometimes it's to the right.
I read here EPSG 3857 or 4326 for GoogleMaps, OpenStreetMap and Leaflet, that:

Google Earth is in a Geographic coordinate system with the WGS84 datum. (EPSG: 4326)
Google Maps is in a projected coordinate system that is based on the WGS84 datum. (EPSG 3857)

Thus, I went ahead and transformed the satellite coordinates EPSG: 4326 to EPSG 3857 and sent those to my database. It gives me the same derivations.
I assume that since the Earth is round and a map is flat those derivations exist, but I am not sure.
Plus I somehow need to iron out these differences. Is there a way I can get coordinates of a building from a satellite layer to map the outlines of the same building on the map layer? If so, how could I do it?
(Had this on SO (https://stackoverflow.com/questions/66366572/) first but removing it there since it might be better suited here)


texts: Grouping lat/lng into grid of 1 km squares

texts:
We have a lot of lat/lng that we would like to group into a grid of ~1 km squares.

1) Grouping lat isn't difficult. That's just grouping them by lat +/- 0.02, i.e. 14.X will be grouped like this:

14.00-14.02
14.02-14.04
14.04-14.06
etc.

2) Problem is lng, since the length of each lng varies depending on where you are looking.

Does anyone know a formula for converting lng to km that handles the varying length of lng?

I know this is not 100% precise 1 km squares, for it's good enough for our purpose.


texts:
I'm trying to map data that I uploaded with an csv file. Originally this data set only included the country code, e.g. BD, for Bangladesh. I tried to geocode it by admin region  and selected the column where the country code is stored. It gives me an error message.

Then I updated the data file with more specific information on the location: address, city and country. But this also did not work.

What do I do wrong? Should try it with long and lat values?

Is there a basic tutorial for newbees like me?


texts: Can subtypes be re-ordered or a new subtype inserted in a specific position?

texts:
I have a feature class with sub-types that need to be modified and re-ordered as new data becomes available.  Right now I have 20 subtypes but could potentially have close to a 100. I use the subtypes numbers for labeling so the subtypes and the subtypes need to be in a logical order, grouped by similar items.

Having to insert a new item in position 3 for example, requires me to rebuild the entire list from position 3 on; adding a new item in position for and bumping up the remaining 17 items.
http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/an-overview-of-the-subtypes-toolset.htm
This will get really cumbersome really fast.  (This features class is exported to a feature service and accessed by many users on our web GIS.)
I don't see a way to add then re-order subtypes without having to rewrite them all, or insert a new subtype at a specific position, or export subtypes to a file, modify this list externally (notepad) then import subtypes.

Is there a workaround or some other method I am not seeing that would make it easier to change the order of subtypes?


texts: Is it acceptable to store OSGB36 data as geometry in SQL Server?

texts:
I've set myself a project of getting some geo data into SQL Server but before I get started doing any of the big jobs I wanted to see if I was on the right track.

All of my data is for the UK and is therefore based on Easting & Northing coordinates (OSGB36), not Latitude & Longitude (WGS84).

My understanding is that I cannot create a geography field based on Eastings & Northings in SQL Server without first re-projecting in a GIS package (I have QGIS).

I can however convert Easting & Northings based data into a geometry field with no problems.

I have managed to re-project some point data to WGS84 in QGIS and create a geography field from this, however when I use the QGIS MSSQL link to bring the resulting point data into QGIS, it really doesn't seem to be able to marry it with layers in OSGB36 (despite having on-the-fly re-projection activated).

So my question is: Following the path of least resistance (and maximising compatibility with other UK datasets), is it acceptable to store OSGB36 data as geometry in SQL Server, despite actually being geography?


texts:
I am working with QGIS (newbie). I have a raster map with 0.5degrees resolution. I created a 0.5degrees width grid layer and I compute zonal statistics from the raster.
My problem is that it seems to me that the computed zonal statistics are shifted 0.5degrees, it is, they do not correspond to the actual zones in the raster layer.
I think that posting some images will make things clearer: http://imgur.com/a/iaUkx
* The grey-scale map is the raster file and the blue-colored map is the zonal statistics displayed as a gradient. As you can see, the blue tiles are shifted one tile to the left.

How can I solve this problem?


texts: How to use sf:gdal_utils with extent option (-te)?

texts:
I want to build a Gdal VRT dataset with several bands and cropped to a specific region. For that, I'm using sf::gdal_utils but somehow I can't get the option -te to work. It either complains about an unknown option or does nothing.
How to correctly set the -te option?
Here goes a reprex:
# bug in gdal_utils and buidvrt
# trying to use sf::gdal_utils buildvrt with -te option

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(raster)
#> Loading required package: sp

# just to check extent
#tm <- raster(system.file("tif/L7_ETMs.tif", package="stars"))

# crop file using gdal VRT -te option

# Unknown option -te
# also tried with options te, using just one string instead of c()...
# can't figure out how to use this
gdal_utils(util = 'buildvrt', source = system.file("tif/L7_ETMs.tif", package="stars"),
           destination = 'tm_crop_r.vrt',
           options = c("-te", "290000 9119000 295000 9120000"))
#> Warning in CPL_gdalbuildvrt(source, destination, options, oo, quiet): GDAL Error
#> 6: Unknown option name '-te'

# check the size of the VRT dataset
readLines('tm_crop_r.vrt', n = 1)
#> [1] "<VRTDataset rasterXSize=\"349\" rasterYSize=\"352\">"

# using gdal CLI works

system(paste("gdalbuildvrt -te 290000 9119000 295000 9120000",
             "tm_crop_cli.vrt /home/daniel/R/x86_64-pc-linux-gnu-library/4.0/stars/tif/L7_ETMs.tif"))

readLines('tm_crop_cli.vrt', n = 1)
#> [1] "<VRTDataset rasterXSize=\"175\" rasterYSize=\"35\">"

sessionInfo()
#> R version 4.0.2 (2020-06-22)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.5 LTS
#>
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
#>
#> locale:
#>  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C
#>  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8
#>  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8
#>  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base
#>
#> other attached packages:
#> [1] raster_3.3-13 sp_1.4-2      sf_0.9-6
#>
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.5         knitr_1.29         magrittr_1.5       units_0.6-7
#>  [5] tidyselect_1.1.0   lattice_0.20-41    R6_2.4.1           rlang_0.4.7
#>  [9] stringr_1.4.0      highr_0.8          dplyr_1.0.2        tools_4.0.2
#> [13] grid_4.0.2         xfun_0.16          KernSmooth_2.23-17 e1071_1.7-3
#> [17] DBI_1.1.0          ellipsis_0.3.1     htmltools_0.5.0    class_7.3-17
#> [21] yaml_2.2.1         digest_0.6.25      tibble_3.0.3       lifecycle_0.2.0
#> [25] crayon_1.3.4       purrr_0.3.4        codetools_0.2-16   vctrs_0.3.2
#> [29] glue_1.4.1         evaluate_0.14      rmarkdown_2.3      stringi_1.4.6
#> [33] pillar_1.4.6       compiler_4.0.2     generics_0.0.2     classInt_0.4-3
#> [37] pkgconfig_2.0.3

Created on 2020-09-18 by the reprex package (v0.3.0)


texts: Select characters using regular expression in QGIS

texts:
I'm trying to seperate characters from database to seperate fields. Database is like:

999-10-1-1
999-10-11-111
999-10-111-1111
999-10-22-1
999-10-1111-2221

and I want something like in seperate fields:

999 10 1 1
999 10 11 111
999 10 111 1111
999 10 22 1
999 10 1111 2221

Basically I want to pick numbers that are seperated by the (-) character.


texts:
How do I write a SLD file based on more than one Value field for each Name. I would like to symbolize by two fields "Type" and "Aabningsgrad".
Some features I would like to see with the combination Type = Ventil and Aabningsgrad = 0, and others I would like to show with combination Type = Ventil and Aabningsgrad = 100.

How do I write that?

I looked at this:
SLD based on two attribute values
But unfortunately it didn't do the trick for me.


texts:
I'm trying to rename a group of rasters to "I1986_" + [The 1st 3 letters of the raster]. For example, if a raster is named "Raster1", the output would be "I1986_Ras".

I've tried including str() functions in a couple places as well as .join in the tool's out data parameter but end up generating other errors.  If I don't use rList variable in the Rename function and use the rasters variable instead, the script will rename the first file in the directory but then say i1986_xxx already exists.  Any ideas?

arcpy.env.workspace = DIRECTORY
rList = arcpy.ListRasters()
for rasters in rList:
    arcpy.Rename_management(rList, "I1986_" + rList[0:3])


texts:
I am trying to view my DEM of terrain for a project in the 3D viewfinder, but having some difficulty. I can see it as flat terrain, but when I turn on DEM for the interpolation it disappears in the viewfinder.
Does it have to be on a certain projection?
I am running the most current version of QGIS on my machine.


texts: QGIS find and replace several unicode symbols in a single field, including replacing an apostrophe

texts:
I work with languages that use Unicode symbols in their placenames. I often receive layers with labels that display incorrectly, but they are messed up in the same way all the time eg: the place name shows Ñ but it should read K̲, or ñ should be k̲.
I have the code for a single find and replace, but I would like to batch find and replace all the problem letters at once. How do I string together several letter replacements?
This is the code I have to replace a single letter:
replace("field",'find_value','replace_value')

But I would like it to fix all of these values at once:
  replace("Label",'Ñ','K̲')
  replace("Label",'Ý','X̲')
  replace("Label",'Ç','G̲')
  replace("Label",'ç','ɢ̲')

I also need to replace a regular apostrophe ' (U+0027) in the middle of the word word with a right single quotation mark (U+2019) style apostrophe. This is my code that doesn't work:
replace("Label",''','’')


texts: How to extract second largest value via Field Calculator?

texts:
I am trying to find the largest and the second largest value in a row from my shapefield table. The largest is easy to find with the max() expression. The second largest is quite easy to find in Python via:

a = [VALUE_1,VALUE_2,VALUE_3,VALUE_4,VALUE_5]
a.sort()
print a[-2]

If I use the same code in the field calculator I get the message that VALUE_1 is not defined. But I stated the input in the expression.


texts: Assigning features to spatial clusters created by Oracle SDO_SAM.SPATIAL_CLUSTERS

texts:
SDO_SAM.SPATIAL_CLUSTERS creates spatial clusters - how to detect which feature was assigned to which cluster? Is the only way to perform an additional spatial join between cluster geometries and original features?
For comparison: Spatial clustering functions in PostGIS like ST_ClusterKMeans return a cluster ID for each feature, not the cluster geometry.


texts: Setting fill style of vector layer in OpenLayers 3?

texts:
I am trying to apply a pre-defined fill-style for my vector layer but I face an error. The firebug says the following:

SyntaxError: missing } after property list

style:myStyle

I tried to put a "}" afther the property list, but the debugger keeps saying the same. I also tried to check my Javascript code with JSHint, but it also says the same problem. Last but not least I tried to look for some solution in this site and stackoverflow too, but I could not find any for my case :\
I am of the opinion that I put the "style: myStyle" line on a wrong place...

Can anyone help me?

My source code:

<!doctype html>

    <link rel="stylesheet" href="http://openlayers.org/en/v3.10.1/css/ol.css" type="text/css">

    <title>OpenLayers 3 example</title>

My Map


texts: Understanding concept of ArcGIS for Flex?

texts:
I am a complete newbie to GIS.

I am trying to download topographic data from a website that displays a map that can be zoomed in and out.  The website is setup with:

  This application is using ArcGIS API for Flex version 2.3.1.  Build
  date 2011-04-26.

I don't understand how the topo data is setup on the server, is this vector data or raster? I can zoom to the most detailed level and get individual 256x256 .jpg tiles but that can't be how the data is setup on the server because that would mean hundreds or thousands of tiles for each zoom layer. I can get very large .gdb files but they don't seem to contain the correct set of files because they won't open in ArcMap.


texts: Why isn't QGIS detecting CRS from .prj file?

texts:
I have a number of 1 km hexagonal grids that cover various counties in the United States in a postgreSQL/postGIS database. Each grid has the CRS EPSG:3857, and the counties layer has EPSG:3857. When viewing the grids with the counties in QGIS, everything looks grand.

But ... in order to share these grids with colleagues, I've had to export them to shapefiles using ogr2ogr. Viewing these in QGIS, each grid looks nudged up about 20 km or so, and QGIS automatically sets the CRS to EPSG:3395 (which is not the project CRS).

When I export the postGIS tables as shapefiles from QGIS, the .prj file looks exactly the same as the ogr2ogr exported shapefiles, but the postGIS exported tables are displayed correctly. I noticed that QGIS creates a .qpj file when exporting shapefiles from QGIS, so I've come to the conclusion that QGIS is ignoring the .prj and looking for a .qpj instead. Why can't it read the .prj without a .qpj? Other shapefiles (such as those from the US Census) don't have a .qpj but QGIS displays these correctly.

I've come up with a workaround by saving a default.qpj and creating a new .qpj from this for every file that export using ogr2ogr, but this seems messy and obviously not reproducible as it only works for EPSG:3857.

Sidenote: I'm using QGIS 2.0.1.

EDIT:

Here's the ogr2ogr command I used:

ogr2ogr -f "ESRI Shapefile" /home/matt/data/hex_grid_1 PG:'dbname=mydb user=matt' hex_grid_1

Contents of the .prj:

  PROJCS["WGS_84_Pseudo_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],PARAMETER["standard_parallel_1",0.0]]

Contents of the .qpj:

  PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","3857"]]

EDIT:

The problem was solved by converting the EPSG:3857's to EPSG:2163 in all my scripts. I'm still not sure what the problem is though since the grids displayed correctly in QGIS when originally loaded from a postgreSQL table (with EPSG:3857).

My workaround proved crude as I thought it did, since my colleague was unable to use the file in ArcGIS, which didn't read the .prj or the .qpj properly.


texts: Agent Analyst tool interface will not open on ArcGIS for Desktop?

texts:
For me the Agent Analyst tool interface won't open.

I have ArcGIS 10.1 sp1 installed on my computer.

I have followed these steps so far:

added the refresh.config.exe text file (ArcGIS 10.1 using Repast's Agent Analyst Extension - refresh.exe does not work)
changed values 10.0 to 10.1 in files agent_analyst.bat and agent_analyst_run.bat in agent analyst installation folder
3- added the ArcGIS 10.1 installation path "for example C:\Program Files (x86)\ArcGIS\Desktop10.1" to environment variables ARC_JAVA & ARC_JAVAW & ARCOBJECTS_JAR (Which GIS software or plugins can perform a Cellular Automata Simulation?)

I restarted and still have had no luck (just a black box flashes up when I create a new agent analyst tool or try to edit the tool).

Is there some other step I need to do or have I done one of the above steps incorrectly?


texts: Formatting label expression with newline using Python parser?

texts:
I'm having trouble with syntax in 10.2.1, and have have tried just about everything.  With arcpy, I am attempting to access a label.expression of an ArcMap layer and modify it to:

lbl.expression = '[FIELDNAME] + "\n" + "taco"'

I can't get the "\n" part to force the rest of the expression to a new line.


texts: Zoom on a PostGIS vector layer in a new project with PyQGIS?

texts:
I load into a PostGIS database, 2 shapefiles of the QGIS sample data :

alaska.shp
storagep.shp

with the shell code :

shp2pgsql -s 2964 storagep.shp public.storagep > storagep.sql
shp2pgsql -s 2964 alaska.shp public.alaska > alaska.sql
psql -h localhost -U postgres -d postgres -f storagep.sql
psql -h localhost -U postgres -d postgres -f alaska.sql

Under QGIS 3.6.1, I execute in a new project the following code in the QGIS Python editor :

def loadPgLayer(instance, vl, path, name, data_filter=""):
    """load a postgis vector layer"""
    vl.setDataSource(path, name, "postgres")
    if vl.isValid() and vl.setSubsetString(data_filter):
        instance.addMapLayer(vl)

canvas = iface.mapCanvas()
instance = QgsProject.instance()

uri = QgsDataSourceUri()
uri_params = dict(aHost="127.0.0.1",
                  aPort="5432",
                  aDatabase="postgres",
                  aUsername="postgres",
                  aPassword="")
uri.setConnection(**uri_params)

vl_alaska = QgsVectorLayer()
vl_storagep = QgsVectorLayer()

name = "alaska"
uri.setDataSource("public", name, "geom")
path = uri.uri(False)
loadPgLayer(instance, vl_alaska, path, name)

name = "storagep"
uri.setDataSource("public", name, "geom")
path = uri.uri(False)
loadPgLayer(instance, vl_storagep, path, name)

# set the project crs
instance.setCrs(QgsCoordinateReferenceSystem(2964))
# zoom on vl_storagep extent
rect = canvas.mapSettings().layerExtentToOutputExtent(vl_storagep, vl_storagep.extent())
canvas.setExtent(rect)
# for update view
canvas.refresh()

This code load the two PostGIS vector layers but doesn't zoom on the extents of vl_storagep layer.

But, after executing this code above, in the QGIS Python console I enter this piece of code below again :

canvas.setExtent(rect)
# for update view
canvas.refresh()

and it zooms on the desired extents.

If I try this kind of code with Shapefiles Vector Layers, everything is good, but as soon as it's PostGIS Vector Layers, this malfunction (?) appears.

For information :

OS : Ubuntu 18.04
QGIS : 3.6.1
PostGIS :

SELECT postgis_version();
--> 2.4 USE_GEOS=1 USE_PROJ=1 USE_STATS=1


texts: Areal Interpolation - Error 40066 - cannot distinguish between polygons

texts:
I have data for some countries but not for others and would like to predict the missing values from neighboring countries weighted by their area, distance etc.

To do so, I am using the geostatistical wizard to interpolate information between polygons. Since I have not even succeeded in implementing areal interpolation using the ArcGIS interface, I did not dare to use ArcPy which is why I cannot provide you with any code. Here is what I did:

I downloaded the world country geodatabase available from ESRI (https://www.arcgis.com/home/item.html?id=3864c63872d84aec91933618e3815dd2), merged it with my data and then used the tool Areal Interpolation from the geostatistical wizard.

At step 2, I received the error 40066: cannot distinguish between polygons which presumably occurs because some polygons (islands) are much smaller than others. As suggested in the ESRI help, I used the tool Delete Identical to identify and delete those polygons that ArcGIS could not distinguish between. Apparently, this was successful but when I tried to continue with the Areal Interpolation I got the same error 40066 again.

Did anyone have similar problems?


texts:
I can't install any plugins on QGIS since updating to 2.18 (I am running High Sierra). I am connected to the QGIS Plugin Repository yet receive the error message:

  Plugin installation failed:Failed to unzip the plugin package.
  Probably it's broken or missing from the repository. You may also want
  to make sure that you have write permission to the plugin directory.

I have downloaded the plugin (openlayers) and extracted to the .qgis2/python/plugins/ directory but receive the error message: No QGIS plugins found. I am not using a proxy so don't think that is the error.

Any suggestions?


texts: PostGIS get geometries within a radius of n meters using WGS84

texts:
select *
from zones z
where st_intersects(z.geom, st_Buffer(ST_SetSRID(ST_Point(2.336031, 48.863172), 3857),1));

It looks like the 1 in the st_buffer is being interpreted as degrees, I need to find geometries within n meters of a given point.

I did find a reference to a new function (st_Buffer_Meters) but that didn't seem to do anything different.


texts: Changing qgis2leaf custom popup width and height?

texts:
I'm using qgis2leaf on QGIS 2.6.

I'm actually working on a map. I need biggers popup. I succeed to change the popup width and height in the leaflet.css file. But the text inside do not fill the whole popup, but just a tiny part on the left. I'm afraid I'm not clear, so you can see my problem using the link below (on the map, clik for example on "THERMES DE CARACALLA").
http://jbouffand.free.fr/spip.php?article320

Can someone help me?


texts: Why is arcpy.da.SearchCursor code giving AttributeError: 'int' object has no attribute 'split'?

texts:
I am new to python and am trying to get a string field in a geodatabase table to split and create a new row. Here is the code based on the solution from Splitting strings into new rows, retaining geometry:

import arcpy

inputTable="PathtoTable/XYZ.mdb/EasementsScratchSelectedRecords"
inputTableAgain="PathtoTable/XYZ.mdb/EasementsScratchSelectedRecords"

searchC =  arcpy.da.SearchCursor(inputTable, ("OBJECTID","DocumentNameorNumber","SecondaryIdentifier"))

insertC = arcpy.da.InsertCursor(inputTableAgain, ("OBJECTID","DocumentNameorNumber","SecondaryIdentifier"))

for row in searchC:
    A = [x.split('Point') for x in row[:-1] ]

    for i in range(min(len(A[0]),len(A[1]),len(A[2]))):
        newRow = [A[0][i] , A[1][i], A[2][i], row[3]]
        insertC.insertRow(newRow)
del insertC

I keep getting the following error:

line 17, in
    A = [str(x.split('Point')) for x in row[:-1] ]
AttributeError: 'int' object has no attribute 'split'


texts: Legend not showing properly in LayerList

texts:
In my web app using ArcGIS API for Javascript 4.18, I am bringing in 3 layers served from a mapServer REST Service using MapImageLayer. I also have a LayerList widget to display these layers. My problem is that while the layer display correctly, the legend doesn't show up in the collapsed sublayer. I does show up in the overall collapsed menu for the entire layer, but that's not where I want it. I want each sublayer to have its own legend next to it:

My MapImage Layer is configured as below:
var Layers_right = new MapImageLayer({
    url: mapServer,
    Title: "Layers",
    sublayers: [
        {
            title: "Streamline",
            id: 0,
            visible: true
        }, {
            title: "Effective Flood Hazard Area",
            id: 1,
            visible: true
        }, {
            title: "Depth 10%",
            id: 9,
            visible: false
        }
    ]
});

My layerList is configured as below:
 var layerList_right = new LayerList({
        view: popupView2,
        listItemCreatedFunction: function (event) {
            const item = event.item;
            if (item.layer.type !== "group") {
                // don't show legend twice
                item.panel = {
                    content: "legend",
                    open: false
                };
            }
        }

What do I have configured incorrectly here? How can I get the legend to show up when I click the collapse menu next to the sublayer?


texts: Select every other n rows in ArcGIS for Desktop?

texts:
I am trying to select n rows then skip n rows and then select the next n rows until the completion of all rows in a feature class in a file-geodatabase.  Below is an example where n =3:

Select
Select
Select
Not selected
Not selected
Not selected
Select
Select
Select

I basically need to make a selection like this and populate the rows that are selected.


texts: What can help gain contribution points for GISP certificate?

texts:
Didn't see this discussed before.

What earns contribution points towards GISP certification?

Does volunteering to answer on forums help contribute?


texts: Recreating this cartographic style in ArcMap

texts:
I am trying to recreate this symbology as shown in the image below. I came across this image on ESRI's service area analysis help page.

From my previous experience and intuition, I can say that roads are symbolized as white and road are on top of the service area polygons. But I am unable to understand

how to give that white glow effect in ArcMmap. Especially the polygon in red, it has this very nice white border with another polygon in beige and we move out towards the boundary of this red polygon its color transition from white to darkish red. I tried to achieve that with the gradient color symbology in ArcMap but failed miserably.
how to give this shadowy effect that has been given to the second polygon shown in beige/yellowish sandy color

I am using ArcMap but from my photoshop experience, I can say that it has something to do with strokes. I am wondering how would that be achieved in ArcMap


texts: Differing checksum on two otherwise identical geopackages

texts:
I found that two exactly identical Geopackages do not share the same checksum as one would probably expect (and is vital when syncing files using rsync for example). Here is a reproducible example:
Get some sample data:
$> wget https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/gadm36_ZWE_gpkg.zip
$> unzip gadm36_ZWE_gpkg.zip

Now use ogr2ogr to copy the data to two different files with exactly the same parameters:
$> ogr2ogr -f GPKG test1.gpkg -nln testlayername gadm36_ZWE.gpkg
$> ogr2ogr -f GPKG test2.gpkg -nln testlayername gadm36_ZWE.gpkg

If we check the checksums they happen to differ:
$> md5sum test*.gpkg
ad2fe21516a94ab9f2c0bca37c171028  test1.gpkg
240b65bac70687d8fc65c70c9c4dbb46  test2.gpkg

So after some experiments, I found that it must have to do something with the "Geo" component, as
ogr2ogr -f sqlite test[1|2].sqlite -nln testlayername gadm36_ZWE.gpkg

would generate the same chechsum for both outputs while
ogr2ogr -f sqlite test[1|2].sqlite -nln testlayername gadm36_ZWE.gpkg -dsco "SPATIALITE=YES"

would not generate the same chechsum for both outputs.
Am I completely off here? Does anybody know what causes the difference in checksums and is there any way to prevent that when using GPKGs?


texts: How estimate or get the value of PDOP from HDOP value?

texts:
I have a gps device that are reporting me the HDOP value and the number of satellites around, but I have as requirement to get the PDOP also, the GPS device doesn't give that value.
Can I get PDOP value from HDOP and the number of sats?


texts: Saving focal statistics in memory using ArcPy?

texts:
I'm currently working on a project where I have written a script that is suppose to do several things. one of these is creating a "quick" focal statistics, from which I extract the values to some points. Afterwards I wont need the focal statistics output, so I was wondering if it was possible to save the output "in memory". Right now, each time I run:

        focalstats = arcpy.sa.FocalStatistics(inDSM, NbrCircle(7, "MAP"), "MINIMUM")

the line creates an output file called "FocalSt_dsm_x" (x is a number) in my environment workspace.

Of course, I could delete the file each time, but it seems a bit troublesome.


texts:
I am using Pygrass in a a Python Grass GIS stand alone application.
I am using GRASS GIS 7.0.3 on a Windows machine and Python 2.7.4.
I am initializing a new Grass environment. I can work in this environment, and apply grass functions (e.g. load shapefiles,...)

I would like to access the Vector data directly using this example:

from pygrass.vector import VectTopo
municip = VectTopo('boundary_municp_sqlite')
municip.is_open()
False
municip.mapset
''
municip.exist()  # check if exist, and if True set mapset
True
municip.mapset
'user1'

Unfortunately I get the following Error during the command municip.exist():

ERROR: MAPSET is not set

But I can work in the environment, the variables Location and Mapset("PERMANENT") are set. I also used g.region to set the region.
I tried to use the module from grass.pygrass.gis import Mapset.
This also gives me an error when I try to access the function:

from grass.pygrass.gis import Mapset
m = Mapset()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-59-73310da5f35b> in <module>()
----> 1 m = Mapset()

C:\Program Files (x86)\GRASS GIS 7.0.3\etc\python\grass\pygrass\gis\__init__.py
in __init__(self, mapset, location, gisdbase)
    265         self.gisdbase = gisdbase
    266         self.location = location
--> 267         self.name = mapset
    268         self.visible = VisibleMapset(self.name, self.location, self.gisd
base)
    269

C:\Program Files (x86)\GRASS GIS 7.0.3\etc\python\grass\pygrass\gis\__init__.py
in _set_name(self, name)
    290
    291     def _set_name(self, name):
--> 292         self._name = _check(name, join(self._gisdb, self._loc), "MAPSET"
)
    293
    294     name = property(fget=_get_name, fset=_set_name,

C:\PROGRA~2\QGISPI~1\apps\Python27\lib\ntpath.pyc in join(a, *p)
     71             b_wins = 1
     72
---> 73         elif isabs(b):
     74             # This probably wipes out path so far.  However, it's more
     75             # complicated if path begins with a drive letter:

C:\PROGRA~2\QGISPI~1\apps\Python27\lib\ntpath.pyc in isabs(s)
     55 def isabs(s):
     56     """Test whether a path is absolute"""
---> 57     s = splitdrive(s)[1]
     58     return s != '' and s[:1] in '/\\'
     59

C:\PROGRA~2\QGISPI~1\apps\Python27\lib\ntpath.pyc in splitdrive(p)
    123     """Split a pathname into drive and path specifiers. Returns a 2-tupl
e
    124 "(drive,path)";  either part may be empty"""
--> 125     if p[1:2] == ':':
    126         return p[0:2], p[2:]
    127     return '', p

TypeError: 'NoneType' object has no attribute '__getitem__'

I am not sure if both errors are connected somehow.


texts: ArcPy Error - 'continue' not properly in loop

texts: Is there a way to get the continue function to work? Here is my code:


texts: Maximo Spatial: Convert UTM to lat/long with Python 2.7.0 (using only the standard Python library)

texts:
Maximo Spatial is a work order management system in which I can write scripts in Python 2.7.0.

I want to write a Python 2.7.0 script in Maximo that will convert UTM coordinates to latitude/longitude decimal degree coordinates. Ideally, the transformation error would be less than 1 meter.

The catch:

It is only possible for me to use standard python libraries in my Python scripts in Maximo Spatial. I can't import external libraries.

Is there a way to convert UTM coordinates to lat/long using out-of-the-box Python 2.7.0?


texts: How to add a legend to rasterlayer in Qgis2Web plugin while exporting maps?

texts:
Is it possible to add a legend for Rasterlayer directly from the Qgis2Web Plugin or do I need to add the code manually?
I added an image of an example (DEM raster) and a Vector Layer.

It doesn't matter if I choose leaflet or OL3, there is no option for any legend in this plugin as it seems.
Also there seems to be no option for showing the raster value as a popup like there is for vector layers.

When exported the in the layer list it appears a legend for the vector layer, but for the raster it's just the name.

If there is really no option for adding a legend via this plugin, could anybody help with adding a code manually?
Which of the resulting file I need to change?


texts: Simple thematic mapping of shapefile using Python?

texts:
I wish to visualize geographical data in Python, without using ArcGIS/ArcPy, and make a map.

On the internet I found how to make a thematic map using Python:

Here is some code that I tried:

import shpUtils
import matplotlib.pyplot as plt

shpRecords = shpUtils.loadShapefile('C:\\Users\\shapefile.shp')

for i in range(0,len(shpRecords)):
    x = []
    y = []
    for j in range(0,len(shpRecords[i]['shp_data']['parts'][0]['points'])):

    tempx = float(shpRecords[i]['shp_data']['parts'][0]['points'][j]['x'])
    tempy = float(shpRecords[i]['shp_data']['parts'][0]['points'][j]['y'])
    x.append(tempx)
    y.append(tempy)
    plt.fill(x,y)

plt.axis('equal')
plt.title("Testing")
plt.show()

However, when I run this, it gives me random colors.

If I want to visualise one certain column of my shapefile, how can I implement this using similar code?

This is very unclear in the link provided above where he only discusses the usage of colors...

Would I perhaps need extra modules to accomplish this, like descartes and PySAL?


texts: Storing ArcObject inside BLOB or XML field of geodatabase?

texts:
Is it possible to store an ESRI object inside a blob field inside a geodatabase? If not, perhaps a serialized XML representation?

I'm thinking about storing IFeatureRenderers inside these fields.

Any tips? Anyone has ever done this?


texts: Excel Sheet to Table - Allow User to Select

texts:
I have a model flow that allows a user to input an excel worksheet which then converts a sheet within the wksht to a table and joins that table to a feature class. Currently, the Excel to Table tool I am using will only convert the sheet within the worksheet that is either the first sheet in the list in the worksheet or will convert the particular sheet that is specified by a name in the tool. My issue is, there is no standard naming convention for the particular sheet that will need to be converted and one excel workbook can have multiple sheets within. I only want one particular sheet within the wkbk to be converted to a table and would Ideally like the user to be able to select that sheet. Is there any way to allow user input on this? Here is an image of my current model below:

Where the Sheet(optional) input is, I would like the user to be able to select the appropriate sheet but I cannot make this a model parameter.


texts: How to control the zoom amount in Cesium.Camera.flyTo ?

texts:
The following code flies the camera to a particular destination. But how do I control the zoom on the particular destination?

   viewer.camera.flyTo({
      destination : Cesium.Cartesian3.fromDegrees(-117.16, 32.71, 15000.0),
      duration : 20.0
    });

Changing Cesium.Camera.defaultZoomAmount doesn't seem to work for me.


texts: Python Crashes when Running ArcPy's ExportToPDF(), ExportToPNG(), ExportToJPEG()

texts:
I am attempting to loop through a folder of ~400 .mxd's and save each to a .pdf in another folder.  All tested arcpy.mapping export functions (pdf, jpeg & png) kill python without a chance to debug or any indication why.  Running script in PythonWin gives error message "PyWin32 has stopped working."  When run from python prompt within ArcMap, program just closes.

Script behaves inconsistently, sometimes makes 10 maps before crashing, other times only 1 or 2.  All mxds were created from the same arcpy.mapping script and all are valid (can open in ArcMap and export without issues).  There is an ArcGIS Online basemap layer in each map (Oceans), I tried unchecking this layer in the first few maps and it did not have an effect on the script, python still crashed.

Tested script on 2 different laptops with same result.  Both running ArcGIS Desktop 10.2.2, Python 2.7...

Windows 7 Pro 64-Bit SP1, Intel Core i5, 8GB Dual-Channel RAM, Intel HD Graphics 4000

Windows 7 Pro 64-Bit SP1, Intel Core i7, 16GB Dual-Channel RAM, Intel HD Graphics 4000, 2048MB NVIDIA Quatro K2000M

Saw a few similar posts here and on other message boards, but no working solutions.

I had no problem generating 400 maps using arcpy.mapping and cannot believe that I cannot save them out to a file.  Is this a bug?

import arcpy, os

ws = r"C:\temp"
map_folder = os.path.join(ws, "maps")
pdf_folder = os.path.join(ws, "pdfs")

arcpy.env.workspace = map_folder
arcpy.env.overwriteOutput = True

# generate list of map documents in folder to loop through
map_list = arcpy.ListFiles("*.mxd")

###
def exportAISMap(mxd_path, out_path):
    mxd = arcpy.mapping.MapDocument(mxd_path)
    # arcpy.mapping.ExportToJPEG(mxd, out_path)
    # arcpy.mapping.ExportToPNG(mxd, out_path)
    arcpy.mapping.ExportToPDF(mxd, out_path)
    print "Exported map file: " + str(out_path) + "\n"
    # mxd.save()
    del mxd
###

print "Saving out " + str(len(map_list)) + " map documents\n"

for map_file in map_list:
    print map_file

    mxd_path = os.path.join(map_folder, map_file)

    pdf_file = map_file.replace(".mxd", ".pdf")
    pdf_path = os.path.join(pdf_folder, pdf_file)

    exportAISMap(mxd_path, pdf_path)


texts: Popup with geojson data loaded using `Leaflet.FileLayer`

texts:
I'm trying to create an interactive map using the Leaflet.FileLayer plugin. This plugin allows me to load geojson files with description but I don't know how to show this description in Popup when data is loaded from separate file. Normally, I call data using this code:

var miejsca_pamieci_l = L.geoJson(miejsca_pamieci, {
    pointToLayer: function (feature, latLng) {
      return new L.Marker(latLng, {
        icon: new myIcon({
          iconUrl: 'icon/cemetery-24.png',
          iconSize: [24, 24],
          iconAnchor: [12, 24],
          popupAnchor: [0, -24]
        })
      })
  },
onEachFeature: function( feature, layer ){
    layer.bindPopup (feature.properties.Opis, {  maxWidth: "auto"})
    }  });

But how to change the code to show popup using FileLayer? This is my script in html code:

  (function (window) {
    'use strict';
    function initMap() {
        var L = window.L;
        var osm = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
            attribution: 'Map data © 2013 OpenStreetMap contributors'
        });
        var map = L.map('map', {
            center: [0, 0],
            zoom: 2
        }).addLayer(osm);
        var style = {
            color: 'red',
            opacity: 1.0,
            fillOpacity: 1.0,
            weight: 2,
            clickable: false
        };
        L.Control.FileLayerLoad.LABEL = '';
        L.Control.fileLayerLoad({
            fitBounds: true,
            layerOptions: {
                style: style,
                pointToLayer: function (data, latlng) {
                    return L.circleMarker(
                    latlng,
                    { style: style }
                    );
                }
            }
        }).addTo(map);
    }
    window.addEventListener('load', function () {
        initMap();
    });
}(window));


texts:
I am trying various ways to have a programatic way of combining various bands of a tiff file into the red, green, blue channels.

I tried the CQL way, eg see below

This is my SLD

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
 xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
 xmlns="http://www.opengis.net/sld"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>custom_raster</Name>
    <UserStyle>
      <FeatureTypeStyle>
        <Rule>
          <Name>rule1</Name>
          <Title>Opaque Raster</Title>
          <Abstract>A raster with 100% opacity</Abstract>
          <RasterSymbolizer>
            <Opacity>1.0</Opacity>
            <ChannelSelection>
              <RedChannel>
                <SourceChannelName>${env('red',0)}</SourceChannelName>
              </RedChannel>
              <GreenChannel>
                <SourceChannelName>${env('green',1)}</SourceChannelName>
              </GreenChannel>
              <BlueChannel>
                <SourceChannelName>${env('blue',2)}</SourceChannelName>
              </BlueChannel>
            </ChannelSelection>
          </RasterSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

I call it using ENV=red:1;blue:3;green:5 but I am getting the following message

Caused by: java.lang.IllegalArgumentException: Band number -1 is not valid.
...
Caused by: java.lang.NumberFormatException: For input string: "${env('red',0)}"
...

It doesn't seem to allow CQL expressions where I am placing them... yet doesn't complain when I validate. Any idea?


texts:
I'm having trouble setting up TauDEM to work in the processing tools of QGIS 2.0.1.

I download TauDEM from here and select "TauDEM 5.2 Command Line Executables for 32 bit Windows PC". I extract them to a directory called "C:\work\taudem".

I then go into QGIS processing menu > "Options and configuration" under TauDEM I set the "TauDEM command line tools folder" to "C:\work\taudem".

When I execute a TauDEM tool I keep getting errors saying:

  Oooops! The following output layers could not be open.


texts: QGIS feature form : Force attribute to save as uppercase

texts:
For example I have feature form with colum name NAME (with names of contries).

I want that all counties names are in capital letters.

As users enter the names in lowercase, uppercase, mixed ... I am setting all nemes to uppercase on this way:

attributre table -> Field calculator -> Expression ->  upper("NAME")

and this is working fine, but the idea is that I don't have to repeat it every time, is there any way to save ('field properties') in the project so that no matter what the user enters, it will be capitalized and saved into table.