texts: Populate an attribute table field using the 'startswith function

texts:
I am hoping that someone may be able to help me.  I am new to ArcGIS, and even newer to Python.

Working in ArcGIS 10.1, I have done is created a new text field in an attribute table, which I have called "Ag".  I want to populate this field with the text from another field (called Address) which is in the same attribute table BUT only if that text starts with the word "Field".

I have written the following function but it is not working, and I am not sure where I have gone wrong, or if I am on the wrong path completely.

def Reclass(agri):
    if (!Address!.startswith("Field")):
      return !Address!
    else:
      return "NA"

I have searched the forums and the ESRI documentation for a solution, but cannot seem to find what I am looking for.  I hope that I have given all the relevant information for someone to be able to help me with my conundrum, but please let me know if there is more (or better) information I should supply.

Many thanks in advance!


texts: Changing attribute values based on consecutive features

texts:
I'm trying to determine the limiting capacity of a pipe based on the capacities of the pipes downstream of it. Basically, if Pipe A is downstream (ahead of) Pipe B and Pipe B has a capacity that is lower than Pipe A's, then Pipe A's limiting capacity would be = Pipe B's capacity. Each pipe has an upstream and downstream manhole number (INLETNODE and OUTLETNODE in the code I've attached) that I'm attempting to use as a way to move through the pipes consecutively. To start off, each pipes Pipe Capacity field and Limiting Capacity field are the same.

Here's the code I've come up with so far, as of right now the limiting capacity field doesn't get edited at all:

ds_pipe_dict = {}
with arcpy.da.SearchCursor(s_pipes, ["INLETNODE", "OUTLETNODE", "PIPE_CAPACITY", "LIMITING_CAPACITY"]) as cursor:
    for row in cursor:
        if row[0] not in ds_pipe_dict:
            ds_pipe_dict[row[0]] = [row[1], row[2],row[3]]

fields = ["INLETNODE", "OUTLETNODE", "PIPE_CAPACITY", "LIMITING_CAPACITY"]
with arcpy.da.UpdateCursor(s_pipes, fields) as cursor:
    for row in cursor:
        pipe_cap = row[2]
        lim_cap = row[3]
        if row[2] == row[3]:
            ds_mh = row[1]
            pipe_cap = row[2]
            lim_cap = row[3]
            next_lim = lim_cap

            if ds_mh in ds_pipe_dict:
                if lim_cap > next_lim:
                    lim_cap = ds_pipe_dict[ds_mh][2]
                    ds_mh = ds_pipe_dict[ds_mh][0]
                else:
                    ds_mh = ds_pipe_dict[ds_mh][0]

        row[3] = lim_cap
        cursor.updateRow(row)


texts:
I want to use a non CartoDB-installed typeface (one from typography.com) for the place labels (closer than 14z) on the map at http://outgoingnyc.com. So far if I try to do it in the css, it fails (blank map).

Here's how I brought in the data:

 cartodb.createLayer(map, layerSource)
    .addTo(map)
    .done(function(layer) {
      layer.setInteraction(true);
      layer.on('error', function(err) {
          console.log('error: ' + err);
        });

    var subLayerOptions = {
      sql: "SELECT * FROM " + outgoingTable,
      cartocss: $("#markerStyle").text(),
      interactivity: 'year_opened, cartodb_id, name, year_opened, year_closed, audience, audience_2, audience, genre, genre2, source, pub_desc',
      infowindow: true
      } //end subLayerOptions

Here's the cartoCSS link (abbreviated):

<style type="cartocss/text" id="markerStyle">
  #outgoing_masterv2{
    marker-fill-opacity: 1;
    marker-line-color: #FFF;
    marker-line-width: 0;
    marker-line-opacity: 1;
    marker-placement: point;
    marker-type: ellipse;
    marker-allow-overlap: true;
    }
...
</style>

I'm guessing that the way I've done this makes me only able to access the CDB pre-installed type. Is there some way I don't know about to change this though?

The font-family is 'Ideal Sans A', 'Ideal Sans B'

font-family: 'Ideal Sans A', 'Ideal Sans B';


texts: Why do colour assignments in the legend not correspond to the display?

texts:
In QGIS 2.18.12 (qt4) on Linux, I repeatedly observe strange legend behaviour when new layers are loaded or generated using the Raster Calculator.

The legend will show colours for two classes (values). In the display window, the colours are apparently reversed.

See the following image:

At 1, we see that in the legend, nan is black and 0 is white.

At 2, I am clicking on the black field in the mgrid layer.

At 3, the value is 0. If I click outside the black field, the value is nan, which is the opposite of the expected behaviour.

If I reverse the gradient in Properties | Style | Band rendering, then the legend is "corrected" (nan is transparent, 0 is black) but in the display window, there are no discernable fields at all:

Updating the min/max statistics (with load min/max values) in the band rendering settings fills both fields with 0 and results in a legend that looks like this:

The image does not change (it is as it is in the first image above).

What causes this, and how do I fix it?


texts: Find nearest point for each point in same layer using ArcGIS network analyst?

texts:
My question is similar to this How to find nearest point for each point in current layer using ArcGIS Desktop?
But I am looking to use arcgis network analyst. I tried using the origin destination matrix and used the same layer for both origins and destinations but I got the nearest point being itself.


texts: There is web based bicycle routing planner that considers accidents but I cannot find it anymore, can someone help me?

texts:
I remember it uses Google as underlying map and has also a function to rank preferences of the path using forms in which the user can rate routes characteristics, I also remember that is for Los Angeles (but I am not sure about this). I came across this web map a couple of months ago but I did not save it as favourite. Now I have been searching continually for two days in many ways but I can't find it anymore.

Does someone know what I am talking about and please provide me the page URL?


texts: Including elevation information while temperature interpolation

texts: I am working on ArcMap 10.5 for temperature interpolation. I am using IDW for interpolation. I need to compare the results with and without elevation data. However, I am unable to understand how to add elevation and temperature interpolation layer? Basically, I want an interpolated temperature raster that reflects elevation effect as well.


texts:
how can convert dwg to sqlitedb or mbtiles to use it as map in galillo for mobile or any other extension that can work in gps program for mobile
dwg contain map i n stereografic reference and how can georeference to wgs84 and view it in correct place


texts: TIGER/LINE 2013 Feature, Edge and, Address Range Data Sets

texts:
I am looking at the Tiger 2013 data sets, in particular the Edge, Feature and, Range tables.

It seems that the Edge tables have all the data elements that are in the Feature and Range sets. Is this true for the entire US sets? If the Edge tables do in fact contain all the data elements that are in the Feature and Range tables, why have 3 separate tables?


texts: Remove a Row in an attribute table using UpdateCursor?

texts:
I'm trying to create a piece of code to remove rows where the field "State_Name" equals Hawaii or Alaska from my shapefile's attribute table. I've tried two different codes and they both seem to run without error in the ArcGIS Pro Python window, however the attribute table doesn't change.
First attempt:
with arcpy.da.UpdateCursor(outputft, ["State_Name"]) as cursor:
    for row in cursor:
        where="State_Name is HAWAII or ALASKA"
        cursor.deleteRow()

Second attempt:
with arcpy.da.UpdateCursor(outputft, ["State_Name"]) as cursor:
    for row in cursor:
        if row[0] == 'HAWAII':
            cursor.deleteRow()

I've closed, reopened, refreshed the attribute table as well but it doesn't seem to work. I'm new to Python.


texts: Calculate distance between a coordinate and a county in GeoPandas

texts:
I have a single coordinate and a GeoPandas Dataframe representing US counties whose geometry is polygons that form the boundary of counties. Each of these polygons is also given through lat/lon pairs.

I need to find the distance, in miles, between my single coordinate (in lat/lon) and a county. When I use the distance function like this:
point.geometry.distance(county.geometry), the distance is calculated through use of Euclidean distance. Since the points are given in terms of lat/lon, the Euclidean distance does not give the real distance between points.

How do I find the distance between a coordinate and a county using GeoPandas?


texts: How can I color code national vs subnational boundaries in QGIS?

texts:
I'm trying to make a regional map of West Africa at the subnational level. I can't seem to figure out how to distinguish the subnational and national boundaries - right now, they're all black, and I'd like the subnational boundaries to be in a different color than the national boundaries.

Does anyone know how I can do this in QGIS?


texts: creating 2 dragable layers at the same time in OpenLayers?

texts:
I have been been working with this for a while now. Whenever you add 2 seperate draggable layers, openlayers will only activate the last one added. I wish it could be as simple as the code below. But that doesn't work

var dragFeatureend = new OpenLayers.Control.DragFeature([vectorLayer1, vectorLayer2], {autoActivate: true, onDrag: onCompleteMove});

I found a solution by making the layers selectable and calling events when each is clicked. But is there any work around other than making them selectable? Thanks!


texts: Using intersection tool with geodatabase and vector layer in QGIS

texts:
I want to use the intersection tool in QGIS version 2.18 for two layers:

Corine land use data downloaded from: http://land.copernicus.eu/pan-european/corine-land-cover/clc-2012 (a geodatabase, vector format)
buffer layer in vector format.

When I try to use the intersection tool, I get an error:

  "Input layer A contains invalid geometries (feature 1255). Unable to
  complete intersection algorithm."

Layer A is the Corine layer.

How can I get rid of this error and use the intersection tool properly?


texts: Complex polygon not displaying properly in QGIS?

texts:
I append two images of the same shapefile which contains a single polygon with several "lakes". In SAGA it displays correctly with all the lakes as transparent but in QGIS the larger lakes contain the colour fill.

Is there something I'm not doing correctly in QGIS 2.6.1, or is there an option for changing the way it imports a polygon?

QGIS

SAGA


texts: How can I use line vectors to alter my DEM?

texts:
I'm trying to use vector lines to change the height of the areas of my DEM that they overlay.  They trace rivers, and i'm hoping to assign arbitrarily high values to the rivers so that my least-cost path analyses using r.walk don't constantly try to follow the riverbeds.

What's the best way to do this?


texts: Landsat level-2 NDVI valid range with saturation

texts:
I have Landsat level-2 NDVI tiles I am using in ArcMap 10.4.1.; however, the level-2 NDVI data in ArcMap ranges from -5,189 to +20,000.

The spectral image product guide (which I've linked) says that the valid range for level-2 NDVI is -10,000 to 10,000 with a scale factor of 0.0001, which makes sense. The saturate value is 20,000. I don't understand what the saturate value is and why the upper bound of the level-2 NDVI layer is 20,000.

I understand that to get NDVI values I need to multiply the level-2 layer by 0.0001. However, I'm unsure what to do about the saturate values since they are above the valid range. I'm guessing that the saturate value indicates pixels that are unusable; in other words, once I multiply the level-2 layer by 0.0001, I would then throw out all pixels >1. Is this correct?

https://landsat.usgs.gov/sites/default/files/documents/si_product_guide.pdf


texts: Indonesia Level 3 Administrative Shapefile

texts:
I've been looking for Indonesia shapefile for level 3 or more administrative boundary and I still haven't found it yet.
I found in (gadm.org) but it only level 2.
So anyone know where I could find it?
I would prefer a single Indonesia shapefile instead of downloading each province shapefile one by one.
Also it would be great if it has properties like the name of the area, etc.


texts: Unable to establish Oracle enterprise connection using ArcMap

texts:
When I run the following command

sqlplus user/mysecretpass@(DESCRIPTION=(ADDRESS=(PROTOCOL=TCP)(Host=10.1.250.81)(Port=1521))(CONNECT_DATA=(SID=gisdb)))

I am able to establish connection and result is as follows

I am also able to connect to same database using SQL Developer but when I try to create database connection from ArcMap 10.6.1 I get the following error

I have also followed the instructions here, but I am unable to solve this issue. Can anyone let me know how to solve this problem?


texts: What is the projected coordinate system equivalent for ETRS89?

texts:
ETRS89 (EPSG: 4258) is the EU-recommended frame of reference for geodata for Europe 1. This coordinate system is geographic. However, in the software world, we need to transform it into a projected coordinate system in order to take measurements.
Is the ETRS89-extended / LAEA Europe (EPSG: 3035) the right equivalent?


texts: Im builder a Buffer based on dealer sales.  The inputs are Dealers points and sales points. I need the buffer to increase distance until 40% of sales volume is reached. Im not sure sure is I should attempt this  For and while or should I be using get value and create a new field instead?


texts: Unable to obtain CRS/Projection info in rasterio for a projected raster

texts:
I have a 9-band sentinel-2 surface reflectance image that I combined in ArcGIS. When I load it into a python script, dst.profile indicates that "CRS: None". However, when i run a gdalinfo command it indicates it has a projected coordinate system as below:

Coordinate System is:PROJCRS["WGS 84 / UTM zone 17N",
BASEGEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]]],
CONVERSION["UTM zone 17N",
    METHOD["Transverse Mercator",
        ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",-81,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["Scale factor at natural origin",0.9996,
        SCALEUNIT["unity",1],
        ID["EPSG",8805]],
    PARAMETER["False easting",500000,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]],
CS[Cartesian,2],
    AXIS["(E)",east,
        ORDER[1],
        LENGTHUNIT["metre",1]],
    AXIS["(N)",north,
        ORDER[2],
        LENGTHUNIT["metre",1]],
USAGE[
    SCOPE["unknown"],
    AREA["World - N hemisphere - 84┬░W to 78┬░W - by country"],
    BBOX[0,-84,84,-78]],
ID["EPSG",32617]]

