texts:
Can somebody help me please to georeference a map with Python?
Maybe I'm not setting right the row and column numbers?
This is my code:

import shutil
from osgeo import gdal, osr

orig_fn = 'G:/Egyetem/UBB Geomatica/Disertatie/EMSR365_AOI01_GRA_PRODUCT_r1_RTP01_v1.jpg'
output_fn = 'G:/Egyetem/UBB Geomatica/Disertatie/output2.tif'

# Create a copy of the original file and save it as the output filename:
shutil.copy(orig_fn, output_fn)

# Open the output file for writing for writing:
ds = gdal.Open(output_fn, gdal.GA_Update)

# Set spatial reference:
sr = osr.SpatialReference()
sr.ImportFromEPSG(32631) #My projection system

# Enter the GCPs
#   Format: [map x-coordinate(longitude)], [map y-coordinate (latitude)], [elevation],
#   [image column index(x)], [image row index (y)]

gcps = [gdal.GCP(280000, 4590000, 0, 417, 499),
        gdal.GCP(310000, 4590000, 0, 4354, 499),
        gdal.GCP(310000, 4560000, 0, 4354, 4437)]

# Apply the GCPs to the open output file:
ds.SetGCPs(gcps, sr.ExportToWkt())

# Close the output file in order to be able to work with it in other programs:
ds = None

This error occurs:

Here is my original map:
https://we.tl/t-dn6uIilxJq


texts: How can I call a .NET / C# tool (i.e. a .DLL file) from a model workflow or a Python script?

texts:
I've developed Add-in tools (with Windows dialog GUIs) for ArcMap which work fine, but I would also like to be able to call these directly (passing in parameters so I won't need to display the dialog) from either a Model Workflow or a Python Script.

Is this possible ?

As a workaround I've tried using inter-process coms (via a named pipe) with a worker task / thread and this works right up to the point where I attempt to store changes - when ArcMap itself goes legs-up (I'm guessing maybe a lock issue?).


texts: Seting 0 values in GeoTIFF to nodata using GDAL?

texts:
Is there a way to set all the 0 value pixels on a raster to nodata using GDAL?

I would like to use fillnodata to fill some holes, but holes are set to 0, so I need to delete them first


texts: How to translate tens of images from a directory as a batch with gdal_translate?

texts:
For example i have a dataset of images in *.cub format like

$path: /work/120614/mg_1164/
f01.img.cub_E.cub
f02.img.cub_E.cub
'
'
f0n.img.cub_E.cub

for converting a single image I use gdal_translate f01.img.cub_E.cub fo1.tif
for whole the datasets.HOW? batch file?


texts: Convert the_geom to text and then transform

texts:
I have this query that returns the_geom as text. However how could I also apply a transform to the same query. I can do this separately with the second example. Though I would like to do this to all the rows in my table, and in that example I copied the st_astext output into the transform query.

select  sitename,  st_astext(the_geom) from sites_tbl;

sitename   |                st_astext

Belfast A    | POINT(145592 530021)

ST_AsText(ST_Transform(ST_SetSRID('POINT(145592 530021)'::geometry,27700), 4326));


texts: Enforce PostGIS constraint that geometries do not overlap

texts:
I'm trying to figure out what is the most efficient way to create a constraint on a table that enforces the property that no two geometries within a table overlap with one another.

Is there an elegant way to enforce this with a table constraint, or would it be better to first just do a select query to check if the row to be inserted overlaps with any geometry already in the table.  This doesn't give me the same confidence in the integrity of the table, but at least should be fast as you are only comparing a single geometry to an indexed set of geometries.


texts:
I'm working in Django 1.9 and Postgres 9.4 with PostGIS extensions. I have a KML file, and I want to use it to update my geo fields. This is my model:

class PCT(models.Model):
    code = models.CharField(max_length=3, primary_key=True)
    boundary = models.GeometryField(null=True, blank=True)
    objects = models.GeoManager()

This is my current code:

    from django.contrib.gis.gdal import DataSource
    from django.contrib.gis.utils import LayerMapping
    ds = DataSource(options['filename'])
    layer_mapping = {
        'code': 'Name',
        'boundary': 'Unknown'
    }
    lm = LayerMapping(PCT, options['filename'],
                      layer_mapping, transform=False) # , unique='code')
    lm.save(strict=True)

It works OK, but if I run it a second time, it won't update any existing PCT fields - instead it creates new ones. This is a problem when I want to update the boundaries. (I can't just delete and re-create the fields, because there are foreign key dependencies.)

How can I use my KML file to update existing files?

@alex suggested trying this, but I'm getting an error:

    ds = DataSource(options['filename'])
    for layer in ds:
        for i in layer:
            code = str(i['Name'])
            geom = i.geom
            print geom

            geos_geom = geom.geos
            print geos_geom

            boundary = GEOSGeometry(geom)
            print boundary

            pct = PCT.objects.get(code=code)
            if pct:
                pct.boundary = geos_geom
                pct.save()

This throws an error at the boundary = GEOSGeometry(geom) line: TypeError: Improper geometry input type: <class 'django.contrib.gis.gdal.geometries.Polygon'>.

The output from the print statements up to that point is as follows:

POLYGON ((-1.407325 54.593612 0,-1.423288 54.603131 0,-1.431455 54.595297 0 ...
SRID=4326;POLYGON Z ((-1.4073249999999999 54.5936120000000003 0.0000000000000000

What am I doing wrong?

I tried simply passing a string to GEOSGeometry:

boundary = GEOSGeometry(str(i.geom))
pct = PCT.objects.get(code=code)
if pct:
  pct.boundary = boundary
  pct.save()

Perhaps this solves the issue? But now I get a different error on pct.save():

django.db.utils.DataError: Geometry has Z dimension but column does not

I can see by looking at the print statements that there is indeed a z column, but it's always zero, is it possible to get rid of this somehow?


texts: Splitting string ID code into various parts

texts:
I have a series of identification codes that I need to split out. The format of these codes is [region(letter)][district(number)] - [place(number)][subdistrict(letter)].

An example of some codes includes S22-201, TT100-12, and V6-1B. Often there is no subdistrict, and all points fall within the same larger district (so no As or Cs or whatever at the end of the string.

I can do parts of it, like splitting at the hyphen.

!Original_ID!.split('-')[0]

and then extracting the district

!Split_ID![1:3]

But it seems like two steps for this are unnecessary, and only works when I know the specific number of characters in the string, which isn't realistic for a large data set.

I'd like to be able to grab each piece at once:

letters on the left of the hyphen
numbers on the left of the hyphen
numbers on the right of the hyphen
letters (if any) on the right of the hyphen.

I'd need the numeric fields to be integers (or I guess possibly floats in some rare cases maybe).

I am still not doing something correctly. I may need to start smaller and brush up on my Python before I do this, I just assumed this would be a good place to start learning. Here's where I am at, in the Python window in ArcMap.

with arcpy.da.UpdateCursor("Wet_Sub",['Flag_ID','District','Split_ID']) as uCur:
for sRow in uCur:
    OrigID  = sRow[0].split('-')[0] # first element in the Original_ID
    charRng = range(len(OrigID)) # a range to iterate over
    Chars   = ''
    Numbers = ''
    for Idx in charRng:
        if OrigID[Idx].isnumeric():
            Numbers += OrigID[Idx]
        else:
            chars += OrigID[Idx]
    sRow[1] = float(Numbers)
    sRow[2] = Chars
    uCur.updateRow(sRow)

"Wet_Sub" and 'Flag_ID' are the names of the feature class and actual original field. I also tried to follow along with user2856's suggestion. It looks like I may need to be using both of those code blocks, one pasted into another, but I wasn't sure how to fit them together and what parts to change/remove (e.g. "etc... from code block above").


texts: Why does arcpy's Alter Field return "The tool is not licensed"?

texts:
According to the help file, the ALTER FIELD command is available at any license level.

ArcGIS for Desktop Basic: Yes

Why am I getting a license error when trying to rename a field alias on a sample file geodatabase featureclass?

> arcpy.Describe(u"e:/test.gdb/testFC").fields[2].name
> u'oldfieldname'

> arcpy.AlterField_management(u"e:/test.gdb/testFC", "oldfieldname", "newfieldname"  ,"newaliasname")

> ExecuteError: Failed to execute. Parameters are not valid.
> ERROR 000824: The tool is not licensed.

This is a brand new, standard, empty featureclass which I just created for the purposes of this test.


texts: How to Point Cluster With Dynamic Data In ArcGIS API

texts:
ArcGIS API is using this Demo to illustrate how Point clustering works with JSON data. In my case I am generating some point graphics from MySQL Database by a jQuery AJAX request like following code: As you can see I am storing the data in project and then access to project.Longitude, project.Latitude to create the points

request.done(function( mapInfo) {
    for (var i = 0; i < mapInfo.length; i++) {
            var project = mapInfo[i];
            var point = new Point(project.Longitude, project.Latitude);
            var pointSymbol = new SimpleMarkerSymbol();
            pointSymbol.setStyle(SimpleMarkerSymbol.STYLE_CIRCLE);
            pointSymbol.setSize(12);
            var pointInfoTemplate = new InfoTemplate();
            pointInfoTemplate.setTitle("Project Details");
            pointInfoTemplate.setContent('//JSON Data');
            var pointGraphic = new Graphic(point, pointSymbol).setInfoTemplate(pointInfoTemplate);
            points.push(pointGraphic);
}
for (i = 0; i < points.length; ++i) {
        map.graphics.add(points[i]);
}
});`

How can I bind this points with Cluster Points in ArcGIS API sample?


texts:
There are plenty of maps around, printed and otherwise, with numbered grids that do not directly equate to real-world coordinate systems.

Typically these grids will have 0,0 in the bottom-left corner - at least, in the most of cases that I have personally seen.

This seems to be at least an informal convention. My question is whether this convention is more established?

Are map-makers just gravitating towards something that they think most people will instinctively understand, or are there additional reasons for them to follow this approach?

Background:
All developer-created maps and most community-created maps in a particular computer game follow this convention. However a very small number of community-created maps have a Y-axis that goes in reverse, so 0,0 is in the top-left corner.

Not only do players often miss this, and thus misread coordinates, but even in-game tools misread these coordinates - both being one grid out on the Y coordinates.

For example, a point on the map, might be denoted as (0125, 0125) on the grid,  but read as (0125, 0135) by the in-game tools.

So as a matter of curiousity, I'm trying to find out which of the following cases is most accurate:

The map is 'wrong' - the convention is that 0,0 is bottom-left
The game tools are 'wrong' - there are additional, competing conventions that use top-left as 0,0, and said tools ought to cater for this eventuality
Neither are wrong - there are no firm conventions and though the many maps might tend towards 0,0=bottom-left, there are lots of other approaches commonly in use.

Note: A whole host of maps have A-Z+ as the x-axis and 0-n as the y-axis, with A1 being at the top-left, but I consider these as separate - these are intended to identify grids squares but not 6 and 8-figure coordinates.

I suspect the answer will be #3 but I will be interested to hear of any background/historical/supporting evidence.


texts: Does Identify tool of ArcGIS Desktop have Python script behind it?

texts:
Is there a Python script for "Identify" tool in ArcMap 10?
I need to create a customized script tool for that and add it to a toolbar.
The attribute table of my feature class has many fields. Out of that the custom identify tool should display the values of only two fields. i need to add that as a button in my tool bar.. when the user clicks any feature, the values of the two fields should be displayed.


texts: Starting animation after X seconds in OpenLayers 3?

texts:
Following the documentation,

A normal animation is working fine.

But an animation with a delay (parameter: start added to previous example) is not working.

Am I doing something wrong?


texts: Issue loading PostGIS out-db rasters (either in QGIS or MapServer)

texts:
I'm playing with a few different options of storing raster data via PostGIS to visualize (eventually in a web app).

The out-db or out-of-database storage seems the best option, but I can't seem to get it to work outside of PostGIS. Functions within PostGIS (like ST_Histogram) output fine results, but taking it outside fails.

First I tried MapServer since that's how I eventually want to handle it.

LAYER
    NAME test
    TYPE raster
    STATUS ON
    DATA "PG:host=website.com port=1234
          dbname='raster_dev' user='username' password='password'
          schema='public' table='rasters_outdb' where='\'rid\'=1'"
    PROJECTION
        "init=epsg:3310"
    END
    PROCESSING "NODATA=252"
    PROCESSING "RESAMPLE=NEAREST"
END

Results in this error:

msDrawMap(): Image handling error. Failed to draw layer named 'test'. drawGDAL(): Unable to access file. GDALDatasetRasterIO() failed: IReadBlock failed at X offset 0, Y offset 0

If I change table to rasters_indb (the in-db storage option), the map/png loads fine.

So I tried loading it up in QGIS. Again, in-db works fine. Out-db, I get a blank, black raster, with nan for min/max values -- however it does seem to occupy the right location and size. Progress?

I did a lot of googling and I seem to get mixed results on whether it works, whether GDAL can handle out-db or not, and most of it outdated. So I'm not sure where I stand here. Am I trying something that just isn't supported yet or am I doing something above wrong?

I'm on QGIS 2.12, MapServer 7, PostgreSQL 9.4, PostGIS 2.1


texts: Select by location based on selected features with ArcPy?

texts:
Is there a way with Python to do Select by Location with the Use selected feature option?

I'm trying to automate a process and am running into barriers.

I need to be able to select the 2 lines within the same feature class that intersect a point, but I need to do this 1 at a time so I can capture the field attributes for each line. The 2 lines will then be dissolved and eventually appended back into the original line feature. Basically, I'm trying to merge the lines and keep the attributes. The exact tool I need is the Merge tool in the Editor drop down, but apparently that isn't accessible with Python for some dumb reason.

There's also an Unsplit tool but of course it's only for Advanced, which I don't have.

What I'm trying to do:

From the Points feature class, get a list of the OBJECTID values and then use Select by attributes to select each point one at a time. But then I need to be able to do Select by Location on the Line using only the Selected point feature.

Then once the 2 lines are selected, I'll copy their field attributes and dissolve them and then append them back into the original line and add the values back.

I think I got it to work finally. The only thing I haven't figured out yet is how to delete the original lines that are being replaced when I append the newly dissolved lines. Because I think the only way I can use the Select by Location is when a layer file and not a feature class, so I'm not sure how to go about it yet.


texts: Remove optional tag for parameters user script processing toolbox QGIS

texts:
I have sucessfully written my first QGIS processing toolbox.
But have this bizzare optional tag on my parameters and i would like to remove these tags.

##Fichier_csv=file
##Dossier_de_sortie=folder
##nomodeler

import sys
import re
import os
import csv
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
import processing

reload(sys)
sys.setdefaultencoding('utf8')
final_list = []
f_csv = Fichier_csv
op = Dossier_de_sortie
with open(f_csv, "rU") as s:
    reader_tables = csv.DictReader(s,delimiter=';')
    liste_tables = list(reader_tables)
    l_tables = [i['table'] for i in liste_tables]
    l_tables = set(l_tables)
    for table in l_tables:
        l1 = [line for line in liste_tables if line['table'] == table]
        geometry = l1[1]['geometry']
        layer_mem = QgsVectorLayer(geometry+"?crs=epsg:2154",table,"memory")
        QgsVectorFileWriter.writeAsVectorFormat(layer_mem,os.path.join(op,table),"utf-8",None,"ESRI Shapefile",False)
        layer_path = os.path.join(op,table)
        shp = QgsVectorLayer(layer_path,table,'ogr')
        pr = shp.dataProvider()
        shp.startEditing()
        for field in l1:
            if field['type'] != 'Double':
                pr.addAttributes([QgsField(field['field'],getattr(QVariant,field['type']),field['typeName'],int(field['len']))])
            else:
                pr.addAttributes([QgsField(field['field'],getattr(QVariant,field['type']),field['typeName'],int(field['len']),int(field['precision']))])
        shp.updateFields()
        shp.commitChanges()
        del layer_mem,shp,pr


texts: How to modify spatialite table column names with QGIS?

texts:
I am trying to modify column names in a SpatiaLite table.
To do that, I run the following SQL statements in the DB Manager :

-- Rename existing table
ALTER TABLE compteur_routier RENAME TO compteur_routier_sav;

-- Create new table with modified column names
CREATE TABLE compteur_routier (
id INTEGER PRIMARY KEY,
geom POINT,
id_compt TEXT,
type_compt TEXT,
angle INTEGER,
id_tmp INTEGER
);

-- Insert data from renamed table
INSERT INTO compteur_routier(id, geom, id_compt, type_compt, angle, id_tmp)
SELECT id, geom, id_compteur, type_compteur, angle, id_tmp FROM compteur_routier_sav;

-- Drop renamed table
DROP TABLE compteur_routier_sav;

This works perfectly in the DB manager, column names are modified and datas are kept.
However, when I load the resulting table in QGIS and open the attribute table, I get old column names and datas replaced by column names.

Table before :

Table after :


texts: Difference between st_length() and st_m()

texts:
My question is about the differences between st_length() and st_m() in postgis.

The reason I ask this is because in a MultiLineString I can´t calculate st_m(), but I can st_length().


texts: OpenLayers removing/adding marker (airplane): update the orientation based on lat/long (next point in database)

texts:
The code below I have is working great. I'll try to explain it in a few lines (btw you can view it here http://codeshare.io/9Z1BE):

-there's a marker on the map. The lat long is returned as json from a database. There is only ever ONE marker on the map at any given time.

-Say marker is at position A. Stays there for 2 seconds, goes to position B. B to C. C to D. So on and so on. This A-B-C-D represents a flight path: ie, 'A' would be the origin, and 'D' would be the destination (say LAX to JFK).

-Now I would like to have the orientation adjust accordingly based on where the marker is going. Say you were flying from Sydney to Hong Kong: you're moving north. The marker/airplane should have a north orientation.

-After googling a bit it seems it's better to use a vector feature than a layer feature so you can use the rotation property. But how exactly do I implement this with the setTimeout function below?

    var map;
    var updateMarker = "";
    var markers;
    function init(){
      map = new OpenLayers.Map('map', {
        controls: [
          //controls here
    ],
         numZoomLevels: 10
    });
    markers = new OpenLayers.Layer.Markers( "Markers" );
    map.addLayers([layer1]);
}

function planeDelay() {
setInterval(function() {
$.ajax({
  url: "latlong.php",
  type: "POST",
  success: function(msg){
    var json_coords = eval('('+msg+')');
    if(updateMarker != "")
    {
            markers.removeMarker(updateMarker);
    }
    updateMarker = new OpenLayers.Marker(new OpenLayers.LonLat(json_coords.longitude, json_coords.latitude));
    markers.addMarker(updateMarker);
  }
});
}, 2000);
}


texts: Creating point at highest value in raster using ArcGIS Desktop

texts:
I am using ArcGIS Desktop 10.1.

I am trying to find the location of the highest value pixel in my Direct Solar Radiation Raster. I can see the extent of the data in my raster but don't know where the highest value pertains to.

Is there a way to simply select pixel with highest value or create a point at highest value?


texts: Convert Point to Arrow (SymbolizedPoint)

texts:
So I managed to get a point to show up at the center of my line, I'm hoping to convert the point into an arrow pointing in the direction of my line

public static final SimpleFeatureType POINT = createPointType();
public static final SimpleFeatureType LINE = createSolidLine();

private static SimpleFeature directedPoint(){
    Coordinate start = CoordConverterUtils.getStart;
    Coordinate end = CoordConverterUtils.getEnd;

    Coordinate mid = (new LineSegment(start,end)).midpoint());

    GeometryFactory gf = JTSFactoryFinder.getGeometryFactory();
    SimpleFeatureBuilder fb = new SimpleFeatureBuilder(DIRECTED_POINT);
    Point p = gf.createPoint(mid);
    fb.add(p);
    return fb.buildFeature(null);
}

private static SimpleFeatureType createPointType(){
    try{
        SimpleFeatureTypeBuilder builder = new SimpleFeatureTypeBuilder();
        builder.setName("DirectedPoint");
        builder.setCRS(DefaultGeographicCRS.WGS84);
        builder.add("the_geom", Point.class);
        return builder.buildFeatureType();
    catch(Exception e){
        //warning
        return null;
    }
}

Is it possible to convert the Point into a SymbolizedPoint that's an arrow pointing in the direction of LINE - and return it as a SimpleFeature ?  Every time I've tried to make a symbolized point I can't seem to translate the midpoint into its geometry.


texts: "Warning, inconsistent extent!" error after using project raster for MODIS data

texts:
I'm using MODIS satellite data that I've converted into GeoTIFFs (using extract subdataset). ArcMap automatically sets the geographic coordinate system to GCS_Unknown_datum_based_upon_the_Clarke_1866_ellipsoid. After some digging around on the MODIS site, I think, but am not entirely sure, that the starting coordinate system is supposed to be GCS WGS 1984. My goal is to get these into Lambert Azimuthal Equal Area projection and ultimately clip them to show the arctic region. I've tried a few different methods, but when the full image is projected, I always get the inconsistent extent error. For the record, I'm using ArcMap 10.3.1.

What I've tried:
 - Project raster as is
 - Define projection (set to GCS WGS 1984) and then project raster (Define projection gives me a warning that it already has a projection since ArcMap has technically given the images a coordinate system)
 - Project raster into GCS WGS 1984 and then project raster into LAEA
 - Going into ArcCatalog to manually change the spatial reference (no way to edit it)

For projecting into LAEA, I've done it both specifying an output cell size and using the default output cell size.

The images always look correct once they're clipped (all of the land area looks to be there and in the general shape it should be) and I don't get the error when I add the clipped images to the data frame, but I'm assuming there's some small distortion or issue with these final images.

I'm processing upwards of 800 images so I'm trying to find a solution that can be done either through batch processing or model builder.


texts: Counting pixels in masked raster in Google Earth Engine results in error 'Remote request too large'

texts:
Following on from a previous question (Points to raster (image) in Google Earth Engine) I have created an Image of observations for a particular species, masked such that we are only considering pixels with a value of one (i.e. the species was observed in that pixel).

The extent of the image is necessarily large, because the observations were spread out, which means that the whole image contains many, many pixels but most are zeroes and should therefore be ignored in calculations.

All I want to do is count the number of pixels with a value of one. I tried the following code:

// Mask
var ebdMasked = ebdCountImg.mask(ebdCountImg.gt(0));

// Count pixels
var count = ebdMasked.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  maxPixels: 1e13
  });
print(count);

Which produced the error:

Dictionary (Error)
Remote request too large (3355443200 > 83886080) for output: [<Object>].

I actually know in this case that the answer is 256, so in no way should this be too large for output. I assumed something was wrong with the mask meaning that all pixels were being considered, so I exported the masked Image to take a look. It seemed fine and, lo and behold, when I read the exported ebdMasked back in, and try again, it works as expected:

// Try again with exported layer
var ebdMasked2 = ee.Image("users/rasenior/ebdImg");

// Mask
var ebdMasked2 = ebdMasked2.mask(ebdMasked2.gt(0));

// Count pixels
var count2 = ebdMasked2.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi,
  maxPixels: 1e13
  });

print(count2);

Produces:

Object (1 property)
   first: 256

I'm sure I'm missing something very simple!

Full code is here:
https://code.earthengine.google.com/6608ad2279f249324ed730a56605a025


texts: Free Tide prediction service - World (a freeware program would work as an alternative)

texts:
Is there a free service where I could input a GPS location and it would return me the current tide situation at the position?
As an alternative - some library / program that does the exact thing one could integrate? For my app I am currently using the public-domain, accurate NOAA predictions, but these are limited to the US. :(
Several apps provide that data and I am curious what their source is.


texts: geostatistics across different projections?

texts:
What are some best practices for calculating local geostatistics across the boundaries of different projections?

I am thinking of using the State Plane Coordinate System to represent a national-level dataset. However I will be conducting analysis at the local metropolitan area level. I'd like to hear how people have handled analysis of zones that cover multiple projections (like the Kings County state plane in the Seattle region).

EDIT: the following was previously a comment:
We are applying convolutions on rasterized land use data to learn significant land use features. It seems like area will be the most important feature to preserve, but distance is a close second. My hope was that by dividing in the country into 120+ projections, that SPCS would essentially make such distortions insignificant. But that of course introduces inconsistencies when convolving over the boundary between two different projections


texts: Why the indexes for the attributes are not being used when GeoServer asks PostGIS?

texts:
I am working on a WMS to publish data with GeoServer (2.11.0). Tables come from PostGIS (PostgreSQL 9.6). I have indexed geometries on the tables and i have also indexed attributes which are being used in the SLD.

The SLD is quite simple, this is an example of one rule:

<Rule>
      <Name>Secondary_tunnel_stor-1</Name>
      <ogc:Filter>
        <ogc:And>
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>highway</ogc:PropertyName>
            <ogc:Literal>secondary</ogc:Literal>
          </ogc:PropertyIsEqualTo>
          <ogc:PropertyIsEqualTo>
            <ogc:PropertyName>tunnel</ogc:PropertyName>
            <ogc:Literal>yes</ogc:Literal>
          </ogc:PropertyIsEqualTo>
        </ogc:And>
      </ogc:Filter>
      <MinScaleDenominator>68001</MinScaleDenominator>
      <MaxScaleDenominator>272000</MaxScaleDenominator>
      <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#666666</CssParameter>
          <CssParameter name="stroke-width">3</CssParameter>
          <CssParameter name="stroke-dasharray">8 4</CssParameter>
        </Stroke>
      </LineSymbolizer>
      <LineSymbolizer>
        <Stroke>
          <CssParameter name="stroke">#faeae6</CssParameter>
          <CssParameter name="stroke-width">1.5</CssParameter>
        </Stroke>
      </LineSymbolizer>
    </Rule>

But after panning through the WMS I can see on PostGIS that the only index being used is the geometry index:

Why the indexes for the attributes are not being used?

The query sent from GeoServer to PostGIS is like:

SELECT
"id_pkey",
"natural",
"landuse",
"waterway",
encode(
    ST_AsBinary ( CASE WHEN ST_HasArc ( "the_geom" ) THEN "the_geom" ELSE ST_Simplify ( ST_Force2D ( "the_geom" ), 7803.3542335301645, TRUE ) END ),
    'base64'
) AS "the_geom"
FROM
"myschema"."mytable"
WHERE
(
    "the_geom" && ST_GeomFromText ( 'POLYGON ((11888876.271085668 4251249.343989519, 11888876.271085668 11795995.306196207, 17000073.294047713 11795995.306196207, 17000073.294047713 4251249.343989519, 11888876.271085668 4251249.343989519))', 3857 )
    AND ((
            "natural" = 'wood'
            AND "natural" IS NOT NULL
            )
        OR ( "landuse" = 'forest' AND "landuse" IS NOT NULL )
        OR ( "natural" = 'glacier' AND "natural" IS NOT NULL )
        OR ( "waterway" = 'glacier' AND "waterway" IS NOT NULL )
        OR ( "natural" = 'sand' AND "natural" IS NOT NULL )
        OR ( "natural" = 'grassland' AND "natural" IS NOT NULL )
        OR ( "landuse" = 'residential' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'commercial' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'retail' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'industrial' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'farmland' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'farmyard' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'orchard' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'vineyard' AND "landuse" IS NOT NULL )
        OR ( "landuse" = 'grass' AND "landuse" IS NOT NULL )
        OR ( "natural" = 'water' AND "natural" IS NOT NULL )
        OR ( "waterway" = 'riverbank' AND "waterway" IS NOT NULL )
    OR ( "landuse" = 'reservoir' AND "landuse" IS NOT NULL )
OR ( "waterway" = 'canal' AND "waterway" IS NOT NULL )))


texts: Viewshed analysis with mutliple observers in ArcGIS Modelbuilder

texts:
I want to perform a viewshed analysis using multiple obeservation points. Each observation point has an observer height of 1.7 m and a buffer visibility distance of 5 km. In order to generate a buffer of 5 km I populated the my testpoint dataset with an attribute table RADIUS2 with 5000 meter, corresponding to the buffer viewshed limit. I used ArcGIS 10 to build a model builder that fairly does what I want. See Figure 1.

I tested the model by using three points: Point 1&2 in proximity of each other and therefore generating visibility value of 2. Then I generated a third point at distance bigger 5 km from the other two points in order to check if visibility buffer of 5 km in my model works. As you can see from figure 2 is does not, because point 3 should have value of 1 (blue), instead it is still influenced by the viewshed of point 1&2. Can somebody help me figure out what is wrong in my modelbuilder?


texts: Heatmap.js is not displaying on OpenLayers

texts:
I tried to follow the example on the heatmap.js but for some reason it is not displaying the heat map layer, and the layer control is not working. Maybe I am missing something stupid. Can someone point me out what goes wrong?

$(function() {
 var testData={
max: 46,
data: [{lat: 33.5363, lng:-117.044, count: 1},{lat: 33.5608, lng:-117.24, count: 1},{lat: 38, lng:-97, count: 1},{lat: 38.9358, lng:-77.1621, count: 1}]
};
var transformedTestData = { max: testData.max , data: [] },
data = testData.data,
datalen = data.length,
nudata = [];

// in order to use the OpenLayers Heatmap Layer we have to transform our data into
// { max: , data: [{lonlat: , count: },...]}
while(datalen--){
nudata.push({
lonlat: new OpenLayers.LonLat(data[datalen].lon, data[datalen].lat),
count: data[datalen].count
});
}
transformedTestData.data = nudata;
var map = new OpenLayers.Map('map_canvas',
                {
                    allOverlays: true,
                    projection: new OpenLayers.Projection("EPSG:900913"),
                    displayProjection: new OpenLayers.Projection("EPSG:900913"),
                    units: "m",
                    maxResolution: 156543.0339,
                    maxExtent: new OpenLayers.Bounds(-20037508, -20037508, 20037508, 20037508),
                    controls: [
                        new OpenLayers.Control.Navigation(),
                        new OpenLayers.Control.PanZoomBar(),
                        new OpenLayers.Control.LayerSwitcher(),
                        new OpenLayers.Control.ScaleLine(),
                        new OpenLayers.Control.MousePosition()
                    ]
                });
 var layer = new OpenLayers.Layer.OSM();
// create our heatmap layer
var heatmap = new OpenLayers.Layer.Heatmap("Heatmap Layer", map, testData,
                {visible: true, radius: 15},
                {isBaseLayer: false, opacity: 0.3, projection: new OpenLayers.Projection("EPSG:4326")});
//map.addLayer(layer);
map.addLayers([layer,heatmap]);

map.zoomToMaxExtent();
heatmap.setDataSet(transformedTestData);
// Handler for .ready() called.
});


texts: How to create schematic attribute as associated field using ArcObjects?

texts:
I'm trying to automate the complete creation of a schematic template with ArcObjects SDK for Java (10.2).

I've managed to create the GDB, the schematic dataset and the schematic diagram class (template).

Then, I created a schematic element class associated with a feature class taken from the same GDB, which works.

I would like to share a field "NAME" between the feature class and the schematic element class, in order to label the nodes of a schematic diagram.

So, I tried that:

...
 ISchematicElementClass nodeSchematicElementClass = schematicDataset.createSchematicElementClass("NODE", esriSchematicElementType.esriSchematicNodeType);
 nodeSchematicElementClass.setSchematicDataSourceByRef(schematicDataSource);
 nodeSchematicElementClass.alterAssociatedObjectClass(schematicDataSource, nodeFeatureClass.getObjectClassID());
 nodeSchematicElementClass.setGeometryType(esriGeometryType.esriGeometryPoint);
 nodeSchematicElementClass.setExternalQueryEvaluationMode(esriSchematicExternalQueryEvaluationMode.esriSchematicQueryBuildEvaluation);
 nodeSchematicElementClass.setQueryString("SELECT * FROM NODE");
 schematicDiagramClass.associateSchematicElementClass(nodeSchematicElementClass);

 UID associatedFieldUID = new UID();
 associatedFieldUID.setValue("{7DE3A19D-32D0-41CD-B896-37CA3AFBD88A}"); // taken from the file ...\Desktop10.2\bin\configuration\CLSID\esri.clsid.ecfg
 ISchematicAttribute schematicAttribute = nodeSchematicElementClass.createSchematicAttribute("NAME", associatedFieldUID);

then I get the following ArcObjects error in the last line:

  AutomationException: 0x80040a28 - Additional field can't be added : NAME in 'Esri Schematics'

I get the same error with any other field name.

Does anybody have a solution? Am I doing something wrong with my attribute field association?


texts: Batch Exporting Layouts to PDF in ArcGIS Pro using ArcPy?

texts:
Is there a way to batch export layouts in an ArcGIS Pro document to a specific format?

For example exporting all layouts to PDF in a file directory.

This was possible with multiple MXDs using python with Desktop, however I haven't found a way to do this with Pro.


texts: add features to a featurelayer and save changes in the server

texts:
I am developping a web application with arcgis javascript api. I am using arcgis server 10.2.
I declared a featureLayer, the URL of this featureLayer is a feature service on my server.
Using a geoprocessing service, I am able to create features from the results of the service. I am able to add the features to the featureLayer I declared and show them in the map, but the problem is the changes made are not saved on the server.
this is the code sample:

 function displayResults(results, messages) {
    var flayer = new esri.layers.FeatureLayer("http://hamri-pc:6080/arcgis/rest/services/edt/FeatureServer/0", {
      mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
      outFields: ['*']
    });
    var features = results[0].value.features;
    var theExtent = null;
    for (var f = 0, fl = features.length; f < fl; f++) {
        var feature = features[f];
        var attribs = feature.attributes;
            feature.setSymbol(normalpictureMarkerSymbol);
            feature.setInfoTemplate(infoTemplate);
                            flayer.add(feature);
                      flayer.applyEdits(feature,null,null);
                            map.addLayers([flayer]);
        }
    }

Does anyone have an idea of what should I do to save the features in the server?


texts: Is there a GIS layer file for the USGS Quad Boundaries in NAD27 coordinate system?

texts: USGS topographic maps currently use the NAD83 coordinate system and all the layer files for the quad boundaries I can find do too. I'm wondering if there exists a nationwide GIS layer for the boundaries of the old NAD27 coordinate system that align with the neatlines for topographic maps pre-1983.


texts: What unit is the topographic openness (SAGA GIS)?

texts:
The SAGA GIS tool "Topographic Openness" calculates the positive and negative topographic openness and is based on Yokoyama et al. (2002) who do not specifically say the unit of the final result. They refer to intermediate angles as > or < 90° (which are then added together and divided by the number of angles), but the maps I get in SAGA have values between roughly 0.6 and 1.7 which would equal 35 to 97 degrees if the unit were radians. So which is it? Or is it unitless?

Maybe the source code helps (I don't understand much of it, unfortunately...): https://sourceforge.net/p/saga-gis/code/ci/master/tree/saga-gis/src/tools/terrain_analysis/ta_lighting/topographic_openness.cpp

Literature

Yokoyama, R. / Shirasawa, M. / Pike, R.J. (2002): Visualizing topography by openness: A new application of image processing to digital elevation models. Photogrammetric Engineering and Remote Sensing, Vol.68, pp.251-266. (Link)


texts: Leaflet making features gone when zoom out

texts:
I am wondering about an option to make all or at least majority of placemarks (depends on their internal zoom settings) invisible when zooming out.

It can be something alternative to the leaflet markerclusters:

https://github.com/Leaflet/Leaflet.markercluster

where the placemarks are clustered as we are zooming out.

My goal is to make them completely dissapeared when my zoom level drops. If you don't know what I mean you can open some tile maps i.e. OpenstreetMap or Google. Having a zoom 1 level you have got a blank map, whereas your zoom level increase, a new features are coming up.

However i don't know, why this attidude doesn't apply to the markers or another features created by an user. It refers both to Google, Google Earth, Leaflet and another interactive mapping tools.

This situation applies to the polygons too, although due to their square features, they tend to be smaller when zoom out.

Is there some way to make it gone, at least in the Leaflet map software?
rise an opposite situation as below:

Features disappear from Leaflet after zooming in


texts:
I have a map. I have external data that I can join to the map by name of the feature. The data contains several variables, say abc, def, xyz1, xyz2 and xyz3.
I can display a variable by selecting it in the symbology of the map layer. I created a print layout with a legend which works for all variables xyz*.
How can I (batch) generate separate maps (in png) for all variables xyz*?


texts: SpatiaLite not complete data on pt_addresses?

texts:
I have a problem.
Imported OSM in SQLite using spatialite_osm_map
Samples query

$select id,name,min(st_distance(geometry,st_point(22.908082,45.763242),0)) from roads limit 1;
564|Strada Muntenia|26.944248068482

Is ok, is my street

select id,name,min(st_distance(geometry,st_point(22.908082,45.763242),0)) from roads limit 1;
564|Strada Muntenia|26.944248068482

Same is ok, is my street again

select id,sub_type,name,min(st_distance(geometry,st_point(22.908082,45.763242),0)) from pg_building limit 1;
33042016|yes|Bloc D2|6.50190449190996

Is ok, my address is "Strada Muntenia Bloc D2"

I want to get the rest of data for my address, my city is Hunedoara

select id,country,city,postcode,street,housename,housenumber from pt_addresses where city like 'Hun%';
NULL

?!?!, no city Hunedoara in pt_addreses ??!?!

select id,country,city,postcode,street,housename,housenumber from pt_addresses where housename like '%D2%' or housenumber like '2A%' limit 1;
701843031|RO|Deva|330051|Aleea Militarilor||2a

Something, but wrong city and street.
I must understand Hunedoara is missing from pt_addresses ?


texts: How to dynamically change text of checkable icon?

texts:
In my QGIS plugin, I have a checkable QIcon which changes image depending on its state.

Is it also possible to change its tooltip text (i.e. if icon is checked, the text would read "myPlugin ON"; unchecked would read "myPlugin OFF")?

At the moment, it only reads "myPlugin". This is the code I use:

def add_action(
    self,
    toolbar,
    icon_path,
    text,
    callback,
    checkable=False,
    enabled_flag=True,
    add_to_menu=True,
    add_to_toolbar=True,
    status_tip=None,
    whats_this=None,
    parent=None):

    icon = QIcon(icon_path)
    action = QAction(icon, text, parent)
    action.triggered.connect(callback)
    action.setEnabled(enabled_flag)
    action.setCheckable(checkable)

    if status_tip is not None:
        action.setStatusTip(status_tip)

    if whats_this is not None:
        action.setWhatsThis(whats_this)

    if add_to_toolbar:
        self.toolbar.addAction(action)

    if add_to_menu:
        self.iface.addPluginToMenu(
            self.menu,
            action)

    self.actions.append(action)
    return action

def initGui(self):
    icon = QIcon()
    icon.addFile( ':/plugins/myPlugin/icons/Off.png', state=QIcon.Off )
    icon.addFile( ':/plugins/myPlugin/icons/On.png', state=QIcon.On )
    self.pluginButton = self.add_action(
        self.toolbar,
        icon,
        text=self.tr(u'myPlugin'),
        callback=self.myFunc,
        checkable=True,
        enabled_flag=False,
        parent=self.iface.mainWindow())

I am using QGIS 2.18.2 for Win7 64-bit.


texts: Copy/Paste attribute table qgis without the wkt_geom column?

texts: Is there a way to copy the attribute table on qgis and paste it on excel without the wkt_geom column?


texts: Looking for examples of bayesian modelling and smoothing in a GIS

texts:
I am looking to assemble clear examples of Bayesian modelling / smoothing in a GIS. I have yet to see GIS being used to conduct Bayesian style analytical models and would like to see what is currently out there. I am interested in both the application of Bayesian modelling and its current state of inclusion in traditional and/or non-proprietary (opensource) software.

At the moment I am very interested in Bayesian modelling for transportation planning, specifically transit ridership and modal split/share estimations.


texts: Conflicting sizes for dimension 'time' with rioxarray

texts:
I've been using rioxarray to convert netcdf files to geotiffs, and up till now, it's worked great.
However, when I open a ERA5 reanalysis netcdf dataset, rioxarray appears to find that the data size and the 'time' dimension are different (72144 vs 32768):
>> data = rio.open_rasterio(filename)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/me/opt/anaconda3/lib/python3.8/site-packages/rioxarray/_io.py", line 756, in open_rasterio
    result = DataArray(
  File "/Users/me/opt/anaconda3/lib/python3.8/site-packages/xarray/core/dataarray.py", line 344, in __init__
coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
  File "/Users/me/opt/anaconda3/lib/python3.8/site-packages/xarray/core/dataarray.py", line 152, in _infer_coords_and_dims
    raise ValueError(
ValueError: conflicting sizes for dimension 'time': length 32768 on the data but length 72144 on coordinate 'time'

When I open the same file in xarray, however, I find that the 'time' dimension and the data length are equal (72144):
>> xdata = xr.open_dataset(filename)
<xarray.Dataset>
Dimensions:    (latitude: 53, longitude: 89, time: 72144)
Coordinates:
  * longitude  (longitude) float32 20.0 20.25 20.5 20.75 ... 41.5 41.75 42.0
  * latitude   (latitude) float32 55.0 54.75 54.5 54.25 ... 42.5 42.25 42.0
  * time       (time) datetime64[ns] 2010-01-01 ... 2020-12-31T23:00:00
Data variables:
    t2m        (time, latitude, longitude) float32 ...
>> xdata.t2m.shape
(72144, 53, 89)

Why might rioxarray think the 'time' dimension in my netcdf is too short for the data, when xarray thinks the dimension and the data are the same size?


texts: Seeking PostGIS equivalent to GeoPandas union overlay?

texts:
What is the PostGIS / PostgreSQL equivalent of the GeoPandas Overlay Union operation?

res_union = geopandas.overlay(df1, df2, how='union')

Say I have two tables, 1 green, 1 red with each with two rows with geometry (polygons):

Then, the objective is to get:

which is a table with 7 rows/ features. See the GeoPandas docs for more info.

I've tried spatial joins and using the ST_Union and ST_Intersection to no avail. I managed to reproduce the equivalent of

res_intersection = geopandas.overlay(df1, df2, how='intersection')

Which is the equivalent of:

using the following SQL command:

SELECT
    table1.letter,
    table2.number,
    ST_Intersection(table1.geom,table2.geom)
FROM
    table1,
    table2
WHERE
    ST_Intersects(table1.geom,table2.geom)

One solution is to create a table for the intersections and one for the symmetrical differences and than outer join. However I don't know how to create the Symmetrical Difference:

SELECT
    test.hybas_nld.pfaf_id,
    test.gadm_nld.gid_1,
    ST_SymDifference(test.hybas_nld.geom,test.gadm_nld.geom)
FROM
    test.hybas_nld,
    test.gadm_nld
WHERE
    ST_????(test.hybas_nld.geom,test.gadm_nld.geom)

with the result of:

Possible duplicate with the exception that I would like to keep the polygons where only geometry of table 1 or 2 exists (symmetrical difference).

What would be the simplest and fastest way to get to the final result (Union)?


texts: How to create sub-menus in QGIS 3 plugin

texts:
I'm trying to create a sub-menu in a QGIS plugin like in the image. I found this Add extra item to drop down menu plugin. "Layers menu from project" plugin but I'm using QGIS 3 and PyQT5 and I don't know how to adapt it.
This is my code:

def initGui(self):
  icon = os.path.join(os.path.join(cmd_folder, 'logo.png'))
  self.action = QAction(QIcon(icon), 'Select Feature', self.iface.mainWindow())
  self.action.triggered.connect(self.run)
  self.iface.addPluginToMenu('&Select Feature', self.action)
  self.iface.addToolBarIcon(self.action)
  self.first_start = True

def unload(self):
  self.iface.removeToolBarIcon(self.action)
  self.iface.removePluginMenu('&Select Feature', self.action)
  del self.action

EDIT: I solved in this way, but I'm wondering if there's a better way

def initGui(self):
  icon = os.path.join(os.path.join(cmd_folder, 'logo.png'))
  self.action = QAction(QIcon(icon), 'Select Feature', self.iface.mainWindow())
  self.action2 = QAction(QIcon(icon), 'Select Feature2', self.iface.mainWindow())
  self.action.triggered.connect(self.run)
  self.action.triggered.connect(self.run2)
  self.iface.addPluginToMenu('&Select Feature', self.action)
  self.iface.addPluginToMenu('&Select Feature', self.action2)
  self.iface.addToolBarIcon(self.action)
  self.first_start = True

def unload(self):
  self.iface.removeToolBarIcon(self.action)
  self.iface.removePluginMenu('&Select Feature', self.action)
  self.iface.removePluginMenu('&Select Feature2', self.action2)
  del self.action
  del self.action2


texts:
I have a layer of points that are spread out in different areas.  I created a polygon layer for each area, and tried to count the number of points in each polygon using Vector > Analysis Tools > Points in Polygon.  QGIS 2.0.1 is crashing every time.  I have already checked and both layers are using the same CRS (WGS 84).
Is this a bug, or am I doing something wrong?


texts: How to assign unique field name for "Extract Multi Values To Point" operation using arcpy in loop?

texts:
I have 1000 of raster file for which I am trying to extract multi values to point in a loop so that I have a series of values for each coordinate.

Here I have one raster for each week, and I want to name the field accordingly in the point shapefile's attribute table while extracting the values.

How to supply the variable name in the command line shown below

arcpy.gp.ExtractMultiValuesToPoints_sa(path/to/pointshapefile,"raster_2012_1.tif 2012_1","NONE")

Above code line works fine when "raster_2012_1.tif" is open in window and "a specific name "2012_1" is to be assigned. However, in case of varying raster files in loop and varying field name variables are required to be assigned.

I tried to assign variables as follows,

arcpy.gp.ExtractMultiValuesToPoints_sa(path/to/pointshapefile, path/to/raster.tif field_name,"NONE")

Here I am updating field name in each loop, but it shows following error

Parsing error SyntaxError: invalid syntax (line 5)

where (line 5) is the shown code line.

I guess this problem must be occurring because of the number of arguments in the command. How to overcome this problem?

Edit: Sample code

import arcpy
shapefile = 'C:\Users\Administrator\Desktop\DSBShape\Point.shp'
tif = 'Raster.tif'
header = 'Header'
arcpy.gp.ExtractMultiValuesToPoints_sa(shapefile, tif header, "NONE")


texts: Is it possible to convert the USNG grid to a perfect square 10mx10m fishnet at CONUS scale?

texts:
I realize this is a ridiculous question, so please forgive in advance.
Would it be possible to generate a 10m x 10m grid at CONUS scale, without any non-square polygons at the Grid Zone Junctions (GZJs)? The desired result would have all GZJs removed from the USNG, and could be served to the web as a vector tile layer. Instead of slivers or polygons at the GZJs, the desired result would look like a fishnet containing perfect squares throughout.
I don't think just generating a fishnet would accomplish this result, because generating 10m x 10m fishnet vector polygons would be significantly computationally burdensome. My question then, is do you think this gird be produced at CONUS-scale without creating some kind of monster vector shapefile? I can't use rasters to solve this problem, as the grid cells are designed to be attributed with multiple variables.
Again, I understand this doesn't really make sense, but I was wondering if anyone had any bright ideas as to how this might be accomplished.


texts: Using 1:n table relationships in rule-based symbology

texts:
I would like to use the value of a field/attribute in a parent table as a filter in an expression used to display features in the child table.

Specifically, I have a 1:n relationship in my project:

Sectors (parent table with polygon geometry, a key "FID" and a field/attribute "TYPE")

Trees (child table with point geometry, and a secondary key "SECTORS_ID")

I would like to use this relationship in the rule-based symbology rendered so that, for example:

Trees whose parent Sectors have "TYPE" = 'Project' are displayed as a blue circle

Trees whose parent Sectors have "TYPE" = 'Parking' are displayed as a red circle

My question:

How do I build my expression in the rule-based renderer so that I can differentiate between Trees related to a Sector of a specific type?

The condition is not spatial (refFunction plugin is not used), and it should only be based on the value of the parent attribute "TYPE".

I'd rather avoid the following, but: how would I create a virtual field "SECTOR_TYPE" in my child table that would be filled with the value of the "TYPE" field located in the corresponding parent feature?

I am using QGIS 3.10.1


texts:
I know that you can do an https request in python, passing in a cert or bundle to requests.post(), with the verify argument.

Is there any way to do something similar in ArcPy?

I'd like to fully automate some code, but can't get a permanent certificate yet.

Failing that (since I'm on windows), is there a way to have it check the system trusted certificates?


texts: flow direction and flow length in arcmap

texts:
I am trying to replicate a flow length calculation in arcmap 10.2. I have an elevation raster with a pixel range from 0 to 3038.96m (32 bit). When I run the flow direction calculation, the pixel range drops to 0-255 (8 bit). I must be losing information that will affect my flow length calculation (I have a flow length calculation that has been completed that I'm trying to match)--how do I prevent this? I've tried searching in environment settings but I haven't found how to change my pixel range.

I'm fairly new to GIS.


texts: Can I use a 4 Band composite band tif file to make blend out borders around an image?

texts:
In ArcGIS 10.6 I have a 3 Band tif image which has large areas which I do not want to display. I have created polygons of these areas and converted them to a raster tif in order to create a 4th band to use as an alpha channel to make areas transparent. I used the Composite band tool to create the 4 band raster but when I activate the alpha channel using this band the whole image becomes transparent and not just the defined areas.

The 4th Band consists of just 2 values 0 and 1. Is ArcGIS assuming that the extent of the 4th band is the transparent area and not just the pixels with a value of 1? How to I set the raster to only make pixel value 1 the transparent pixel.

I am aware than I can set a nodata color in ArcGIS but I need to try to do this using an alpha channel for a client.


texts: How to load stored composer from project file in standalone QGIS application?

texts:
I would like to load a .qgs file and get a composer from it, but I can't figure out how.  Is there a way, or do I have to read the xml myself and individually set the appropriate values for the composer?

My question is somewhat related to this question, only I don't want to specify composer attributes myself, just load them from a .qgs file.


texts: Extracting paddy transplant dates (like SoS in TIMESAT)

texts:
I am trying to get transplant dates of paddy in Northern India using this code.

The idea is to find 15% of the maximum of an SG fitted curve of MODIS EVI, and find the Day of Year which corresponds to when the ascending limb had that fraction, like the Start of Season definition in TIMESAT while setting the method to seasonal amplitude.

The approach currently is to subtract this fraction from the images, and find which values lie closest to it by running quality mosaic, thus also retaining the Day of Year.

function percentCompEVI(collection,studyArea,percent) {
  var pc15 = amplitude.select('EVI').multiply(percent);
  var SoS = qmMin.select('EVI').add(pc15.select('EVI')).addBands(qmMin.select('DayOfYear'))
  print ('pc15 name test',pc15);
  print ('SoS value', SoS)
  // Compute distance of every pixel from the 10th percentile in that location
  var withDist = collection.map(function(image) {
    return image.addBands(image.select('EVI')
    .subtract(SoS.select('EVI')).abs().multiply(-1)
    .clip(FC).rename('quality'))
  })
  return withDist.qualityMosaic('quality').abs()
  //Retain the corresponding DayOfYear value
  return withDist.select('DayOfYear').abs()
}
var percent=0.15; //range: 0-1
var percentCompImage = percentCompEVI(data_,FC,percent);

The goal is to do what TIMESAT does. i.e., add the fraction of

[Value @ t+1  -  Value @ t] / [15% fraction - Value @ t]

to the image index that corresponds to the image at the time t.
I am assuming this approach would return Day of year values other than the distinct 8-day intervals; but I am finding it hard to apply on GEE.


texts:
How can I get a 33% equal split for 3 map panes - I keep getting the middle map being smaller than the left and right ones for some reason.  Theres no header or footer - just 3 maps.  If anyone can help that would be great.

Code below:

#mainWindow {
    width: 100%;
    height: 100%;
}

 #map0, #map2 {
    width: 33%;
 }

#map0, #map1, #map2 {
    margin: 5px;
    border:solid 4px #224a54;
    -moz-border-radius: 4px;
}


texts: Minimizing number of nodes in polygon using QGIS?

texts:
I'm building a web application that will be transferring polygons from server to client. Initially, the shapefile is coming from an external source that I have no control over.

Often, I find polygons with way too many points to describe a straight line. The result is, the size of the payload increases and therefore the fetch takes longer. This obviously isn't an issue with only a few polygons, but it will add up as more shapes are added to my dataset and I'd like to optimize bandwidth.

Is there a plugin, or native QGIS function that will take some kind of tolerance/variance as input, and produce (approximately) the same polygon with much fewer actual points?


texts: Creating random points in QGIS with values from Attribute Table

texts:
I'm trying to create a visual map in QGIS 2.4 of the approximate distribution of journey origin points using a shapefile of UK postcode sectors, layer joined to a table with a column giving the number of journey origins in each sector. Using the Vector > Research Tools > Random Points menu tool, I can select the Input boundary layer and click the 'Use value from input field' radio button, but the pull-down menu to select the required field does not seem to work.
If I then specify the Output Shapefile name and target location , I get an error message - "Errors: The following features IDs couldn't be calculated due to conversion error: 335, 335".
To confirm, QGIS does display random points if I specify the density or a specific number of points using the first two 'Stratified Sampling' Designs for individual polygs.


texts: Combining Raster and Polygons using ArcGIS Desktop

texts:
I am struggling with GIS (Esri).

I have the following overlaying data:

1) Several Polygons of suitable areas (m²) for the deployment of Photovoltaic
2) A TIF-raster (asc also available) of possible Photovoltaic Power Output (kWh/kWp) in different areas which is based on average sunshine duration per year.

I would now like to combine both 1) and 2) in order to see which Polygons/areas exhibit the highest possible Photovoltaic Power Output.

How can I combine both layers?

I'm clueless. For example, the Power Output for the upper polygon would be much higher than for the lower ones (see picture).


texts: Importing MaxSea 3D Database to PostGIS?

texts:
I am currently working on a project, where I combine several MaxSea 3D databases, which contain depth measurements made by the echosounder onboard on a ship. The idea is to combine these databases to create a detailed seafloor map around the Faroe Islands.

I do not know if any of you have worked with this program before, and I am not able to find anything on this site related to MaxSea.

The file format is proprietary, and the files generated are named 3DDB.PTZ. By inspecting the file, is see that it is a password protected Microsoft Jet database, but with some tools, I have been able to read from the file. The database contains a table named 3DTable which has the following structure:

╔═══════╦══════════╦═══════════╦══════════════════╦═════════════════╦═════════════╦════════════╗
║ Level ║   Lat    ║    Lon    ║ WeightDepth_CMap ║ ValueDepth_CMap ║ WeightDepth ║ ValueDepth ║
╠═══════╬══════════╬═══════════╬══════════════════╬═════════════════╬═════════════╬════════════╣
║     7 ║ 88021504 ║ 220475264 ║                  ║                 ║           1 ║ -48.64608  ║
║     7 ║ 88021504 ║ 220475392 ║                  ║                 ║           1 ║ -48.64608  ║
║     7 ║ 88021504 ║ 220508928 ║                  ║                 ║           3 ║ -19.2024   ║
║     7 ║ 88021504 ║ 220509056 ║                  ║                 ║           2 ║ -13.35024  ║
║     7 ║ 88021632 ║ 220474880 ║                  ║                 ║           1 ║ -48.64608  ║
║     7 ║ 88021632 ║ 220475008 ║                  ║                 ║           1 ║ -48.64608  ║
║     7 ║ 88021632 ║ 220475136 ║                  ║                 ║           1 ║ -48.64608  ║
║     7 ║ 88021632 ║ 220509056 ║                  ║                 ║           3 ║ -21.0312   ║
║     7 ║ 88021632 ║ 220509184 ║                  ║                 ║           3 ║ -21.9456   ║
║     7 ║ 88021632 ║ 220509312 ║                  ║                 ║           1 ║ -7.3152    ║
║     7 ║ 88021760 ║ 220461696 ║                  ║                 ║           1 ║ -70.77456  ║
║     7 ║ 88021760 ║ 220474624 ║                  ║                 ║           1 ║ -51.75504  ║
║     7 ║ 88021760 ║ 220474752 ║                  ║                 ║           1 ║ -51.75504  ║
║     7 ║ 88021760 ║ 220509952 ║                  ║                 ║           2 ║ -14.6304   ║
║     7 ║ 88021760 ║ 220510080 ║                  ║                 ║           5 ║ -33.8328   ║
║     7 ║ 88021760 ║ 220510208 ║                  ║                 ║           2 ║ -13.53312  ║
╚═══════╩══════════╩═══════════╩══════════════════╩═════════════════╩═════════════╩════════════╝

The column Level ranges from 5 to 21, in this particular dataset. However, I do not know how to process the data, as I do not know what format/projection the columns lat/lon is, and how to make sense of this.

In the project (it is my Bachelors project in Software Engineering), I will focus on the data quality, and how I can combine different datasets from ships and trawlers to create a detailed map. I have reference data from a research ship belonging to the Faroese Marine Research Institute, taken with a calibrated echo sounder of a higher quality.

To assist me in the process of combining the data, I intend to create a open source ptz2pgsql utility, similar to the shp2pgsql utility, and I would appreciate it, if someone of you could help me make sense of the data I have available at hand. I am able to send a sample of the data to you (a PTZ database). I have also contacted MaxSea to ask if it is legal to reverse engineer their proprietary file format.


texts: How to correctly rasterize a polygon which crosses 180E longitude?

texts:
I want to rasterize a polygon wich crosses 180E lon with gdal_rasterize.

The polygon is in the WGS84 projection.
I used the following command:

gdal_rasterize -burn 1 -of GTiff -ot Byte -ts 1500 1500 in.shp out.tif

But resulting raster is looks like this:

I tried to use different projections etc. but result is the same.

How can i do it in 'proper' way? Any solution will be good, not just gdal.


texts: How can i topologically sort 2D triangles and put in an array?

texts:
I have a set of 2D triangles in a Geometry array which together forms a
connected strip when placed them according to edge adjacency. However, in
the array they are not in topological order based on the edge adjacency.
How can i topologically sort them and put in the array?


texts: Where can I find High Quality Road Network Information

texts:
I am carrying out a research on how complex city road networks get and I need to get accurate map information.

I started with openstreetmaps however I ran into a couple of issues: some of the dumps contain duplicate information - for example two ways which have the same name tag have distinct way ids and this skews my representation.

Moreover, when I converted their data into a graph, the graph was not connected, infact about 30% of the roads are not connected to the giant component. Could this be a mistake in my representation or something to do with openstreetmaps?

I have also tried dumps from the US geological survey however roads are split into segments without enough information to recombine them again. This makes it impossible to get full road information.

Can you help with road network data suggestions?


texts: Split label style in GeoServer CSS or SLD

texts:
I would like to split the style of my labels, which consist of 2 attributes.
I already managed to separate the lines, but I need to style them differently as well, as the second line should be bold.

Currently my map looks like this:

And I would like it to look like this:

I would prefer to do it in CSS, but I also prepared the style in SLD, so I would be happy for a solution in either of them.

My code looks like this in CSS:

* {
    label: [IDENT][LOWERLIMIT];
    -gt-label-auto-wrap: 1;
    font-weight: normal;
    font-fill: black;
}

and like this in SLD:

<Rule>
  <TextSymbolizer>
     <Label>
       <ogc:PropertyName>IDENT</ogc:PropertyName><![CDATA[
]]>
       <ogc:PropertyName>LOWERLIMIT</ogc:PropertyName>
      </Label>
      <Font>
       <CssParameter name="font-fill">black</CssParameter>
       <CssParameter name="font-weight">normal</CssParameter>
      </Font>
  </TextSymbolizer>
</Rule>


texts: Rotating view without losing labels in QGIS?

texts:
Using QGIS 2.8 Wien and a beautiful shapefile from OpenStreetMap.

When I rotate the map (either in the print composer or in the main QGIS panel) the labels with the street name disappear. I've tried different solutions from this post, but it didn't work in my case.

Is there a way to rotate the map and keep the labels?


texts:
I am looking for help to modify mouse click so that I can click a line segment without having to be directly over it.

I have been looking at Leaflet.AlmostOver but have been having difficulty getting it integrated.

Has anyone done this and have any suggestions on how I might resolve my issue?


texts: Applying layer file not formating labels in ARCGIS 10

texts:
I'm applying layer files in ARCGIS 10 using the ApplySymbologyFromLayer_management command via a python script tool.  The symbology part works fine but the lables are not formatting.

It is my understanding that applying a layer file should also format the lables as well.  Has anyone else had this issue?  Any ideas would be appreciated.

Here is the code I am using.

import arcpy, os
from arcpy import env

## Create parameters for modelbuilder user specified fields.
symbolize_layer = arcpy.GetParameterAsText (0)
symbol_lyrFile = arcpy.GetParameterAsText (1)

arcpy.AddMessage("**** Step 1 ****")
#  apply symbology
arcpy.ApplySymbologyFromLayer_management (symbolize_layer, symbol_lyrFile)

arcpy.AddMessage("**** Step 2 ****")
# refresh the view
arcpy.RefreshActiveView()
arcpy.RefreshTOC()

arcpy.AddMessage("**** OK now were Done ****")


texts: Exporting satellite image raster causes change in pixels using ArcGIS Pro

texts:
I have a WorldView-2 satellite raster which I re-project from WGS84 to WGS84/UTM using Raster Project. I wish to export this from outside the geodatabase so I can share it with colleagues. However, when exporting the re-projected raster to my desktop (via right-clicking on it from the contents pane -> export), the exported raster has a changed pixel value range and white stripes appear where there was previously no data. I used the default settings when exporting.
The white background is caused by areas beyond the original raster extent (i.e., areas of no-data in the original raster) being assigned high-pixel values (65535) when exported. To resolve this, I've tried to set the "NoData Value" to zero. This removes the white background stripes. However, it also causes the exported raster's pixel range to begin from 1 rather than 0.
How can I export my raster image without creating pixels beyond the original extents or changing the pixel value range?


texts: ArcObjects: Getting the properties of a layer that has lost its data source (VB.NET)

texts:
When you save an ArcGIS map document (MXD) and then change the data source of a layer in the MXD, you can no longer see the layer in the MXD and you must redefine its data source.

But suppose you need to discover the previous data source's feature class name. In ArcMap, you can find it under the layer's Properties -> Source Tab, but I want to access this information programmatically.

I tried the code below but it doesn't work for this kind of layer. Can anyone suggest another method for retrieving the properties of a layer whose data source has been lost?

ctype(mLayer,IFeatureLayer).aliasname


texts: Convert a number in attribute table to binary

texts:
I would like to convert whole integers from decimal to binary in ArcGIS using the field calculator.

I am looking for the equivalent of DEC2BIN excel command.

All values are 255 or below so 8 bit is fine.

It must have 0 padding to the left. So for example 4 must be 00000100.

I can do it by exporting the attribute table to Excel and using the Dec2Bin command in Excel and rejoin the data but I would like it in a model so want to use the field calculator.

Dec2Bin (Value, 8) works in Excel.

Attached is how I would like it to look when complete.

So basically change the value to and 8 digit binary value stored as a text field with leading 0s so all field have 8 values.

I found this code in arcscripts but cannot get it to function.

Dim valI as integer
Dim valB as string
Dim x as integer
valI = [fieldName]
for x = 7 to 0 step -1
  if valI >= 2^x then
    valB = valB & 1
    valI = valI - 2^x
  else
    valB = valB & 0
  end if
Next
valB = valB & "00000000"

fieldName =
valB


texts: I have a shapefile with 634k features that when I open in QGIS is filtered to 6 squares. How can I expand it to get all the points, I need to make a topo map from it.


texts: Error "No module named PyQt5.QtCore" in QGIS standalone application

texts:
import os, sys
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\qgis\python')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\qgis\bin')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\Python37')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\Python37\Scripts')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\Qt5\bin')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\Python27\Scripts')
sys.path.append(r'C:\Program Files\QGIS 3.6\bin')
sys.path.append(r'C:\Windows\System32')
sys.path.append(r'C:\Windows')
sys.path.append(r'C:\Windows\System32\wbem')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\Python37\Lib\site-packages\pywin32_system32')
sys.path.append(r'C:\Program Files\QGIS 3.6\apps\Python37\Lib\site-packages\numpy\.libs')
import qgis

I would like to write a standalone QGIS application. My Python file is on the desktop and I would like to access the QGIS functions. At the beginning I had the error import qgis module not found, which could be solved by the first line of sys.path.append. Now I get the error:
from PyQt5.QtCore import *
ModuleNotFoundError: No module named 'PyQt5.QtCore'

Although I have appended many different paths to the environment variables, I keep getting this error.
How do I solve this error and use the QGIS functions in a standalone application?
Error trying the anwser:


texts: How to disconnect a signal when untoggling plugin icon?

texts:
Following on from this question,  I've created a signal whereby when the plugin icon is toggled (or checked), the user can change the visibility of a layer which is connected to a function (print_me). This will print a simple message.

However, I've tried to disconnect this signal when the plugin icon is untoggled but I receive the following error:

  TypeError: 'function' object is not connected

Even though the function is connected as the message is still printed. How could I solve this?

class Example:
    def __init__(self, iface):
        ...

def run(self, checked):
    #################
    if not self.pluginIsActive:
        self.pluginIsActive = True
        if self.dockwidget == None:
            self.dockwidget = Example_dockDockWidget()
        self.dockwidget.closingPlugin.connect(self.onClosePlugin)
    #################
    root = QgsProject.instance().layerTreeRoot()
    group = root.findGroup('Group')

    def print_me():
        print 'Message!'

    if checked:
        for layers in group.children():
            layers.visibilityChanged.connect(print_me)
    else:
        for layers in group.children():
            layers.visibilityChanged.disconnect(print_me)

Edit:

Although a solution I posted involved putting the print_me function outside the class, I need this function to be inside the class. This is because, in my plugin, the function makes references to the QGIS interface such as printing values to the plugin's dockwidget.

If I incorporate a try method such as the following:

if checked:
    for layers in group.children():
        layers.visibilityChanged.connect(print_me)
else:
    for layers in group.children():
        try:
            layers.visibilityChanged.disconnect(print_me)
        except TypeError:
            pass

The error is ignored but if I toggle the icon again, the connection is made again. So it runs the print_me() function twice! I don't understand why a connection can be made when the icon is checked but when unchecked, it cannot be disconnected...


texts: Adding points to every "eligible" cluster in PostGIS (ST_ClusterDBSCAN)

texts:
I'm currently using ST_ClusterDBSCAN in PostGIS to identify similar points for conflation.
(Using clustering to see where multiple entries are within a given range and then manually reviewing & deduplicating.)
This is working well most of the time, but I noticed on the docs for the function it says:

Note that border geometries may be within eps distance of core geometries in more than one cluster; in this case, either assignment would be correct, and the border geometry will be arbitrarily asssigned to one of the available clusters.

So I think I'm currently missing the case where, e.g. there are points A, B and C evenly spaced 50 meters apart — the function could give me results of:

(A, B) / (C)

or

(A) / (B, C)

With both cases equally valid from PostGIS point of view. For my use-case though, if B & C were the duplicates in #1, I currently wouldn't find them.
Is it possible to somehow "duplicate" the border geometry here (i.e. something akin to a LEFT OUTER JOIN vs an INNER JOIN), so it gets added to any 'eligible' cluster to give something like:

(A, B) / (B, C)

I don't think that's possible from ST_ClusterDBScan alone from reading the docs, but wondering if there's a better approach in this case.


texts: Updating Location of Point Features using XY Fields?

texts:
I have changed the coordinates of couple points in the attribute table and I want to update the location. I have found the script that is working in field calculator.

Does anybody know how to write it using a Python script in ArcGIS Desktop (not in Field Calculator)?

def XYsetVALUE( shape, X_value, Y_value):
   point = shape.getPart(0)
   point.X = X_value
   point.Y = Y_value
   return point

__esri_field_calculator_splitter__
XYsetVALUE ( !SHAPE!, !X_COORD!, !Y_COORD! )

COORDINATE FIELD TYPE ------ STRING !!!!


texts: Cannot save map from composer as image or pdf

texts:
I have created a map and made a nice lay-out in the printer composer (in QGIS 2.4), the weird thing is, that if I try to save my map as an image or as a pdf immediately after entering the name in the "save as-dialog box" my map disappears from the composer interface, and then it is nowhere to be found.  Luckily I did save the template, and only after restarting QGIS can I load the template again.

I don't understand what is going on here and would be very happy if someone can tell me this.


texts:
How would you label this map logically?

I will preempt the question, with the statement it is for a fire truck, so we really have to worry about seeing the address numbers at 3 am with one contact in. I have tried a few different options and seek any suggestions.


texts: ArcGIS/arcpy: Geometry object comparison methods not working. Do they take z values into account?

texts:
I have an arcpy script that makes calculations based on an input PolygonZ feature class.  What I need to do in several places in the script is test whether a point falls within a polygon in 2D.  In these cases I have a 3D geometry object for an input polygon and a 2D geometry object for a test point that I create on the fly.  I use the pntgeom.within(polygeom) comparison.  The problem is that this statement always evaluates to False.  So, I am wondering if the comparison is actually working in 3D (and implicitly giving the point a Z value of 0), in which case the point will not ever be "within" the polygon.  The polygons are 3D surfaces raised above the ground.  The whole procedure is intended to determine which of two possible sizes the current surface possesses.  I create a test point that would fall only within the larger size, and, if the point is not within (2D) the current surface polygon, then the size must be the smaller one.  Specifically my code chooses between two different arc radius lengths.  Here is the relevant code snippet:

pntTest = arcpy.Point()
pntTest.X = coords_A[0] + (12200.0 * uCLA[0])
pntTest.Y = coords_A[1] + (12200.0 * uCLA[1])
pntGeom = arcpy.PointGeometry(pntTest)
testpnt = pntGeom.centroid
arcpy.AddMessage('  Test point coords:  X = %s   Y = %s' % (str(testpnt.X), str(testpnt.Y)))
where = '"OISSURTY" = \'%s\' AND "RWYDESG" = \'%s\'' % (surfInfo[0], surfInfo[2])
arcpy.MakeFeatureLayer_management(oids, 'oidsurfs', where)
cnt = arcpy.GetCount_management('oidsurfs').getOutput(0)
arcpy.AddMessage('  No. featurs in layer "oidsurfs" = %s' % str(cnt))
surf_sc = arcpy.SearchCursor('oidsurfs', "", spatref)
shapefieldname = desc_oids.shapeFieldName
for row in surf_sc:  #should be only 1
    surfGeom = row.getValue(shapefieldname)
del row
del surf_sc
arcpy.Delete_management("oidsurfs")
len = surfGeom.length
arcpy.AddMessage('  surfGeom length = %s' % str(len))
if pntGeom.within(surfGeom):
    radius = 10000.0
else:
    radius = 5000.0
arcpy.AddMessage("  Distance of opnt from rwy end = %s" % str(dist))
arcpy.AddMessage("   Inner arc radius (conical surface) = %s" % str(radius))

I have checked the coordinates of the test point and confirmed that it lies within (2D) the current polygon.  This also confirms that it is a geometry object.  I also call the length property on the polygon geometry to confirm whether it too is a geometry object; it is.  The current data I am testing has all polygons being of the larger size, so the code should choose the larger radius (which is used in a later calculation).  The problem is that the geometry object comparison evaluates to False:  no errors - just doesn't work.  I thought these were 2D comparisons, but now I'm beginning to wonder.  I just don't get it.  Any ideas?

ArcGIS-ArcView 10.0
Windows 7
8 GB RAM


texts: Passing ModelBuilder parameter to update cursor gives ERROR 999999?

texts:
I have two point layers.  I trust the geometry (read: x and y coordinates) from one of them, but not the attributes data.  I trust the attribute data but not the geometry from the other one (the attribute tables have different schema).  The basic process (shown in a model builder screenshot below) is to spatially join the point layers, join the result to the points with good data/bad geometry and update their shape field with a script. The script produces one of those damnable "999999: Error executing function" errors.

If I manually step through every stop of the process and feed the update cursor a ready-made feature class it works fine, the points move as intended.  It's just when it's set up to run in a model that the script throws the error.  I have a feeling that I may be passing the parameter incorrectly, but I can't say for sure.

I've tried moving the "fc = GetParameterAsText()" line around (including in the update cursor parameters) to no avail.

I will also point that I shamelessly took the python from a post on Geonet, located here.

The Model:

The Script:

import arcpy
from arcpy import env

env.workspace = r"C:\Temp\TempGDB.gdb"

#Get the feature layer from the previous step in the model
fc = arcpy.GetParameterAsText(0)

pnt = arcpy.Point()

#Update the shape field with the new coordinates
rows = arcpy.UpdateCursor(fc)
for row in rows:
    pnt.X = row.X_Coord
    pnt.Y = row.Y_Coord
    row.shape = pnt
    rows.updateRow(row)

del row, rows

The Script Tool Properties


texts: Image resolution and image scale of air photo

texts: I was wondering if 2 aerial photos of different scales can have the same resolution. My understanding is they cannot as fine resolution images will only cover a small area and vice versa?


texts: Is there a way to "disperse markers" across multiple point layers?

texts: I have a bunch of point layers which often have overlapping features. The disperse markers geoprocessing tool works within a layer but is there a way to do it across more than one layer, so none of the points are overlapping? Using ArcMap 9.3.


texts: Projection Discontinuties/NoData value errors in google earth engine export

texts:
Below is an example of a larger process (downloading MODIS data-- multiple products for 23 cities between 2004-2016) that is giving me a bit of a hassle:

#load libraries
import ee
import datetime

#Init earth engine and load a collection
ee.Initialize()
imgs = ee.ImageCollection('MODIS/006/MOD11A2')

#filer by some dates and extract a single "band" and extract one image. I'm sure there is a cleaner way to do this
dat = imgs.filterDate('2004-01-01','2004-01-02').select('QC_Day').toList(1).get(0)

#print some info
print dat.getInfo()

#define a bounding box
bbox = (38.94993166800003, -7.127943276999929, 39.44200277300007, -6.552704129999938)

#Create the bounding box as an earth engine object thing
target_projection = ee.Image(imgs.select('QC_Day').first()).projection()
thecity = ee.Geometry.Rectangle(bbox).transform(target_projection,250) #hoping the 250 is in meters

#For some reason ee sees dat as a computed object and so clip won't work
#trick it
dat = ee.Image.cat([dat])

#print the info again
print dat.getInfo()

#clip the image to the city
dat = dat.clip(thecity)

#make the ask to EE
ask = ee.batch.Export.image.toDrive(image = dat, description = 'hello')

#start the ask
ask.start()

When the download completes and I move to R (my preference) to visualize the image:

#Import the raster package, analgous to import raster
library('raster')

#load the raster
ras = ras("C:/Users/Me/Downloads/hello.tif")
#Set all values that equal 0 to 1 and everything else to 0
zeros = ras == 0
#make a crappy map
plot(zeros)

I get the following map. My problem is the green area on the sides that looks like some combination of a projection + no data issue. The raw MODIS data (downloaded from LPDAAC) does not have this pattern. Thoughts on how to make the green go away?


texts: Where can I find ArcGIS/ESRI sample graphics?

texts: I am looking for sample graphics/pictures that I can use in my program for pinpoints and other markers. Does ArcGIS/ESRI have any that are available to the public? I'm trying to use PictureMarkerSymbol and I would like to use a url to the graphic/picture location if there is one.


texts: Use SDE.ST_GEOMETRY functions in a custom function

texts:
I have a dummy query where I generate an SDE.ST_GEOMETRY point:

select
    shape as srid_test
from
    (
    select
        sde.st_geometry ('point (10.01 50.76)', 4326) as shape
    from
        dual
    )

I've written a function that returns the SRID of the point (it's just a test function--it doesn't have a practical purpose yet).

The function works as expected.

with function test_function (shape_in in sde.st_geometry) return number is
  srid number;
begin
--    srid := sde.st_srid (shape_in);
    select sde.st_srid (shape_in) into srid from dual;
    return srid;
end;

select
    test_function(shape) as srid_test
from
    (
    select
        sde.st_geometry ('point (10.01 50.76)', 4326) as shape
    from
        dual
    )

You'll notice that I use a query to assign the SRID to the SRID variable:

select sde.st_srid (shape_in) into srid from dual;

Instead of using a query to get the SRID, I want to simply use the ST_SRID function:

 srid := sde.st_srid (shape_in);

However, I'm getting an error when I switch from the query to the function:

  ORA-06553: PLS-548: invalid use of operator.

with function test_function (shape_in in sde.st_geometry) return number is
  srid number;
begin
    srid := sde.st_srid (shape_in);
--    select sde.st_srid (shape_in) into srid from dual;
    return srid;
end;

select
    test_function(shape) as srid_test
from
    (
    select
        sde.st_geometry ('point (10.01 50.76)', 4326) as shape
    from
        dual
    )

Why am I getting this error?

(Oracle 18c, 10.7.1 Geodatabase)


texts: Use of ol.source.BingMaps culture parameter

texts:
I have noticed that even though openlayers 3 include the culture parameter when requesting the Bing maps layer, Bing maps itself does not support imagery for other languages except English (supported codes).

Is there a way to set the display language of Bing maps roads and places in openlayers 3 other than English and if so, how? I am interested in the Greek language

Here is a fiddle to try your coding.


texts:
We are attempting to port a Python-based CartoDB authentication sample app to NodeJS:
The example, "Basic auth middleware with CartoDB," was created by CartoDB labs.

We have sucessfully replaced one section of the Python code which authenticates to cartodb and makes SQl API calls by utilizing the "cartodb-nodejs" client module.

However, we are stuck on another section which executes a PUT to CartoDB Maps API endpoint, in order to create a Named Map (we have an on-premise Cartodb instance, so we utilize a custom endpoint URL).

What we need to do is to POST a JSON-formatted template (stored in the "named_map" variable, below) to our on-premise endpoint.

var named_map ={
    "version": "0.0.1",
    "name": "testmap",
    "auth": {
        "method": "token",
        "valid_tokens": ["password"]
    },
    "placeholders": {
        "color": {
            "type": "css_color",
            "default": "red"
        },
        "filter": {
            "type": "number",
            "default": 1
        }
    },
    "layergroup": {
        "version": "1.0.1",
        "layers": [
            {
                "type": "cartodb",
                "options": {
                    "cartocss_version": "2.1.1",
                    "cartocss": "#layer { polygon-fill: <%= color %>; }",
                    "sql": "select * from v1_8_walgreens"
                }
            }
        ]
    }
}

We can accomplish the equivalent task via cURL and a stored/file-based template, "template.json".

 curl -X POST  -H "Content-Type: application/json"  -d @template.json  "https://cartodb.brighterdevelopment.com/user/demo-admin/api/v1/map/named?api_key=API_KEY_GOES_HERE"

Below is the actual Python code section we need to duplicate in Node. NOTE: the endpoint, api_key and named map name are basically passed-in as variables/configs.

requests.put(os.path.join(config.get('cartodb', 'maps_endpoint'), "named", map_name),
         data=json.dumps(named_map),
         params={"api_key": config.get('cartodb', 'api_key')},
         headers={'content-type': 'application/json'})

We have made a few unsuccessful attempts at implementing this with the Node “Request” Package.
We are open to any approach.


texts: How to automatically start a Torque animation using Carto Builder

texts:
I'm trying to figure out how to get my Torque animation to start automatically when my map is loaded. Here is a sample map:

Torque Map with Carto Builder

I know this was possible with the old Carto Editor (see here) but I can't figure out how to do it with the new tool.


texts: OpenLayers 6: blur effect with ImageStatic

texts:
I'm using a static image on my map in a project with Django. The code is inspired to the official example:

var imageSourceExtent = [
  12.1986543459999996,
  44.9052984770000023,
  12.3448225140000005,
  44.9825835529999978,
]

var imageSourceProjection = new ol.proj.Projection({
  code: 'EPSG:4326',
  extent: imageSourceExtent,
});

var imageSource = new ol.source.ImageStatic({
  url: '{% static 'gis-data/raster/clip_NDVI_color4326.png' %}',
  projection: imageSourceProjection,
  imageExtent: imageSourceExtent,
});

var imageLayer = new ol.layer.Image({
  source: imageSource,
});

var outDoors = new ol.layer.Tile({
  type: 'base',
  title: 'Topographic',
  visible: true,
  source: new ol.source.XYZ({
      attributions: 'powered with OpenLayers',
      url: 'https://api.mapbox.com/styles/v1/maxdragonheart/...'
    }),
});

var map = new ol.Map({
  layers: [
      outDoors,
      imageLayer,
    ],
  target: 'map',
    render: 'canvas',
  view: new ol.View({
    center: ol.proj.fromLonLat(ol.extent.getCenter(imageSourceExtent)),
    zoom: 12,
  })
});

I've used QGIS to create the png from a tif and inside QGIS I see the image in this way:

But when I see the same image on the map it appear in this way with a blur effect:

Why happen this? How I can solve?


texts: Google Earth Engine - Minimum of two Linear Equations

texts:
I am trying to compute scaled Land Surface Temperature from Modis per this equation:

Does anyone have any suggestions on how to code this in GGE? I'm quite new to java. I have no trouble creating both equations separately, but am curious as to how I might combine them in one function.


texts: Add r.null functionality to GRASS toolbox algorithms in QGIS

texts: I'd like to be able to use GRASS's r.null operation in the QGIS processing toolbox (not via the GRASS plugin). However, it doesn't seem to be available. How would I go about adding code to allow for the r.null algorithm?


texts:
I want to find optimum antenna location for an area. Every antenna can cover an circular area. These areas diameter is 750meter. So how can i find minimum number of antennas and locations.
NOTE: I use QGIS. Blank areas are not important. Only buildings must be covered.

To better understand here is example images...

Here is building layer

Firstly i generated 50m grid of the layer.

After that i intersected it with building layer. So i can eliminate blank areas.

I stucked in here. Which tool can help me determine location of  antennas (750meter diameter). I know it can not show me %100 optimum locations. Bu i need possible locations. Here is another screen that i made it myself without analysis tools. The result must look like that. Thanks anyway.


texts:
I have a problem with creating tiles from from Mars HiRISE images using gdal/gda2tiles. Can't find info regarding this topic, so, maybe someone already faced that issue?

Image looks like this

Final goal is to show this tiles on the top of CesiumJS, but for now Leaflet/OL/etc map would be also fine.

Input source file has following parameters:

What worries me a bit here is the fact, that

Datum is unknown. Maybe this is the source of all my problems with further processing, but this is just my guess.

QGIS show following SRS:

For tiles generation I use gda2tiles utility.

If I try to run gdal2tiles directly from Python like this

minZoomLevel = 8
maxZoomLevel = 12
options = {'zoom': (minZoomLevel, maxZoomLevel)}
gdal2tiles.generate_tiles(sourceImagePath, outputTilesFolder, **options)

I got following output:

So, for sure something is wrong.

I tried different resampling values, tried to provide also source Spatial reference, but that didn't help.

After googling I also tried before calling gdal2tiles to wrap original image and only then call gdal2tiles.
So, I tried

gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -t_srs EPSG:3857 -ot Int16 -r cubic /MarsData/esp_057760_2025_red.jp2 /MarsData/warped.tif

and then call gda2tiles for warped.tif, but that didn't help.

I also tried to provide input file SRS for wrapping (took it from QGIS)

gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -s_srs '+proj=eqc +lat_ts=0 +lat_0=20 +lon_0=180 +x_0=0 +y_0=0 +a=3393833.2607584 +b=3393833.2607584 +units=m +no_defs' -t_srs EPSG:3857 -ot Int16 -r cubic /MarsData/esp_057760_2025_red.jp2 /MarsData/warped.tif

But this also didn't help. I still see the same output tiles.

Can someone point me where is my problem and how to solve it?


texts:
Suddenly I got a problem opening the project I was working on since a long time.
When I open the project, QGIS crashes.

I think the problem might be in one layer, but as I cannot open the project I neither can delete this layer in QGIS desktop. Is there any option to do this another way?


texts: Creating a Radial Sector Overlay in QGIS

texts:
I need to create a sector grid radiating out from a point of interest. Here is an example that I laboriously put together by hand:

I have found a couple of similar solutions in How to generate dartboard shapefiles using ArcGIS for Desktop? and Draw wind rose with QGIS from PostGIS but neither exactly fit my criteria.

Ideally I need each sector to be an individual polygon (for data analysis purposes) and to be able to create this overlay for different sites using parameters to control the number of sectors and radial distances.

Can anyone suggest or provide a means to generate these shapes in QGIS easily? I rather presume that this is a relatively common task so don’t want to go to great effort to reinvent the wheel but I am struggling to find a solution and being a complete newcomer it would take me a very long time to come up with a solution from scratch.


texts: Getting azimuthal (direction) values of raster pixels?

texts:
I have extracted urban footprints using Landsat ETM+. I am now interested to make rose diagram depicting direction of each pixel to ascertain in which direction urban area has been increased over the study period.

Obviously, this need reference point from where direction of each pixel will be calculated. I can get the results while converting the pixels into polygon and clip the features in four cardinal directions (0-90, 90-180, 180-270 & 270-360 degrees) and plot the resultant data as rose diagram.

But, is there any elegant way to retrieve better results while converting each pixel to assign its direction ranging from 0-359?


texts: Cannot access items array property of MapView.allLayerViews - ArcJS API 4.4

texts:
See below:
require(["esri/Map","esri/views/MapView"], function(Map, MapView) {
    var map = new Map({ basemap: "streets" });

    var view = new MapView({
        container: "viewDiv",
        map: map,
        center: [-97.744483, 30.272705],
        zoom: 15,
    });

console.log(view.allLayerViews.items);

console.log(view.allLayerViews.items[0]);

console.log(view.allLayerViews.length);

I suppose the array calls are accessing the Array prototype, but I don't know how to make it not do that. I would like to be able to loop through the layerViews with forEach().


texts:
I'm running this script to iterate through layers, zoom to them, then export each as a pdf but I keep getting this error:

AttributeError: 'str' object has no attribute '_arc_object'

I'm not skilled enough at python to interpret this error message.

Here's my code:

import arcpy
def zoom_export(zlayer,pdfname,ratio):
    zoom = arcpy.mapping.Layer(zlayer)
    arcpy.SelectLayerByAttribute_management(zoom)
    df.extent = zoom.getSelectedExtent()   #zooming to selected feature layer
    df.scale = df.scale * ratio
    arcpy.RefreshActiveView()
    zoomin = arcpy.mapping.ExportToPDF("C:\Users\TSPA standard\Desktop\gis\Test.mxd", "TESTPDF")
    arcpy.SelectLayerByAttribute_management(zoom, "CLEAR_SELECTION")

def turn_on(layer,tf,mxd):
    mxd = "C:\Users\TSPA standard\Desktop\gis\Test.mxd"
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    layers = arcpy.mapping.ListLayers(mxd, "*", df)
    for x in layers:
        if x.name == layer:
            print x.name
            x.visible = tf
            arcpy.RefreshTOC()
            arcpy.RefreshActiveView()

mxd = arcpy.mapping.MapDocument(r"C:\Users\TSPA standard\Desktop\gis\Test.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
layers = arcpy.mapping.ListLayers(mxd, "*", df)
for layer in layers:
    arcpy.MakeFeatureLayer_management(layer,"blah")
    turn_on(layer.name,True,mxd)
    pdfname = layer.name+"PDF"
    zoom_export("blah",pdfname,1.1)
    turn_on(layer.name,False,mxd)


texts: Making a line from points but preserving classification of points in QGIS?

texts:
I want to classify a gpx file by elevation in QGIS.

This is what I have done:

Imported track_points from gpx file (because only this part of file containes elevations)
Saved that file as shp (just in case)
Classified shp as categorized

Now I would like to connect points as line but I would like to preserve categorization color of points. I have tried the plugins Points2one and Recostruction Line. But I haven't had luck and I also got errors in some cases.

Images below

Left: what I have achieved Right: created line by IDs with Points2one plugin (colors are not preserved)

EDIT

After suggestion given by mgri I end up having this in SAGA:

* I have tried to choose couple of combinations

Result:


texts: Does Leaflet layer 'add' event not fire properly?

texts:
A popup won't open until the data file has loaded and the layer is added to the map - that's pretty obvious. I figured it would help to bind it to the 'add' event - but that never fires:

var l1 = omnivore.kml('data/Niveauregulering01a.kml');
var pop = '<img src="content/img/whatever.jpg" alt="..." style="width:300px;"/>'
l1.on('add', function(){
  console.log('add');
  l1.openPopup();
});
l1.bindPopup(pop).addTo(map);

There seems to be no other event that I can bind it to:
http://leafletjs.com/reference.html#path-add

Any suggestions?


texts: How do you export a 3d-model in QGIS to a file that can be read in Rhino?

texts: I am really new to QGIS and haven't really gotten familiar with all the functions yet. I've made a 3D map/model in QGIS, and would like to export this in a file format so that I can use it as a terrain model in Rhino. How do I do this?