When I run this python snippet

sent_data = rasterio.open(image)
with rasterio.Env():
    sent_data.profile

I get the following result:

{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': None, 'width': 5490, 'height': 5490, 'count': 9, 'crs': None, 'transform': Affine(20.0, 0.0, 399960.0,
   0.0, -20.0, 3900000.0), 'blockxsize': 128, 'blockysize': 128, 'tiled': True, 'compress': 'lzw', 'interleave': 'pixel'}

I have also tried

sent_data.gcps

And also tried using gdal python bindings as given in a few similar questions here, but still not working. I need to be able to obtain the CRS in order to apply it to another output raster so that it can be projected properly.

Any suggestions would be greatly appreciated.


texts: Import non spatial data into SQL Server using ogr2ogr

texts:
I am trying to import a CSV files in to SQL Server 2012 using ogr2ogr. The spatial data imports without any problems, but I cannot work out why the non spatial tables fail. I can import the data using other tools but it would be convenient if I could import everything using ogr2ogr so I can script it more easily.

My VRT file looks like this...

    <OGRVRTDataSource>
      <OGRVRTLayer name="test">
        <SrcDataSource relativeToVRT="1">test.csv</SrcDataSource>
        <FID>id</FID>
        <Field name="id" src="Id" type="Integer"  />
        <Field name="Test" src="Attrib" type="string" width="100" />
     </OGRVRTLayer>
    </OGRVRTDataSource>

And my command line:-

   ogr2ogr -append -f MSSQLSpatial "MSSQL:server=SQL01;database=TESTDB;trusted_connection=yes" test.vrt

The first error I get is:-

ERROR 1: Column ogr_geometry requested for geometry, but it does not exist.

After that there are errors trying to create the attribute field because the table does not exist. If I manually create the table in SQL Server with a single field named "id", then this command works (still with error messages).


texts: Can two people work on the same QGIS project in real time?

texts: I was wondering if QGIS allows for multiple users to edit the same project at the same time from different computers, in real time (similar to how two people can work on the same google doc at the same time).


texts: Using only distinct values to iterate as atlas features in QGIS?

texts:
I want to iterate several maps using the atlas feature of the print composer in QGIS 3.4. The field I want the maps to iterate over has several duplicate values (ie several features belong to the same group and are labelled as such) but at the moment, atlas produces one map for each feature regardless of the fact that some have the same value in the field they're being iterated over.

What I want it to do is produce a map which is zoomed to fit all the features labelled 'A' and another map zoomed to fit all the features labelled 'B', not one map for each feature.

I could merge the features by their labels, but then I'd lose the remaining attributes which do vary and I'd like to present them in a linked attribute table, so they need to be preserved within the layer really.

All I really want QGIS to do is assign the same @atlas_featureid to all the features labelled 'A' and a different @atlas_featureid  to all the feature labelled 'B', my existing set-up will do the rest.


texts: Can anyone point me to a public WPS (Web Processing Service)?  I need to do some testing on it and results in Google and Bing show nothing.


texts: Action after iteration in ArcGIS ModelBuilder?

texts:
I've build a model using the ArcGIS ModelBuilder. The model uses an iterator that clips all the features in a specific geodatabase. The iterator-part is nested in a submodel. I want that geodatabase ("Clipping-InputGDB") to be deleted after all the features within it have been clipped. Is there a way to do that?

Screenshot shows submodel.


texts: GDAL polygonize only certain raster values

texts: I want to polygonize a raster file with GDAL. I only want to rasterize cells with a value of 0. Is it possible to specify that within the polygonize function or do I have to modify my raster before hand?


texts: How to set QGIS temporal controller time range using a timestamp from the raster's filename or TIFFTAG_DATETIME

texts:
I am trying to work with and eventually animate between 20-100 raster layers using the QGIS temporal controller.  It seems silly to manually set the date on each layer in the "temporal" tab of the layer properties when this information is contained both in the layer filename (in the form 20210304_084108 for 8:41 AM today, for example), as well as in the TIFF tag (TIFFTAG_DATETIME=2020:12:21 08:41:08).
I cannot see mention of a way to do this for rasters.  Unlike vector layers, there is no "configuration" option in the temporal properties tab that would allow me to set the start date based on an expression (that could, conceivably, access the filename or TIFF tag).
I do not program myself, but if I did, I imagine there would be a way to do this if I dug into PyQGIS and used the GetMetadataItem("TIFFTAG_DATETIME") of GDAL.
Has anyone solved this?  Since much raster data (from satellites, for example) has the date in the filename and TIFF tags, I imagine this might be a common need.


texts: How can I shift raster layer for making cloud-shadow mask?

texts: I want to make cloud shadow mask in QGIS. I have raster cloud mask and I have to shift it by X,Y coordinates. Is there any easy solution?


texts:
I am exporting OSM data (in .osm) files to use in QGIS. I am interested in having the street names displayed. I know they are in the OSM data set because they are displayed when I browse the data online, but when I view my downloaded files and open the attribute tables for streets, most of them simply have "null" or "label" as values.

What am I doing wrong?

How can I extract data from OSM which includes the street names?


texts: Finding angle between vertices to estimate curvature in the same feature class using python script?

texts:
Background info

What I'm attempting to create is an estimation of curvature, by finding the angle between vertices. The vertices layer was created from using the "Split line at vertices tool" on a roads layer. From there I used calculate geometry to get xy data for the line.

The problem

Create a function which uses a line's start and endpoint xy data to determine the angle between it's neighboring line segment(s). What I'm struggling with is determining what how the lines connect. The angle can be determined using math.atan2(), but I have no idea how find the connecting segments.

x1- x coordinate for the start of the line segment

y1- y coordinate for the start of the line segment

x2 -x coordinate for the end of the line segment

y2- y coordinate for the end of the line segment

Basic Pseudo code

function angle(x1,y1,x2,y2)

    Given x1,y1 in column[observed line segment]

    Find x2,y2 in column[neighboring line segment] == x1, y1 in column[observed line segment]

    If lines segments have matching coordinates find angle between lines.

   return angle


texts: How to highlight track in KML on hover with OpenLayers

texts:
I have KML layer with multiple tracks and I would like to highlight them individually on mouse hover. I tried to implement this like in the example found in OpenLayers examples, but couldn't get it work.

In the example there is following code that adds polygon as feature. In my case, how do I add KML-tracks as feature so that I can highlight them with SelectFeature?

var feature = new OpenLayers.Feature.Vector(
    OpenLayers.Geometry.fromWKT(
        "POLYGON((28.828125 0.3515625, 132.1875 -13.0078125, -1.40625 -59.4140625, 28.828125 0.3515625))"
    )
);
vectors.addFeatures([feature]);

You can find simplified version of my index.html and KML file here: https://gist.github.com/1379118


texts: Intersecting two shapefiles from Python or command line?

texts:
I'm new to working with geo data and trying to address a similar question to the one posted as How to calculate the size of a particular area below a buffer in QGIS

Using command line and/or shapely/fiona, in the most efficient and elegant way.

I created two shapefiles with their respective numerical attributes:

small areas with data on population
police districts with crime stats(single number)

They are independent.

I'd like to recalculate those such that both estimates are present in each of the shapefiles, eg. estimate of crime within small area (take data from 2 and look for all areas in 1 that are contained in it, cut those that intersect and do maths).

I have the maths figured out and stuck looking for elegant ways and efficient data formats to do intersections and cut out the overlaps to appropriate boundaries.

Update following the response of @gene

The adapted code does not finding any intersections in the data, which is structured as follows (and surprising since it's from official sources):

polSHP_upd:  police precincts across SA
sal_subSHP: subset of small regions with population data in Cape Town

The files live here:

https://github.com/AnnaMag/GIS-analyses

What am I missing?

polSHP_upd = 'crime/GIS-analyses/polPres_updated4crimes/polPres_updated4crimes.shp'
sal_subSHP = 'crime/GIS-analyses/sal_sub/sal_sub.shp'

import geopandas as gpd
g1 = gpd.GeoDataFrame.from_file(polSHP_upd)
g2 = gpd.GeoDataFrame.from_file(sal_subSHP)
data = []
# A: region of the police data with criminal record
# C: small area with population data
# we look for all small areas intersecting a given C_i, calculate the  fraction of inclusion, scale the
# population accordingly: area(A_j, where A_j crosses C_i)/area(A_j)*  popul(A_j)
for index, crim in g1.iterrows():
   for index2, popu in g2.iterrows():
      if popu['geometry'].crosses(crim['geometry']): # objects overlaps to partial extent, not contained
         area_int = popu['geometry'].intersection(crim['geometry']).area
         area_crim = crim['geometry'].area
         area_popu = popu['geometry'].area # there is a Shape_Area field in properties: to check precision
         popu_count = popu['properties']['PPL_CNT']
         popu_frac = (area_int / area_popu) * popu_count# fraction of the pop area contained inside the crim
         data.append({'id': (index1, index2) ,'area_crim': area_crim,'area_pop': area_popu, 'area_inter': area_int, 'popu_frac': popu_frac} )
         # data.append({'geometry': crim['geometry'].intersection(popu['geometry']), 'crime_stat':crim['crime_stat'], 'Population': popu['Population'], 'area':crim['geometry'].intersection(popu['geometry']).area})

      elif popu['geometry'].intersects(crim['geometry']):
        pass
        #print("intersect")

      elif popu['geometry'].disjoint(crim['geometry']):
        pass`


texts: ArcGIS Pro SDK forcing process windows to be minimized?

texts:
I'm running a Process from my add-in, and even if I use P/Invoke to force the window to be open and in the foreground, the process windows that appear will be minimized. Weirdly enough, they are minimized only when using the installed version of the add-in. If I'm debugging in VS, process windows are never minimized, and always get focus.
I start the Process as one normally would,
var startInfo = new ProcessStartInfo(fileName, arguments)
{
    WindowStyle = ProcessWindowStyle.Normal
};
var process = Process.Start(startInfo);

// Here I would put my P/Invoke method calls to try to force-show the window,
// which of course didn't end up working

process?.WaitForExit();

And despite this, the windows from the process will appear, and then quickly become minimized. I've attempted to use P/Invoke methods like ShowWindow(), ShowWindowAsync(), setForegroundWindow(), with process.MainWindowHandle and nothing fixes the issue. I've tried process.WaitForInputIdle() before showing the window to ensure the window has been created before I attempt to show it, but this doesn't work either. It's as if Arc is forcefully minimizing the windows. I've also tried to obtain a Control from the MainWindowHandle, but it always comes out null. It's bizarre that there is no issue when running Arc with a Debug configuration in VS, it's only when I use the add-in as a typical user would.
How can I stop ArcGIS Pro from minimizing the process windows?


texts:
after inspecting the gdalinfo of my input file :

Driver: GTiff/GeoTIFF
Files: elevation.tif
Size is 5555, 2525
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-1.630972222222222,35.957916666666655)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -1.6309722,  35.9579167) (  1d37'51.50"W, 35d57'28.50"N)
Lower Left  (  -1.6309722,  35.2565278) (  1d37'51.50"W, 35d15'23.50"N)
Upper Right (  -0.0879167,  35.9579167) (  0d 5'16.50"W, 35d57'28.50"N)
Lower Right (  -0.0879167,  35.2565278) (  0d 5'16.50"W, 35d15'23.50"N)
Center      (  -0.8594444,  35.6072222) (  0d51'34.00"W, 35d36'26.00"N)
Band 1 Block=5555x1 Type=Int16, ColorInterp=Gray

my guess was that this is because the unit are in degrees so 0.0002 give black or something like...

so i have tried to reproject to utm with

gdalwarp -t_srs "+proj=utm +datum=WGS84+zone=30+units=meter" elevation.tif elevationutm.tif

but gave me an error

  ERROR 1: Too many points (2704 out of 2704) failed to transform,
  unable to compute output bounds.

also i have played a little with gdal_translate with the scale parameter

gdal_translate -scale 0 1 0 100 -of png -ot Uint16 elevation.tif elevation.png

wich gave me something that looks good, but here the 100 value parameter are completly based on some random tests, so my question :

how to scale properly my dem to get the correct png output ?

I want to use these stuff in a terrain 3d engine that take a png Uint16 or raw Int16 as an input.


texts: Duplicate project with different database username & password

texts: I have created a QGIS project that connects to a database using a unique username and password. I need to replicate the project with changes to the username and password. Is there a script that can do that for me?


texts: Saving multiple layers from QGIS to SpatiaLite database?

texts:
Each time I save a layer from QGIS as a spatialite it insists on creating a new db ( if you select a existing db you are asked whether to cancel or replace).

I want to package all the layers of a project into a single spatialite db for distribution.

I figure that I can use native sqlite tools to combine the layers from the individual files but it would be nice not to have to do so.


texts: Join Zip Codes Table to Nearest City from Second Table

texts:
I have two tables, siteareas and ziplookup. The former contains large cities in the US; the latter all US zips with lat/lngs. I have a foreign key into the former table in the latter: sitearea_id. I have made the join for all cities in the zip table that area present in the siteareas table. I wish to do so for the remaining zip codes by finding the nearest city from siteareas and updating ziplookup with that city's sitearea_id.

Here's what I have so far:

update ziplookup set sitearea_id = siteareas.id
FROM siteareas
     inner join ziplookup
     on ST_DWithin(siteareas.lnglat, ziplookup.lnglat, 0.1)
where ziplookup.sitearea_id is null
order by siteareas.lnglat <-> ziplookup.lnglat
limit 1 ;

This is obviously not working, but I need a way to order by the closest result and limit to one result. This is probably more of a SQL question, but any help is appreciated.


texts: How to trigger a second drop-down list when selecting a coded value in ArcGIS?

texts:
I am making a drop down list inside an attribute table of car makes, but want to have a list of models from the respective maker populate when I select a make.

Example:
I select Ford in my CAR_MAKE field
My CAR_MODEL field would populate with models that Ford makes.


texts: I'm trying to allocate a box to install geoServer for my company's intranet site. We have a box with Windows 2012 R2 installed but I can't find any specific documentation if geoServer will run on this OS. Can anyone help me answer this?


texts: GEE Asset Batch Manager credential error

texts:
I am trying to use the Python API to upload batch files to Google Earth Engine Assets. I have successfully installed the GEE Python api found here.

When I try and execute the upload command in the command prompt I get the following error:

C:\Python27>geebam upload -u walteral@oregonstate.edu --source C:\Users\Alex\Desktop\tombali\001 --dest users/walteral/Tombali/001 --nodata 222

2018-08-31 08:14:20,947 :: oauth2client.transport :: INFO :: Attempting refresh to obtain initial access_token
2018-08-31 08:14:20,950 :: oauth2client.client :: INFO :: Refreshing access_token
2018-08-31 08:14:25,348 :: root :: WARNING :: Collection users/walteral/Tombali/001 already exists
2018-08-31 08:14:26,726 :: root :: INFO :: Processing image 1 out of 12:
C:\Users\Alex\Desktop\tombali\001\complexity_mat.tif
2018-08-31 08:14:28,301 :: root :: ERROR :: Incorrect credentials. Probably. If you are sure the credentials are OK, refresh the authentication token. If it did not work report a problem. They might have changed something in the Matrix.

I have tried refreshing the token but that hasn't worked. One issue that I had during the install is that pycrypto would not install, instead I used pycryptodome which was recommended here. I don't know if that is part of the problem or not. Is there some way around this credential error so I can upload my files?


texts: Determine value nearest to given coordinates

texts:
I'm new to GIS.

I have a raster that specifies the geographical distribution of temperature. I would like to find the temperature that is closest to a particular set of coordinates. How do a find the raster value nearest to [x, y]?

I'm using ArcGIS Desktop.

I have numerous points specified by longitude and latitude in a table that was imported using Add X-Y Data... from an Excel spreadsheet. These points plot in the correct places on my map. Nevertheless, the points do not necessarily coincide with data in the raster, so when I do an Extract Values to Points I get some -9999 values. Instead, I'd like to know the value of the nearest point where there is data. The output format doesn't really matter, but a column appended to the original table would be ideal. As far as extensions go, I think I have them all, but I can't say for certain.


texts:
I am trying to process some raster layers from Global Forest Change (https://earthenginepartners.appspot.com/science-2013-global-forest/download_v1.4.html, data produced by Hansen et al., 2013). However, even though I am working in a 32 Gb RAM memory station, everything goes extremely slow.

I am aware that the data I am working with is heavy. Each raster has 40000x40000 raster cells and, while some of rasters only weight some 20 Mb, others go as high as 600 Mb.

Then, procedures such as reclassify() and aggregate() take as long as ~20-30 minutes for the 20 Mb layers. Taking into account that I need to download and process some 2x500 tiles... it is going to take ages.

Is there an efficient way to deal with this kind of data? Is there something more efficient than the raster library (which is awesome, by the way) for this particular workflow?

The rasters are downloaded as separate tiles that form a larger grid spanning the whole world, and there is a separate link for each block (e.g.: the first one is https://storage.googleapis.com/earthenginepartners-hansen/GFC-2016-v1.4/Hansen_GFC-2016-v1.4_treecover2000_00N_000E.tif). This is the part of my code that processes each tile:

library(raster)

down_links <- read.table("https://storage.googleapis.com/earthenginepartners-hansen/GFC-2016-v1.4/treecover2000.txt")

x <- 1
temp_tc <- tempfile()
download.file(as.character(down_links[x,]), destfile = temp_tc) # down_links is a table with the links mentioned above)
tc <- raster(file.path(temp_tc))
unlink(temp_tc)

# Reclassify raster
tc2 <- tc
matrix <- structure(list(from = c(0L, 20L), to = c(20L, 100L), becomes = 0:1), .Names = c("from",
"to", "becomes"), class = "data.frame", row.names = c(NA, -2L
))
tc2 <- reclassify(tc2, matrix)

# Calculate area
layer_area <- area(tc2)
tc3 <- tc2 * layer_area

# Aggregate raster to a larger pixel size
tc4 <- aggregate(tc3, fact = 32, fun = sum)


texts: Batch Clip shapefiles by filename or attribute using Python?

texts: I used a Split Layer By Attribute tool I found online to split a feature class into 90 different shapefiles. These were outputted to a folder with the name for each shapefile ending in a 5 digit unique ID relating to the attribute I split the feature class by. I have done the same thing to another feature class that I need to clip to the first mentioned set of shapefiles, based on the same 5 digit code tied to the end of the shapefile name. How can I batch clip these 90 shapefiles using Python based on the similar file name or similar attribute? Code examples needed because I am new to Python!


texts: How to find fluctuant parts in a trajectory

texts:
I am currently doing an assignment on how to compress a trajectory. Given a trajectory T1 in figure 1, and I want to get a compressed version as T2. For segments p0~p3 and p6~p9 in T2, I can obtain them using traditional line simplification algorithms, such as Douglas-Peucker algorithm. However, in my opinion, segment p3~p6 is very important for further analysis, so, I want to preserve this part while compression T1.

How can I do it?

Segment p3~p6 can be called fluctuant part. It is composited with points that with relative slower speed and continuous direction changes. May be there are any other accurate definitions to describe this part.


texts: Mask a Raster layer in QGIS using a vector layer with multiple polygons

texts:
I would like to extract statistics from a raster file based on a land cover classification vector file, however the land covers (features) are not contiguous (see pic).

In the example below, I would like to calculate a 90th percentile of all pixels in the "12" feature of the overlying vector file. I am looking for a single value, not a 90th percentile for each of the polygons of feature "12".

I can't figure how to do this and get one statistic value for each feature type. Is this even possible in QGIS? If not, is there another way to accomplish this?


texts: Python function which recognises when to convert coordinate system

texts:
How do I write Python code for an AWS lambda function which will convert easting and northings to longitude and latitude?
For example, if a CSV is uploaded into an S3 bucket it may contain either eastings/northings or longitude/latitude columns for the location of different features. I'd like to create a Python code which recognises this and knows what to do... e.g.
if in long/lat format = leave and move on to next function
else:
if in easting/northing then convert to long/lat and move on to next function

So far, I have my working Python code which will convert the eastings and northings to longitude and latitude (see below). Next I just neeed to write a function that recognises when to do this conversion.
from pyproj import Proj, transform

inProj = Proj(init='epsg:27700') #British National Grid
outProj = Proj(init='epsg:4326') #WGS84
x1,y1 = 357551,259776 #(x,y coordinates)
x2,y2 = transform(inProj,outProj,x1,y1) # convert to long/lat
print (y2,x2)

edit
The CSV that is uploaded to the S3 bucket has the following headers for ease of reference.

Feature
Easting
Northing
Length
Date
Comments

Feature1
181016.6
32794.39
6.5 metres
17/01/2020


texts: Oracle equivalent to postgis EPSG 4326 "geometry" type

texts:
We have polygons stored in Oracle SRID 8307 (geodetic datatype), but need to perform various operations on this layer in a 2d cartesian/non-geodetic system based on the geographic coordinates.

I.e. this would correspond to casting from "geography" to "geometry" datatype SRID 4326 in Postgis.

My question basically is: What would the equivalent of "select geom::geometry from ..." (geom being a 'geography' datatype) in postgis be in Oracle? (ref: http://postgis.net/workshops/postgis-intro/geography.html#casting-to-geometry )

I've found that converting the data to a spatial type with SRID set to NULL, works. However this seems like a "duct tape" solution, not exactly elegant, so I am hoping that there is a "clean" alternative?


texts: How to create "two-dimensional" color ramps in QGIS?

texts: Is there a way in QGIS to create a "two-dimensional" color ramp for the graduated symbology; i.e. a color ramp that determines a color based on two different attributes? My situation is that I want to control the hue of the color of a polygon using one attribute and the saturation of the color using another attribute.


texts:
Suppose I have layers A and B.  I don't know their projections but by the looks of the map I know they are not the same.

In QGIS 1.6 I would use one of my favorite features:

Vector> Data Management Tools> Export to New Projection

What was nice about that feature is that you could export B to A's projection by selecting "choose projection from existing layer".  Hence you need not know A or B's projection.

In QGIS 1.7 the Data management Tools menu is much reduced and the Export to New Projection is missing.

Yes, I can do projections on the fly, and I can do save as.  The former is not useful if I want to convert B's projection to A's.  The latter does not give option to set projection from open layer.

I have tried to find different ways to access the Export to New Projection dialog without success.  Any tips?


texts: Show/hide Google Earth folder without changing visibility settings of child items?

texts:
I would like to show/hide a Google Earth folder without affecting the visibility settings of items nested within that folder.

For example, several top-level folders might each contain a hierarchy of folders and placemarks, each of which may have had hand-tuned visibility that I don't want to lose when temporarily toggling visibility of folders higher up in the hierarchy.

It appears that when the visibility of a higher-level folder is toggled, its new visibility setting is forced upon all its children.

Is there a way to preserve the visibility settings of children (ie., make them independent) when their parent's visibility is changed?

[Of course I'd expect if their parent is hidden, a child would be hidden too regardless of its visibility setting.  A child item would display if their parent's visibility is set to 'not hidden' AND their own visibility is set to 'not hidden'.]

If that's not possible with the current GE interface, is there a good alternative or work-around?  RadioFolder style does not seem to be a solution as it seems to be handled by the GE interface the same way as regular folders (children settings are altered when parent visibility is altered).  I guess I could convert folders into network links and then manage them independently, manually doing a 'save' (overwrite) to set child settings and a reload to restore the visibility settings (but that would be clumsy). Any suggestions?


texts: Rasterizing a point layer without using an attribute field

texts:
I am following the "QGIS by Example" text book and a task requires that a point layer is rasterized using the following GDAL code

gdal_rasterize -l layer name -burn 1 -a_nodata 0 -te coordinates - tr 10 10 -ot Byte

This should provide a raster layer with the value of 1 for the point locations and 0 elsewhere within the coordinates provided. A tiff file is generated but doesn't show on the map canvas despite my styling attempts. I have tried all sorts of variations in the GDAL command and added an integer field set to 1 to the shapefile to use this but without success. This may not be relevant but I have used the book download file for the next exercise using the raster proximity tool but receive a message stating that the GDAL dataset can't be opened. Is there an issue with settings in QGIS - my version is 2.18.26


texts: kml to postgis with filename in attribute (ogr2ogr)

texts:
I try to load my kml-s into an existing table (the geometry field is geometry type). One kml file can contains multiple geometry and I would like to add these to an existing table (each feature is a new record) with an ID which should be the kml file's name. I have a partly working code, but I want to add the filename to the ID field as well

ogr2ogr -append -f "PostgreSQL" PG:"host='0.0.0.1' port='5433' user='user' dbname='dbname' password='password'" h:\Users\gpapp\kmls\AEY17525.kml -nln cat_projects_geom

this one generates error:

ogr2ogr -append -f "PostgreSQL" PG:"host='0.0.0.1' port='5433' user='user' dbname='dbname' password='password'" h:\Users\gpapp\kmls\AEY17525.kml -nln cat_projects_geom -sql "UPDATEcat_projects_geom SET nameref = 'AEY17525 or filename'"


texts: Automatic numbering of rows in attribute table fails

texts: I'm working on a big data set, with 385'000 features. Each feature is the centroid of a polygon, and I'd like to assign a value to each of them. I've used the @row_numberin the field calculator to do so (QGIS 3.10). It works perfectly until the row 235'262 and then for not apparent reason, the following row is labelled "NULL". Any insight on what I may be doing wrong?


texts: Adding properties to pop ups in leaflet maps

texts:
I have a leaflet map which display pop-ups that give the name of a country when clicked. How can I add additional properties and have them all displayed, similar to this map http://lotrproject.com/map.

For example, how can I make my Iceland pop-up feature display the properties for "name" and "language", rather than just "name".

var countries = [
        {   type: "Feature",
            properties: { name: "Iceland", language: "Icelandic"},
            geometry: { type: "Point", coordinates: [2535,911] }
        },
        {   type: "Feature",
            properties: { name: "Ireland" },
            geometry: { type: "Point", coordinates: [1324, 1580] }
        },
        {   type: "Feature",
            properties: { name: "England" },
            geometry: { type: "Point", coordinates: [1498, 1662] }
        },
        {   type: "Feature",
            properties: { name: "France" },
            geometry: { type: "Point", coordinates: [1608, 1918] }
        },
        {   type: "Feature",
            properties: { name: "Italia" },
            geometry: { type: "Point", coordinates: [1923, 2093] }
        },
        {   type: "Feature",
            properties: { name: "Hispania" },
            geometry: { type: "Point", coordinates: [1374, 2148] }
        },
    ];

    var layerCountries = L.geoJson(countries, {
        // correctly map the geojson coordinates on the image
        coordsToLatLng: function(coords) {
            return rc.unproject(coords);
        },
        // add a popup content to the marker
        onEachFeature: function (feature, layer) {
            if (feature.properties && feature.properties.name) {
                layer.bindPopup(feature.properties.name);
            }
        },
        pointToLayer: function (feature, latlng) {
            return L.circleMarker(latlng, {
                    radius: 8,
                    fillColor: "#800080",
                    color: "#D107D1",
                    weight: 1,
                    opacity: 1,
                    fillOpacity: 0.8
                });
        }
    });


texts: Finding overlapping areas between two or more GeoTIFF files

texts: How can one programmatically understand overlapping areas between two or more GeoTIFF files captured at various timestamps for the same location?


texts: Original polygon is completely contained by new after simplifying

texts:
Is there a polygon simplification method that can create a polygon, such that the original is completely contained by the new one?

The aim is to create the new polygon from as few points as possible while maintaining a specified maximum tolerance. Until now I have been using a method that creates a buffer to a distance of specified tolerance and then uses a regular simplification method but that seems to be very inefficient as it sometimes creates an unnecessary offset. The original does not intersect itself.

My desired output looks as following


texts:
We have 4 people accessing shapefile data simultaneously in QGIS using a WFS from Geoserver.

We use a file as a container of shapefiles and we all have to work on the same shapefiles at the same time.

But we encounter some issues that we don't understand :

When we modify the geometry or the attributes of multiple entities in the shapefiles, if we don't save the shapefile between EVERY MODIFICATION, the attributes or the geometry is cloned and all the objects that we have modified have now the same geometry or the same attributes (the same as the last modified).

Can somebody explain why this is, and what we can do to fix it? Not sure if it is a bug or a misconfiguration.


texts: What is difference between b3dm and gltf?

texts:
What is difference between b3dm and gltf ?

When to use b3dm and when to prefer gltf over b3dm ?


texts: Seeking QGIS alternative for ArcGIS Schematics

texts:
I want to draw a schematic drawing from a network.There is an extension in ArcMap that can do this.
Is there an extension or a way to do this in QGIS?


texts: Work with individual Polygons of a SpatialPolygonsDataFrame (R)

texts:
I'm trying to work with individual polygons of a SpatialPolygonsDataFrame (SPDF). I've seen these two questions: Using sapply on SPDF and a for-loop on a SPDF. The latter I thought may be helpful, but I don't understand the %>%-syntax and don't know how to apply it to conditional code.
I want to iterate over an SPDF, shrink (gBuffer) polygons, as long as they don't get smaller than 1 ha (10'000 m²), then maybe apply other steps to those polygons. Something like this:
# get shp file with multiple polygons, creates SpatialPolygonsDataFrame
shpA = shapefile(file.path(PATH_DATA, 'sample_test_A.shp'))

# iterate over shp and work with individual polygons
for (i in 1:length(shpA)) {
  # THIS DOESN'T GIVE ME A SPATIAL POLYGON I CAN WORK WITH:
  pol_temp = shpA[[1]] # or shpA@polygons[[i]]

  # create buffered polygon and check if it is > 1ha
  pol_buffered = gBuffer(shpA, width=-30)
  if(area(pol_buffered) > 10000){
    # replace polygon
    shpA[[1]] == pol_buffered
  }

  # maybe do some other stuff, like sampling the polygon
  ...
}


texts: Misplacement while adding shapefile with no projection

texts:
I'm stuck trying to add a shapefile of points (digitised from geotagged aerial pictures, I'm assuming GPS tagged) in QGIS 3.4 (Mac). The .prj-file that I opened in notepad gives me a WGS84 GCS.

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

But when I add the shapefile in QGIS, and I click on the info button to see the "coordinates" it's in metrics. The layout of points seems correct just their.

The area is in the Mackenzie Delta in the NorthWest Territories, Canada. I tried pseudo mercator (EPSG:3857), conic Lambert (EPSG:3978 and EPSG:3581), UTM 8 (EPSG:26908), the values are just off (it's like they just "slid").

Any idea how to fix this so the points would be at their actual accurate locations?


texts:
I have several thousand small glacier rasters (in rectangular shape) that need to be clipped to the glacier's extent. I have the glacier polygons stored in a large SpatialPolygonsDataFrame (matching in length and position).

I need parallel computing or else it would take forever, so I would like to use the rslurm package.

My code so far (using the for-loop provided here: https://gis.stackexchange.com/a/342955/145570)

library(rslurm)

fun_clip_rslurm <- function(filelocations, filenames, glaciers){
  for(i in 1:length(filelocations)){
    r = raster(filelocations[i]) # get the ith file name
    g = glaciers[i,] # get the i'th row of `glaciers`.
    r = crop(r, g) # crop
    writeRaster(r, file.path("/data/.../RGI60-13_crop/", filenames[i]))
  }
}

params <- list(filelocations, filenames, glaciers)
names(params) <- c("filelocations", "filenames", "glaciers")

slurm_call(fun_clip_rslurm, params, jobname = "my_crop",
           slurm_options = list(time="1-00:00:00", share =TRUE),
           submit = TRUE)

If I run this, I get the following error message:

sbatch: error: Batch job submission failed: Invalid account or account/partition combination specified
$jobname
[1] "my_crop"

$nodes
[1] 1

attr(,"class")
[1] "slurm_job"

So now I have two questions:
Do I have to give rslurm some kind of permission to run the batch-script on the cluster?
And how can I specify the number of nodes on which to run the calculation?


texts:
QGIS 1.9

Is it possible to destroy a spatial index?  I am noticing that if I generate a spatial index on a set of polylines (within a plugin), split one of the lines, and then recreate the spatial index, the index values within the spatial index are doubled.

For example:

I have 5 lines, so the spatial index has feature 0 - 4.  I split one of the lines and recreate the spatial index.  The spatial index should have indices 0-5, but it has indices 0-9.

Deleting the features from the spatial index does not help with the indexing issue.

Code

class slicer:

    def init(self):
        '''standard stuff here'''
    def activated(self):
        '''standard stuff here'''
        create_temp_polylines()

    def create_temp_polylines(self):
        self.line_layer = QgsVectorLayer('Make the layer')
        self.spat_index = QgsSpatialIndex()
        self.point_layer = point_layer #Grabbed by iterating over layer list
        pt_prov = self.point_layer.dataProvider()
        polylines = {}
        for feat in self.point_layer.dataProvider():
            '''build the line geometries based on a membership attribute'''
        for key, value in polyline.iteritems():
            temp_polyline = QgsFeature()
            temp_polyline.setGeometry(value)
            temp_polyline.setAttributes([attr1, attr2, ...])
            line_layer.startEditing()
            line_layer.addFeature(temp_polyline, True)
            self.spat_index.insertFeature(temp_polyline)
        line_layer.commitChanges()

    def split(self, cut_line):
        '''logic to split a polyline'''
        '''logic to update the attribute table'''
        '''commit the changes to the DB'''
        QgsMapLayerRegistry.instance().removeMapLayer(self.line_layer.id())
        self.create_temp_polylines()

Spatial_Reference Checking Code

'''draw a line on screen (rubberband)'''

def get_intersect(self, rubberband):
    line = rubberband.asGeometry().asPolyline()
    intersections = self.spat_index.intersects(rubberband.asGeometry.boundingBox())
    touches = set()
    for id in intersections:
    id = (id * -1) -1 #Why are they ids negative longs?
    '''convert if to feature object'''
    '''check touches'''
    if touching:
        touches.add(feature.id())


texts: SpatialPointsDataFrame properties and operators in R

texts:
I have created an object of type SpatialPointsDataFrame using the sp package in R. However, I am confused about the @, $, . and [] operators and when to use them to access the different properties of my object. Here is my sample code:

library(sp)
library(rgdal)

#creating a SpatialPointsDataFrame with sample points in UTM
x <- c(15.2, 15.3, 15.4, 15.5, 15.7)
y <- c(50.4, 50.2, 50.3, 50.1, 50.4)
v1 <- c(1.0, 2.0, 3.0, 4.0, 5.0)
v2 <- c("a","b","b","c","a")
attributes <- as.data.frame(cbind(v1,v2))
xy <- cbind(x,y)
locationsDD <- SpatialPointsDataFrame(xy, attributes)
proj4string(locationsDD) <- CRS("+proj=longlat")
locations <- spTransform(locationsDD, CRS("+proj=utm +zone=33"))
plot(locations)

#using the different operators: WHEN TO USE @, $ or [] ?

#all these work!
property1 <- locations$v1
property2 <- locations@data$v1
property3 <- locations@data[,"v1"]
property4 <- locations@data["v1"]

#these also work
property5 <- locations@coords
property6 <- locations@bbox
property7 <- locations@coords[,2]

#these three work only in my special case
property8 <- locations@coords[,"y"]
property9 <- locations$x
property10 <- locations$y

#these don't work: $ operator is invalid for atomic vectors
property11 <- locations@coords$x
property12 <- locations@coords$y

Could anybody help me, when to use the @, $, [] operators? When I try to read the documentation ?SpatialPointsDataFrame I can see the different properties such as coords or bbox but I'm confused which operator @, $, [] to use to access them or modify them.


texts: Setting metadata options for ouput GeoTIFF when using gdalwarp?

texts:
I'm reprojecting and changing the format of a vrt to GeoTIFF and can't seem to control metadata options with gdalwarp.

Other tools like gdal_translate have the -mo option to specify metadata options but gdalwarp doesn't.

In this case, the original .dem (CDED format) have a metadata option that says AREA_OR_POINT = POINT and the output GeoTIFF, using gdalwarp, defaults to AREA_OR_POINT = AREA.

To be honest, I don't know the significance or importance of the AREA_OR_POINT metadata option, but just like matching input to output as much as I can.

The command I'm using is:
gdalwarp -of GTiff \
          -ot Int16 \
          -s_srs "EPSG:4269" -t_srs "EPSG:3005" \
          -tr 15 15 -r bilinear \
          -wm 4096 -multi \
          -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9 \
          -co NUM_THREADS=ALL_CPUS \
          -srcnodata "-32767" \
          mosaic_dem_cded_4269.vrt dem_3005.tiff


texts: Make an object move along a route given by an OSMR RouteStep object

texts:
I want to create an agent based simulation of traffic. I have an object for a car (an agent) that has a starting point and a destination. Using the OSMR API, I request a route between the two locations. I store this route as an attribute of the agent.

At each time step, I want to update the coordinates of the agent as it moves toward its destination along the route.I cannot figure out how to use the response of the API to do this.

The response includes a RouteStep object which describes the steps needed to reach the destination. One step looks like this:

{'intersections': [{'out': 0, 'entry': [True], 'bearings': [344], 'location': [24.707075, 60.157693]}],
'driving_side': 'right', 'geometry': 'qpdnJgrxuCg@VoAv@', 'mode': 'driving', 'maneuver':
{'bearing_after': 344, 'bearing_before': 0, 'location': [24.707075, 60.157693],
'modifier': 'left', 'type': 'depart'}, 'weight': 46.8, 'duration': 30.7, 'name': '', 'distance': 70.3}

Question: How can I make the agent move (given its speed) along the route described in these steps?

Any advice would be appreciated, even just high-level ideas about how to make it work conceptually. This is supposed to be more of a proof of concept so it does not have to work perfectly at all times.


texts: Running Viewshed on one feature at time in shapefile?

texts:
I have over 100 features in a shapefile and I want to create individual viewsheds for each and then run processes on each of those resulting viewsheds in ArcGIS using Python.

Is it possible to write a script using a cursor to run through each "row" and use that individual observer point to run the viewshed?

Then update the cursor to the next row and rerun it until all features have been exhausted?

Or would I have to first Split the shapefile into over 100 files and run each individually?

Also, is it possible to designate the filenames for the resulting viewsheds (or shapefiles) using a specific field?


texts: How to check for the nodata cells in a raster image

texts: Please let me know how to check for the locations of nodata cells in raster images  in ArcMap. The raster images has no attribute table.


texts:
Sorry to ask, but I just don't have time for help menus this week.  That said, does anyone know if CAD is capable of re-projecting shapefiles?  I find myself in a situation whereby I may have to use spatial data without access to ESRI software.  The data are delivered in the GCS_North_American_1983 datum, and my clients site drawing is in NAD83 Z17N.

I fully expect problems.  Can any of you post any related links that you're aware of?  There's no budget to purchase software on this one either (just a heads-up).

Alternatively, is there a way to get CAD to project on the fly?


texts:
For simplicty, referencing this example:

https://carto.com/developers/carto-vl/examples/#example-load-external-geojson-layer

I would like to have this bit of code behave the same as through it were an in-line GeoJSON declaration

fetch('https://libs.cartocdn.com/carto-vl/assets/stations.geojson')
  .then(response => response.json())
  .then(function(data){
      // Define layer
      const source = new carto.source.GeoJSON(data);
      const viz = new carto.Viz();
      const layer = new carto.Layer('layer', source, viz);

      layer.addTo(map, 'watername_ocean');
      layer.on('loaded', hideLoader);
  });

For example (this code works, although for brevity only includes one coordinate):

const stations =
{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "name": "Van Dorn Street",
        "marker-color": "#0000ff",
        "marker-symbol": "rail-metro",
        "line": "blue"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -77.12911152370515,
          38.79930767201779
        ]
      }
    },
  ]
};

const source = new carto.source.GeoJSON(stations);
const viz = new carto.Viz();
const layer = new carto.Layer('layer', source, viz);
layer.addTo(map, 'watername_ocean');

The above is an in-line example which avoids issues relating to global variables versus a variable defined within a function.

What I'm trying to figure out to no avail is how to do  use fetch to define a global variable, with the same ease as defining it in-line.

This code does not work:

let stations;
fetch('https://libs.cartocdn.com/carto-vl/assets/stations.geojson').then(
        function(u){ return u.json();}
      ).then(
        function(json){
          stations = json;
        }
      )

const source = new carto.source.GeoJSON(stations);
const viz = new carto.Viz();
const layer = new carto.Layer('layer', source, viz);
layer.addTo(map, 'watername_ocean');

Generates the following errors:

  Error: [Missing required property] 'data' carto-vl.min.js:6:3429


texts: problem saving GRID to .img keeping original names (export by attributes)?

texts:
I am really new in Python programming. I try to make a loop for Extract by Attributes tool in ArcGIS 10.1. I have a huge dataset of raster (760pc) and I want to extract a specific value. This part of script runs fine. But ArcGIS in this function can only produce a GRID format with 13 characters so it changes my original names (.img file).

How can I keep the original names of rasters from my datasets?

Maybe I have to produce a list of results and convert them to .img file?

When I did it, it has kept the new generated GRID name (Extract_1), not my original name.

Anyone can help?

# Name: ExtractByAttributes_Ex_02.py
# Description: Extracts the cells of a raster based on a logical query.
# Requirements: Spatial Analyst Extension

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = r"d:\UEL\Data\poskod_cez_VI_2\ch_msi_vci"

# Set local variables
inRaster = "ch_n_msi_00_02_10.img"
inSQLClause = "VALUE = 5"
outFolder = r"d:\UEL\Data\poskod_cez_VI_2\export_values"

#Allows overwriting
arcpy.env.overwriteOutput = True

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Get a list of the rasters in the workspace
rasters = arcpy.ListRasters()

# Loop through the list of rasters
for inRaster in rasters:
    # Set the outputname for each output to be the same as the input
    outRaster = outFolder + "\\" + inRaster

    # Process: Extract by Mask
    attExtract = ExtractByAttributes(inRaster, inSQLClause)
# Save the output
attExtract.save(r"d:\UEL\Data\poskod_cez_VI_2\export_values")


texts: Replacing blank values in a field using Calculate Field in ArcGIS ModelBuilder

texts:
I'm using ArcGIS 10.5 and working with geodatabase layers.

My previous step was only a union between my area of study layer (which has no info in its table) and another layer with the soil cover. The resulting layer took the COVER field of the second layer, but the parts where both layers did not overlap are now polygons in which the COVER field is blank, and I need to change it to "Urban".

I could update them manually by selecting and using Field Calculator, but I'd like to be able to do this step inside ModelBuilder so it is automated. I found the sequences to replace strings using Calculate Field, both for VBSCript and Python, but what gives me trouble is the part where I tell the tool what to find, the ones below with the xxxx:

replace([COVER], "xxxx", "Urban")

!COVER!.replace("xxxx", "Urban")

What I don't know is what iI should type to tell the tool that it should look only for blank cells in that particular field (there is not a single NULL). I tried putting a space between the ", not putting anything at all, just writing blank, NULL, typing a "-", etc, but none worked. Many simply don't change anything, and other such as typing a space change the spaces between the words in the rest of my other cells (but not the blank ones), which that should not be changed.

My current sequence (copy/pasted a it is) is the following, maybe something else is wrong?

!COVER!.replace("", "Urban")

What should I write in the formula so that it replaces only the blank cells?


texts:
I am quite new with QGIS and I created some maps with different layers these last weeks (seems easy but took me hours to do so). As usual, I saved everything in a folder yesterday, but today, I couldn't open my project or the shapefiles. It says "datasource is not recognized..." and, as I checked, I realized the .shx files are all missing.
I read in one of the discussions here on Stack Exchange that it may be caused by an antivirus "cleaning" some folders (I have Avira).

Is it possible that Avira has deleted these files?

How can I get them back and avoid that this happens again next time?


texts: Reload Table Cache with Python & ArcGIS 10.1

texts:
I am trying to create a button in arcmap that reloads the table cache. Exactly, like the existing button Table Options > Reload Cache. The button is inside when you open an attribute table. I would like to create a custom button that is in a toolbar I made. I tried to search for tools in customization and it isn't present nor could I find any sample code for this.

I have a custom Python Addin Toolbar with comboboxes that manipulates data in a database and since it doesn't update the data upon selection I would like another way to refresh with having to close the table and reopen and//or selecting the existing tool.

Thanks for your help.


texts: Is it better to use QGIS LTR or current release version

texts:
When downloading QGIS you can download the Current Release or the LTR Version. But what are the main differences between them? Here this question has been asked before but I don't really understand the answers. In the most voted answer the user wrote:

QGIS (long-term release) can be considered static (no new features added) but with cherry-picked patches that are back ported for
stability, enhancement and security.
QGIS (current release) is mainly different from the LTR, that in
addition to back-ported patches, it'll get new features added (from
master version) when they are considered stable enough with frequent
releases (proposed every 3 months).

I understand that the Current Release has the advantage that it has more tools and features, but what are exactly the advantages of the LTR Version? Does it have better "back-ported patches"? What are back-ported patches practically?


texts: Union of shapefiles to generate composite shapefile with overlap frequency count for each polygon

texts:
I'm new to QGIS, making the switch from ArcMap. I have 5 shapefiles of polygons representing 5 years of insect defoliation. I need to overlay the 5 years and create a new composite shapefile that is a union of all 5.
There are some areas where the polygons overlap from year to year and I want the polygons of the new composite shapefile to have an attribute field "Number_of_Yrs" which represent the total number of years of defoliation. The values of this field should be 1 to 5 years.


texts: Dynamically have symbology change as a due date nears?

texts: I need to make a countdown based on symbolization.  Let's say I have several groups of sections with a Date field in it.  I have when the section is due in the date field.  Can I get it to dynamically hatch close lines as the date nears?


texts: Map statistical significance in the difference of two raster maps

texts:
I have two raster maps that are the output of the same model (they describe the same variable) but ran with different parameters and I would like to assess if there is a statistically significance change in the two maps. How could I do that applying a statistical test and mapping also the pixels where the differences are statistically significant, like it is shown on the map below:
https://phys.org/news/2016-03-raincloud-climate.html

I would like to do this in r. Thus far, I have tried to turn my two rasters as data frames and then calculate perform the wilcoxon test or a chisq test, but that results to getting simply the values of the test, such as "df" and "p-value". How could I produce a map such as the one displayed in the link above?


texts: Identifying Parcels on Street Junction in QGIS

texts:
I have the following shapefiles (click here to download them):

parcels polygon
parcels point
streets/row polygon
streets line

I use QGIS to manually identify polygon parcels with street corners (street junctions) on them. I would like to know if there is a plugin or algorithm to provide such a solution to identify all parcel polygons that are on street junctions. That is those polygons that are on the corner of the street line.

For example, the parcels mark in the red box below have corners. So, they are selected for export.

This could be done a manual visual inspection of the parcels. However, I am looking for an automated in QGIS.


texts: How to find the highest points along an embankment and form a ridge polyline- Methodology in ArcGIS 10.3

texts:
I am analysing data from a series of flood embankments which run parellel to a river. I have 80 KM of flood embankment and laserscanning data covering the area with 50cm resolution from which I have created a TIF DEM.

For a web-app project I am required to deliver a polyline and point feature class.These will show the highest points along the embankment with a distance of 10 meters between each point as well as a polyline showing the rigdge of the embankment parellel to the course of the river.

I have 3 methods which could work.

digitize the ridge of the embankment per hand using aerial photos and then convert the polyline to points with a distance of 10 meters between each point. Using ArcGIS Spatial Analyst I could then extract height from the DEM and convert to Points.

The only problem with this method is that if I miss the summit of the line, I am not extracting the highest value.

An alternative workflow would be to use the polyline which estimates the ridge of the embankment and convert to points with distance of 10m between each point. Then using ETGeowizards, create "station lines" (ie a perpendicular line with its center on each point. Using Spatial analyst again, I could then extract the highest value along each station line and thus determine the summit value.
The third method would be along the lines of the folowing thread.(How do you identify ridges on an elevation raster?). This defines the ridge using ArcGIS focal-statistics and the raster calculator. My problem here is the fact that the result is not a point file.

Could anyone point to a better solution or advise me on which method to follow using the ArcGIS Platform?


texts:
I've got a regular polygon file that looks like this:

Unfortunately, this file stems from a v.overlay so that the 'cat' values are are a bit messed up. I would like them e.g. to start in the upper left corner and increase from left to right and top to bottom.

How is that possible?
I already tried renaming and updating the columns, but somehow this is never allowed as is touches (of course) the 'cat' column.


texts: Using ArcGIS Pro to select all data collected "Yesterday" regardless of the time of day

texts:
I am using ArcGIS Pro 2.4.3.
I can upgrade if needed.

I am reading a hosted feature layer from my Organisation's AGOL site in ArcGIS Pro.  This is using the standard Esri Data Store.  This dataset has a date field "CreatedOn" in the format dd/mm/yyyy hh:mm:ss

I want to use the "Select Layer by Attribute" tool to select all records that were captured yesterday.

I need the selection to be "yesterday" rather than a specific date as I want to use the tool in a model in ModelBuilder that will run daily, hence I need this to be a dynamic selection.

I have so far tried the following:

CreatedOn = CURRENT_DATE() -1

This returns all records collected yesterday at the exact time of day that I hit run on the query.

CreatedOn = CURRENT_TIME()- 1

This returns all records collected yesterday at the exact time of day that I hit run on the query.

CreatedOn = CURRENT_TIMESTAMP()- 1

This returns all records collected yesterday at the exact time of day that I hit run on the query.

CreatedOn BETWEEN CURRENT_DATE -1 AND CURRENT_DATE

This returns all records collected in the last 24 hours, not just yesterday. i.e.  If I run this query at 09:37 and a record was collected yesterday at 09:35, it won't return that record.  This will also return everything collected today.

CreatedOn BETWEEN CURRENT_TIME -1 AND CURRENT_TIME

This returns all records collected in the lat 24 hours, not just yesterday. i.e.  If I run this query at 09:37 and a record was collected yesterday at 09:35, it won't return that record.  This will also return everything collected today.

CreatedOn BETWEEN CURRENT_TIMESTAMP -1 AND CURRENT_TIMESTAMP

This returns all records collected in the lat 24 hours, not just yesterday. i.e.  If I run this query at 09:37 and a record was collected yesterday at 09:35, it won't return that record.  This will also return everything collected today.

CreatedOn BETWEEN CURRENT_DATE -1 AND CURRENT_DATE -1

The expression is valid but won't return any records

CreatedOn BETWEEN CURRENT_TIME -1 AND CURRENT_TIME -1

The expression is valid but won't return any records

CreatedOn BETWEEN CURRENT_TIMESTAMP -1 AND CURRENT_TIMESTAMP -1

The expression is valid but won't return any records

It seems that regardless of whether you use CURRENT_DATE, CURRENT_TIME or CURRENT_TIMESTAMP it always takes the time of day the record was collected into account and so either returns what was collected exactly 24 hours ago or what has been collected in the last 24 hours depending on the query.

Is there a way, without calculating an extra field, to truncate the date field on the fly within the query so that it always returns everything that was captured "Yesterday" regardless of the time?

Or is there a different tool I should be using in my model to do this?

My ultimate aim is to have a simple model that can be run once a day that will export all the features from the hosted feature layer (that were added the previous day) into a file Geodatabase.  Whether this is into an existing fGDB or a new fGDB every day doesn't really matter.


texts: Why does my threejs map render outside of the visible extent of the browser?

texts:
I created a 3d visualization using qgis2threejs but the map renders outside the visible extent of the browser window so I have to pull it down into the window each time.

Is there a way to change the rendering height in the .js file or something like that?  If so can you explain what parameter I would need to change?


texts: How to get the geometry value from coordinates

texts:
I have a table with empty geometry column. I have two columns in the same table which has xcoordinate and ycoordinate values. I want to get the geometry value formed by these two coordinates. The geometry is POINT so I am using st_geomfromtext('POINT(x y)') function.

I tried below syntax to update the geometry column the_geom in table 't1':

 update t1 set the_geom=st_geomfromtext('POINT('||t1.xcoordinate||t1.ycoordinate||' )');

However this gives an error that :

"POINT(3458196941 )" <-- parse error at position 17 within geometry

I am not using '||' string concatenation correctly here. Can any one correct me ?


texts: Pairing line segments that follow a similar path

texts:
I have a shapefile with a regional travel model comprised of straight line features (model.shp, blue in the image). These lines roughly match the actual paths of the roadways they represent (road.shp, magenta in the image). I need to copy traffic volume data from road.shp to the model.shp.

Some of the issues I see right away are:

There is no attribute data (roadway name, segment, intersection, etc)
to compare between the two datasets. This means the comparison will
have to be done based on their spatial relationship to one another
Single lines in model.shp are often comprised of
multiple segments in road.shp. This means that many segments
in model.shp will have multiple attribute matches (one
for each segment in road.shp that falls along the path of a feature in model.shp), and the locations of the beginning and end of
each segment will sometimes be very different.
The paths of the lines are not very close to each other. This means that there is no
simple spatial query that can be done to select by the relationship
of their locations to one another.

The locations of the intersections is the areas that are the closest between the two datasets, and I think that is the best way to go about linking the two together. My initial thought is that the best way to go about this is to make a tool to:

Calculate the coordinates of the beginning and end of each line in road.shp.
Create a straight line segment from that coordinate information in a new shapefile (roadstraight.shp, carrying over the traffic count number we want to copy to model.shp
This should result in a dataset comprised of straight lines which should be closer to the model.shp data
For each line segment in roadstraight.shp, conduct a spatial search for each model.shp line segment within X feet (or the closest X number of lines)
Calculate the slope of each line that was selected to determine which line follows the closest path to the segment in question

Does it seem like I am on the right track with this?
tldr: need help to match lines from one shapefile to another based on how similar they are in shape and path.


texts: Editing geojson without separating segment endpoints

texts:
I have a few maps of city sidewalks in Geojson format which I need to manually edit to make corrections. However, in the tools that I've tried, geojson.io and QGIS, segments that share an endpoint are not treated as connected but rather as four separate segments. So when I drag a shared endpoint, it becomes multiple separate endpoints as shown below:

Original map:

After dragging the bottom left point, it separates into four:

I don't want this to happen; instead I need all the segments to remain connected if I drag an endpoint. Is there any existing tool I can use to do this? Or an easy way to program a custom tool?

I realize this may be a limitation of the geojson format itself, since each segment is just represented as a start coordinate and an end coordinate. But what I'm imagining is a tool that can recognize when multiple segments share the same endpoint, so if that point is dragged, the coordinates are automatically changed for all segments that originally had that endpoint.


texts: What's the difference between Integer 16bit, 32bit and 64bit in QGIS?

texts: When creating a new attribute field in QGIS, these three options are offered for integer numbers but what is the difference between Integer 16bit, 32bit and 64bit in QGIS?


texts:
Using the QSpatiaLite plugin in QGIS I composed the following query

SELECT  'AddressPoint_2'.'buildingNu', 'AddressPoint_2'.'postCode' WHERE  'AddressPoint_2'.'postCode' = 'NR13 3AA'  FROM  'AddressPoint_2'

This returns an error:

"The SQL query seems to be invalid near "FROM":syntax error."

A previous query:

SELECT   'AddressPoint_2'.'buildingNu', 'AddressPoint_2'.'postCode' FROM 'AddressPoint_2'

runs OK. I cannot see why the slightly more complicated query is failing on the same use of FROM

Any pointers please?


texts: How to make the Polygon tool the active tool when the Sketch widget opens?

texts:
Using the ArcGIS Server JavaScript API v4.18, I am adding the Sketch widget. My goal is to have the Polygon tool as the only available tool, and to have this tool active when the widget opens.
I've been able to switch off all other tools using the visibleElements option:
visibleElements: {
  createTools: {
    point: false,
    circle: false,
    polygon: true,
    polyline: false,
    rectangle: false
  },
  selectionTools:{
    "lasso-selection": false,
    "rectangle-selection": false,
    "pointer-selection": false /* this has no effect - I don't know the correct name for this tool */
  },
  undoRedoMenu: false
}

However, the standard pointer select tool is shown by default, and is the active tool...

... and I need to click on the Polygon tool to activate it. How can I:
(a) not show the pointer select tool
(b) have the polygon drawing tool active by default?
The option activeTool: 'polygon' does select the polygon tool by default, but not correctly - the select tool is simultaneously selected, so neither tool works as expected:

See the demo at https://codepen.io/slead/pen/ZEBNpbd, based on the Esri Sketch Widget demo.


texts: Python scripting (update cursor, extract previous row values)

texts:
I am new to Python scripting.
I would like to use SearchCursor to iterate in each row using the FID as index. Then I would like to evaluate length based on the following conditions:

If the length >100 then (CAT remains at 0) and proceed to next row.
If the next length is less than 100, then it should go back to the previous (immediate) row (that has a value of >100)  and extract values in that row (FID, Length, CAT), and save them in an empty table, then the script should stop there.

fields = ["FID", "Length","CAT"]

with arcpy.da.SearchCursor(fc, fields) as cursor:

    for row in cursor:
        if row[1] > 100:
           print (row[2])

        elif row[1] < 100:
           #Some code accessing the previous row here
           print (row[0], row[1], row[2]) #to test what is to be printed


texts:
I have installed all most every requirements to configure Mod Tiles. When I am trying to install it by ./autogen.sh its not installing. Its replying

./autogen.sh: 2: ./autogen.sh: autoreconf: not found

Can any body tell me whats happen?


texts: Inconsistency in areas extracted using Tabulate Area tool

texts:
I am extracting data using the tabulate area tool. The data has too many duplicates for each plot. However, when extracting these plots with duplicates separately, I do not see the duplicate values. In other words, I get different values. With respect to duplicate values, these different plots cannot surely have same values at precision 6 decimal places. Which means something is wrong with my batch area extraction of classes by plot.

The general steps I used are as follows:

Classified my image into three class, Spatial reference: WGS_1984 (Output file=MLClass.tif)
Created polygons over the image, with project coordinate system (WGS_1984_UTM_Zone_36S)  (Output file=Field Plots.shp)
Convert Field Plots.shp to raster using: Feature to raster/polygon to raster tools Environmental settings [Snap Raster (MLClass.tif)] (Output file=Feature_Field_Plots)
Tabulate Area (Spatial Analyst)>> Input zone=Feature_Field_Plots; input Class=MLClass.tif. (I selected attributes of interest in each case) (Output file=Tabulate_Feature)

Can someone help figure what is it I am doing wrong?


texts: Is Output Cell Size Environment the same as Resample tool?

texts:
I have a raster catalog that contains hundreds of rasters (10 x 10). I need to perform Cell Statistics SUM on all the raster. End user able to specify under what cell size the want to perform the operation (ex: 50 x 50).

I want to find out if the following two scenario is actually the same?

1st: Set the CELLSIZE Environment variable of Cell Statistics tool to 50 x 50 and perform cell statistics SUM.

2nd: For each of the raster, run Resample (type: NEAREST) tool to resize to 50 x 50 and run cell statistics SUM on all the rasters.


texts: Missing icon within legend if layer has a complex styling with QGIS

texts:
I'm using a complex styling for a points layer in QGIS, including two markers and a Geometry generator to set up a line between points and labels (see style description in this post).

The result is that either in the Layers Panel and within a legend created with the map composer, the layer has no representation (no icon at all!).

Layers Panel (see either obs_tri and site_tri layers) :

Map composer's legend :

Is this a bug or is there a way to get a nice legend all the same?


texts: Is it possible to make offline public transport routing? What services should I use and what steps are? I'm developing in iOS platform.


texts: Exploratory Spatial Regression: How to limit 60 variable candidates?

texts:
I want to perform an Ordinary Least Squares (OLS) regression (global) or even a Geographically Weighted (GWR) regression (local) to identify possible causes for specific phenomena, like crime.

I have around 60 independent variable candidates, that might be explanatory for the dependent variable.

To identify which of these variable candidates to use for in an OLS or GWR model, I use the tool "Exploratory Regression" within ArcGIS. The tool allows to Input multiple variable candidates and model criteria, tries all different combinations and lists the results, which can then be compared based on the AICc, R2 adjusted, VIF, etc.

I am aware how to use the tool, how to interpret the result and how the overall procedure works.

What I am unsure is, how to limit my huge number of variable candidates?
The calculations are either too time consuming or not commutable at all.

One approach, suggested on the ESRI Website is to start with a low minimum and maximum number of candidates to be tested, e.g. starting with 1, then 2, and so on. When I run it with 1 candidate to be tested, I get a significance of 100% or 0% for each variable candidate. If it is done for 2, I get many different percentage values. If those lower e.g. 50% are removed, should I run the model again for 2 tested candidates or already for 3? There appear to be different results. Is there any threshold that could be used?


texts:
I'm attempting to place a custom marker over a point in my ArcMap add in on my map but nothing shows up on the map. What am I doing wrong? The file extracts correctly and I can view it just fine if I open the .emf but the image doesn't appear over my point or even appear at all.

    internal static void IndicatePoint(IGeometry shape)
    {
        //Get and clear the graphics layer
        var graphics = ArcMap.Editor.Map.ActiveGraphicsLayer as IGraphicsContainer;
        graphics.Reset();
        graphics.DeleteAllElements();

        if (shape != null)
        {

            //Compose the symbology
            var symbol = new PictureMarkerSymbol();

            var extractedPath = System.IO.Path.Combine(System.IO.Path.GetDirectoryName(Assembly.GetExecutingAssembly().Location), "SelectionSymbol.emf");

            if (!File.Exists(extractedPath))
            {
                var path = System.IO.Path.Combine(Environment.GetFolderPath(Environment.SpecialFolder.MyDocuments), @"ArcGIS\Addins\Desktop10.4", ThisAddIn.AddInID, @"xivictoolbar.esriaddin");
                using (ZipArchive zip = ZipFile.OpenRead(path))
                {
                    ZipArchiveEntry zipEntry = zip.GetEntry(@"Images\SelectionSymbol.emf");
                    MemoryStream memoryStream = new MemoryStream();
                    using (Stream zipStream = zipEntry.Open())
                    {
                        zipStream.CopyTo(memoryStream);
                    }

                    var bytes = new byte[memoryStream.Length];
                    memoryStream.Position = 0;
                    memoryStream.Read(bytes, 0, (int)memoryStream.Length);

                    File.WriteAllBytes(extractedPath, bytes);
                }
            }

            symbol.CreateMarkerSymbolFromFile(esriIPictureType.esriIPictureEMF, extractedPath);
            symbol.Size = 24;
            symbol.Angle = 0;
            symbol.XOffset = 0;
            symbol.YOffset = 0;
            symbol.Color = new RgbColor() { Red = 255, Green = 255, Blue = 100 };

            //Compose the marker
            var marker = new MarkerElement();
            ((IMarkerElement)marker).Symbol = symbol;
            marker.Geometry = shape;

            //Add marker to map
            graphics.AddElement(marker, 0);
    }


texts: Can single raster pixel hold more than one value?

texts: Can a single pixel of a raster hold more than one value?


texts: How to interpret 6S atmospheric correction results?

texts:
i apply 6S atmospheric correction on Landsat 8 imagery. My study area is in the ocean. After applying 6S atmospheric correction, there is a negative value in my result instead of being between 0 to 1.

Is this result is good or do I need to retry the correction?


texts:
Is there a simple way to automate calls to QGIS plugins from a QGIS Python script? Either from console or a plugin?

I am ideally looking for a way to specify input/output parameters and bypass GUI.

Specifically, I'm trying to use the SimpleSVG plugin in a script that generates a whole bunch of maps.


texts: What could cause issues migrating a Python plugin to another Windows version?

texts:
I am developing a plugin in my laptop and the plugin runs well. So, I installed the same version in other computer and the plugin doesn't run well. Sometimes I develop the code in laptop, sometimes in the other computer. The same language (python), the same version. The windows is different, one is vista and the other is windows 7. The difference can be here?

Can someaone explain me why?


texts: How do I create a raster file from a geopandas dataframe storing lats, longs, values, and geometries

texts:
From a text file with time, lat, long, and cf values, I create a pandas dataframe (cf_df). Then I use the lat, long values to create a geoseries and create a geopandas dataframe from the original pandas dataframe and the geoseries:

#creating geometries
geometry = [geom.Point(xy) for xy in zip(cf_df['long'], cf_df['lat'])]

# creating geodataframe with original cf and new geometries
geo_cf_df = GeoDataFrame(cf_df, crs=crs, geometry=geometry)

Ending up with a geopandas dataframe like this:

        time    lat    long   cf              geometry
0          1  15.25  200.25  NaN  POINT (200.25 15.25)
1          1  15.25  200.75  NaN  POINT (200.75 15.25)
2          1  15.25  201.25  NaN  POINT (201.25 15.25)
3          1  15.25  201.75  NaN  POINT (201.75 15.25)
4          1  15.25  202.25  NaN  POINT (202.25 15.25)
5          1  15.25  202.75  NaN  POINT (202.75 15.25)
6          1  15.25  203.25  NaN  POINT (203.25 15.25)
7          1  15.25  203.75  NaN  POINT (203.75 15.25)
8          1  15.25  204.25  NaN  POINT (204.25 15.25)
9          1  15.25  204.75  NaN  POINT (204.75 15.25)
10         1  15.25  205.25  NaN  POINT (205.25 15.25)
...      ...    ...     ...  ...                   ...

What methods are available for creating a ASCII-based raster that could be imported to ArcMap?


texts: How to create a style file for QGIS from a plain text file?

texts:
Using QGIS 1.7.

I have a plain text file, which lists 115 rgb values against a map code for a geological map. Viz.:

$ head gsv1Msymbology.txt
MAPCODE RED GREEN   BLUE
Oc  143 255 255
WAT 255 255 255
Qa  244 250 202
Qdl 195 239 218
Na  248 255 238
Qd2 227 255 190
Qxw 248 255 238
Qns 255 148 83
Qn  255 202 190
....

Is there some tool to generate a QGIS style file from this? How much of the full style declaration is required? I have very little scripting-fu.

Here is a section from a style file with random colour-fill choices:

<symbol>
  <lowervalue>Daw</lowervalue>
  <uppervalue>Daw</uppervalue>
  <label></label>
  <pointsymbol>hard:circle</pointsymbol>
  <pointsize>2</pointsize>
  <pointsizeunits>pixels</pointsizeunits>
  <rotationclassificationfieldname></rotationclassificationfieldname>
  <scaleclassificationfieldname></scaleclassificationfieldname>
  <symbolfieldname></symbolfieldname>
  <outlinecolor red="0" blue="0" green="0"/>
  <outlinestyle>SolidLine</outlinestyle>
  <outlinewidth>0.26</outlinewidth>
  <fillcolor red="215" blue="97" green="224"/>
  <fillpattern>SolidPattern</fillpattern>
  <texturepath></texturepath>
</symbol>

If I can get away with just

<symbol>
  <lowervalue>Daw</lowervalue>
  <uppervalue>Daw</uppervalue>
  <fillcolor red="215" blue="97" green="224"/>
</symbol>

and somehow merge the remainder, then I can probably do it in a spreadsheet.

Ben.


texts: R leaflet - Bring base layer labels to front

texts:
I want to use a base layer so that my filled polygons are in geographic context.
However, my polygons have simply gone over the top of the base layer and therefore only bits and pieces of place names are visible. The image above shows this.

1) Is it possible to have the place names and detail over the top of my choropleth?
2) If 1 is not possible, can I remove the names and detail to make it cleaner?

leaflet(MyMapb) %>%
setView(-1.88, 55, 6) %>%
addControl(rr, position = "topleft") %>%
addTiles(group = "OSM (default)") %>%
addProviderTiles("MapBox", options = providerTileOptions(
id = "mapbox.light",
accessToken = Sys.getenv('MAPBOX_ACCESS_TOKEN'))) %>%
#Polygon Layers
addPolygons(data = MyMapb,
          fillColor = ~pal1(MyMapb$eoi),
          fillOpacity = 1,
          color = "darkgrey",       ## color of borders between districts
          weight = 1.5,            ## width of borders
          popup = popup1,
          group = "EOIs") %>%
addPolygons(data = MyMapb,
          fillColor = ~pal2(MyMapb$Awards),
          fillOpacity = 1,
          color = "darkgrey",       ## color of borders between districts
          weight = 1.5,            ## width of borders
          popup = popup2,
          group = "Awards") %>%

#Legend Layers
addLegend(pal = pal1,values = ~eoi, group = "EOIs",
        position = "bottomleft") %>%
addLegend(pal = pal2,values = ~Awards, group = "Awards",
        position = "bottomright") %>%

#Layers Control
addLayersControl(overlayGroups = c("EOIs","Awards"),
               options = layersControlOptions(collapsed = FALSE)) %>%
hideGroup("Awards")
addLayersControl(baseGroups = "CartoDB.PositronNoLabels",
               overlayGroups = c("map labels",
                                 "MyMapb"))


texts:
How can I create a layout and import it into .jpg or .png format?
This code:

    from PyQt4.QtCore import *
    from PyQt4.QtGui import *
    from qgis.core import *
    from qgis.gui import *
    # create image
    layer = iface.activeLayer()
    img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)
    color = QColor(255, 255, 255)
    img.fill(color.rgb())
    # create painter
    p = QPainter()
    p.begin(img)
    p.setRenderHint(QPainter.Antialiasing)
    render = QgsMapRenderer()
    # set layer set
    lst = [layer.id()] # add ID of every layer
    render.setLayerSet(lst)
    # set extent
    rect = QgsRectangle(render.fullExtent())
    rect.scale(1.1)
    render.setExtent(rect)
    # set output size
    render.setOutputSize(img.size(), img.logicalDpiX())
    # do the rendering
    render.render(p)
    p.end()
  # save image
    img.save("C:/project_qgis/render.png","png")

give me:


texts: How to only delete vertices where the Polyline doesnt bend

texts:
I am working on pipes and I only need vertices where the lines bend. A lot of my pipe have vertices in the middle of the polyline even though there isn't a bend.

Is there any way to delete them with out me going through every polyline 1 by 1?

I am using ArcGIS Desktop 10.4.1 with an Advanced level license.


texts: How to write quotes in PG connection string in where clause when open Postgis Raster?

texts:
I have some rasters in PostGIS and I want to open just those which match to where clause. How can I do that? I have tried:

ds = gdal.Open("PG:host=<host> dbname=<dbname> table=<tablename> user=<user> password=<psswd> where='observ_date=\'2015-06-26\'::date'")

and some other combinations of quotes, but with no success. Everytime I got an error like this

ERROR 1: Error browsing database for PostGIS Raster properties : ERROR:  cannot cast type integer to date
LINE 1: ...r.alaska_mod13q1_ndvi where observ_date=2015-06-26::date gro...

It means that wrote quotes wrong. So how should I write quotes in where clause in PG connection string correctly?


texts:
I'm a beginner in ArcPy and I want to perform some tasks using ArcPy but I am not getting results.
I have two layers named survey and survey selection in my current MXD. I want to create map document for each layer in survey selection and this is my script:
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
titleElement = arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT")

rows = arcpy.SearchCursor("survey selection")
for row in rows:
    titleElement = row.VendorFarm
    query = '"FID"=' + str(row.FID)
    arcpy.SelectLayerByAttribute_management("survey selection", "NEW_SELECTION", query)
    df.zoomToSelectedFeatures()
    arcpy.RefreshActiveView()

del rows
del row


texts: How to Install GDAL Python Bindings Alongside osgeo4w for a separate Python version

texts:
I am currently attempting to install the Christopher Gholke distrubution of GDAL (and associated bindings) on a work computer as a dependency for a package. This question is related, but I need to run the version of the Python bindings because it is a Python 3 module, so running through the osgeo4w shell is not sufficient.

I am running Windows 7 Enterprise.

When trying to import any variation on:

import gdal #or
from osgeo import gdal # or
import osr # or
import ogr # or
#etc

It produced the error

 ImportError: DLL load failed: The specified module could not be found.

I successfully installed the Christopher Gholke distribution of GDAL on a separate machine with no errors. The machine that had a succesful install did not have osgeo4w installed, which leads me to believe it is the problem.

Is there any way to successfully import gdal from a Python 3.6 shell with osgeo4w without throwing the above error?


texts: GeoServer + PostGIS: multiple tables or SLD filters

texts:
I'm setting up a PostGIS and GeoServer for OpenLayers use and started wondering how to achieve the best performance in selecting and rendering the geometries.

Is it better to create multiple tables / materialized views to the database and then create layers and styles for each, and not use filters in the SLDs, or should all the data be dumped on a single table and then use SLD filters for applying the styles? And if SLD filters are the way to go, should I have multiple FeatureTypeStyles with one Rule each, or multiple Rules in single FeatureTypeStyle?