texts: Calculating area of an EPSG:4326 Polygon

texts:
Possible Duplicate:
  What is the most accurate coordinate system for calculating areas of polygons?

In our Postgres/Postgis database we store polygons in a column defined as lonlat/EPSG:4326, as it may contain data from anywhere in the world.

We have a requirement to calculate the surface area (i.e. in square meters or hectares). The method ST_Area can perform the calculation but returns the answer in the used coordinate systems units. That is, it returns the area as what I believe is "square degrees". As the size of a degree is variable depending on the location in the world, this does not make much sense.

One solution is to transform (ST_Transform) the polygon to another crs before calculating the area, perhaps the appropriate UTM32 zone, but that will require us to dynamically determine the zone first. Alternatively we could always transform to the Google Mercator (EPSG:900913), but we fear for the precision.

Any inputs?


texts: Is it possible to use non-printable characters as delimiter for text file layer in QGIS?

texts:
I am trying to import a Delimited Text File Layer from a .csv file in QGIS. Because this caused a lot of errors, I posted a question here and got the hint towards using ASCII non-printable control characters as delimiters. These are described on Wikipedia and in this blog post.

Since I am using Tweets, for which it is very difficult to find a convenient delimiter because nearly all characters are used in them, I really want to give this a shot. However, I cannot find out how to import a text file into QGIS specifying these characters as delimiters. I tried the different spellings that are given in Wikipedia as custom characters and as regular expression.

0x31
0x1F
^_

None of these work. Interestingly, the lines are detected correctly - only the record delimiters are a problem.

Any ideas?


texts: Formulae to convert WGS84 coordinates (B,L) to spherical coordinates (phi,lambda)?

texts:
How to convert WGS84 coordinates(B,L) to spherical coordinates(phi,lambda)?

I Need Formulae?

the thing I'm wondering seems quite simple but I've been searching through internet and library of my school about 3 days. I've found nothing about this. Maybe someone have done this before?


texts: Converting factfinder2.census.gov map data to KML?

texts:
I've downloaded a map from http://factfinder2.census.gov/faces/tableservices/jsf/pages/productview.xhtml?pid=DEC_10_SF1_P1 (limited to San Francisco census tracts), and I'd like to convert it to KML format. However, after unzipping the shapefile .zip file and running

$ ogr2ogr -f KML pop_by_census_tract.kml 140_00.shp

I get the following error message:

ERROR 6: No translation for Mercator_Auxiliary_Sphere to PROJ.4 format is known.
Warning 1: Failed to create coordinate transformation between the
input coordinate system and WGS84.  This may be because they
are not transformable, or because projection services
(PROJ.4 DLL/.so) could not be loaded.
KML geometries may not render correctly.
This message will not be issued any more.

Source:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Mercator_Auxiliary_Sphere"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",0.0],
    PARAMETER["Standard_Parallel_1",0.0],
    PARAMETER["Auxiliary_Sphere_Type",0.0],
    UNIT["Meter",1.0]]

ERROR 1: Latitude 4553026.920500 is invalid. Valid range is [-90,90]. This warning will not be issued any more
Warning 1: Longitude -13627851.842700 has been modified to fit into range [-180,180]. This warning will not be issued any more

When I open the resulting KML file in Google Earth, I get randomly scattered lines across the entire planet. The example at code.google.com/apis/kml/articles/vector.html works fine, so the installation of GDAL (1.8.0 from macports) isn't totally broken.

http://trac.osgeo.org/gdal/ticket/3962 seems related, but I'm not sure if it's saying that Mercator_Auxiliary_Sphere is just not supported by GDAL, or if it's implying I should pass another command line option to make it work.

I wound up accomplishing my goal by downloading http://www2.census.gov/geo/tiger/GENZ2010/gz_2010_06_140_00_500k.zip, which contains shapefiles that ogr2ogr and shp2pgsql can understand, and joining it with the pure table data from factfinder2. I'm still curious about the format factfinder2 is using, but it's no longer blocking progress for me.


texts: Open QGIS project from PyCharm or other Python IDE

texts:
I have QGIS project and Python script.
I know how use QGIS Python console to run script.
Is there any way to open a QGIS project in PyCharm and to run the Python script there?


texts: Attribute Join creates headers but cells remain empty in ArcMap?

texts:
I'm trying to carry out an attribute join using one CSV and shapefile using the "LSOA11CD" fields as pictured below

The join validates in ArcGIS Desktop, however I get the following error message when I try to join the two tables.

The CSV contains significantly more entries than the shapefile field, however identical entries are present in each of the tables.


texts: What are the low cost options for converting FLEX API based web mapping applications to any other platform?

texts: As flex is reaching the end of its life, I have one project to convert Flex based web mapping application to any other possible solution with minimum development effort and secure way. I believe Dojo based ESRI javascript API (as Siverlight direction is not clear either) will suit the need here (If company's policy allows cloud and then ArcGIS Portal can be leveraged too). But there are some suggestions around Geocortex framework or any microsoft based approach. The belief is Geocortex framework can give a secure and objects based development approach. But in one of the project I worked with Geocortex and multiple web mapping application- I found the development effort was so much underestimated under similar assumptions and project was running into mess. Any suggestion/ opinion/ experience will be helpful here. There is another web mapping application based on leaflet based ESRI API with Angular JS used by the organization and some of the client need can be microscopic too. Will Dojo be able to meet the need here?


texts: Is it possible to position QGIS labels by clicking on polylines (and polygons)?

texts:
QGIS has now few really good options for custom labels but I can't find a way how to easily manage labels along polylines.

Is there a tool which helps me place labels by clicking on a polyline so the label will be anchored to that point on line - to automatically set position above/on the line and proper rotation or even better curved with the line?

Similarly with polygon - click on polygon to exact spot where I want the label to be.

I know there is a workaround to create columns for labels x/y/rotation and move and rotate labels manually but that is often very time consuming and a clunky way compared to some other GIS software.


texts: What are object Names of of Python Console ToolBars?

texts: I have a custom plugin that when loaded it hides all tools and toolbars not specified to load. How ever by doing this it hides all the tools that are present in the python console. Clean Console, Show editor, save script, run script etc. I have tried in many ways to track down what are the names of these toolbars so i could specify them to load with my plugin but just cannot find them. I am referring to the tool bars in the photo. Can anyone help me track these down?


texts: How to create Where Clause from user input via Python list?

texts:
I am trying to perform a query based on values in a Python list. The where clause will then be used in an Update Cursor, so that the rows NOT equal to the values in the list will be deleted.

This is what the list will look like after it is passed into the script from the user (using a multivalue parameter in the script tool), and the semicolon delimited string has been split. This is just an example of many choices the user can make.

userSelectionList = ["40","77","95]

selectionQRY = "RTE_1_NBR <> 'userSelectionList'" <-- I know this doesn't work

Hard coded, it would look like:

selectionQRY = "RTE_1_NBR <> '40' AND RTE_1_NBR <> '77' AND RTE_1_NBR <> '95'"

But the input and number of selection depends on the user's selection.

It should finally be something like this:

userSelectionList = ["40","77","95]
selectionQRY = "RTE_1_NBR <> 'userSelectionList'"

with arcpy.da.UpdateCursor(userI,["RTE_1_NBR"],selectionQRY) as userUC:
    for row in userUC:
        userUC.deleteRow()
del userUC

Is it possible to make the where clause from a user input/dynamic Python list? Thank you!


texts: Jupyter and ArcGIS error: arcpy needs to run within an active ArcGIS Conda environment

texts:
I have installed ArcGIS Pro 1.3 and Anaconda 4.2.0 (in that order). I then followed the instructions from ESRI on Adding the ArcPy Path to Anaconda which made Jupyter find the arcpy module but I get error:

>>import arcpy

arcpy must be run from inside the Conda environment, which ArcGIS
uses to manage the installation of Python.
You can access this environment by launching Python from the links inside the
"ArcGIS > ArcGIS Pro" Start Menu Group:
    "Python Interactive Terminal" starts a Python session.
    "Python Command Prompt" opens a Command Prompt window initialized with Conda.

Or by running the proenv.bat file in the terminal:
     C:\> C:\Program Files\ArcGIS\Pro\bin\python\scripts\proenv.bat

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-50d283eebbb3> in <module>()
----> 1 import arcpy

C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\__init__.py in <module>()
     62 """.format(install_dir)
     63     print(msg)
---> 64     raise ImportError("arcpy needs to run within an active ArcGIS Conda environment")
     65
     66 from arcpy.geoprocessing import gp

ImportError: arcpy needs to run within an active ArcGIS Conda environment

How can I get ArcPy to work in Jupyter Notebook?


texts: How do I specify the CRS in light of the rgdal update?

texts:
My data (df) are in Africa Albers equal area projection:
  x_        y_ t_                  id
  704109. -3029239. 2015-05-12 09:28:00 A
  703680. -3029412. 2015-05-12 10:00:00 A
  703043. -3030028. 2015-05-12 11:00:00 A
  702852. -3030189. 2015-05-12 12:00:00 A
  702638. -3030319. 2015-05-12 13:00:00 A
  701831. -3031793. 2015-05-12 14:00:00 A
  701823. -3031780. 2015-05-12 15:00:00 A
  701818. -3031796. 2015-05-13 04:29:00 A
  701823. -3031796. 2015-05-13 05:01:00 A
  701815. -3031801. 2015-05-13 05:33:00 A

I used to be able to specify this with the move function. But when I run my model I get two warnings which affect my results:

Warning messages: 1: In proj4string(object) :   CRS object has
comment, which is lost in output 2: In FUN(X[[i]], ...) : CRS object
has comment, which is lost in output

This problem seems to arise from updates to rgdal and I'd like to know how to specify my projection.
Reproducible example below. Thanks.
library(tidyverse)
library(move)
library(sf)
library(amt)

track <- data.frame(df)

#' create a move object
loc <-
  move(
    x = track$x_,
    y = track$y_,
    time = track$t_,
    proj = CRS(
      "+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
    ),
    data = track,
    animal = track$id
  )

    #' Now create a dBBMM object
   dbbmm <-
  brownian.bridge.dyn(
    object = loc,
    location.error = 20,
    window.size = 7,
    margin = 3,
    raster = 10000,
    ext = 5
  )


texts: Is there a way to add custom ids to Leaflet markers derived from a GeoJSON file?

texts:
I'm fairly new to JavaScript and I'm currently trying to connect each of my Leaflet markers to a certain portion of text in another div so that when a user clicks the marker, the site will auto scroll to that relevant point in the text.
This is the code I intend to use once each of the markers have their unique ids to be manipulated:
document.getElementById("A").addEventListener("click", function() {
  window.location.href = "#B";
});

So in this case, "A" would be the unique id given to one of the markers, and #B is the id of the text I want to scroll to.
*Edit
I should clarify that the primary issue is actually giving ID's to the markers from the geoJson elements I added onto the map. The relevant code for this is here:
//convert points to circles in layer
function pointToLayer(feature, latlng, attributes){

    var attribute = attributes[0];

    // styling for the circles which are pink with white borders
    var geojsonMarkerOptions = {
        fillColor: "#800000",
        color: "white",
        weight: 1,
        opacity: 1,
        fillOpacity: 0.5,
    };

    var layer = L.circleMarker(latlng, geojsonMarkerOptions);

    //returns layer with circle markers
    return layer;

//function to fill the global attribute array variable at top of page
function processData(data){

    //variable to store the first feature in the data
    var properties = data.features[0].properties;

    //for loop to push each attribute name into array
    for (var attribute in properties){

        //indexOf will only allow attributes with population values to be pushed
        if (attribute.indexOf("Victims") > -1){
            attributes.push(attribute);

        };
    };

    return attributes; //returns attributes array to be used in callback
};

//function to retrieve the data from geojson
function getData(map){

    //load the data and calls functions
    $.getJSON("data/WarCrimes3.geojson");

};

Edit with Whole Code
var map;
var featureId =0;
var attributes = [];

function createMap(){
    // create map and set parameters

    map = L.map('mapid', {
      center: [20.20, 136.25],
      zoom: 2,
      minZoom: 4,
      maxZoom: 8,
      maxBounds: [(25.2, 115.7), (45.7, 155.54)],
    });

    getData(map);

};

function createSymbols(data, attributes){

    L.geoJson(data, {

        pointToLayer: function(feature, latlng){

            return pointToLayer(feature, latlng, attributes);
        }

    }).addTo(map);
};

function pointToLayer(feature, latlng, attributes){

    var attribute = attributes[0];

    var geojsonMarkerOptions = {
        fillColor: "#800000",
        color: "white",
        weight: 1,
        opacity: 1,
        fillOpacity: 0.5,
    };

    var layer = L.circleMarker(latlng, geojsonMarkerOptions);

    layer.featureId = featureId++;

    console.log(layer.featureId)

    return layer;

};

function processData(data){

    //variable to store the first feature in the data
    var properties = data.features[0].properties;

    //for loop to push each attribute name into array
    for (var attribute in properties){

        if (attribute.indexOf("Victims") > -1){
            attributes.push(attribute);

        };
    };

    return attributes; //returns attributes array to be used in callback
};

//function to retrieve the data from geojson
function getData(map){

    $.getJSON("data/sample.geojson", function(response){

      createSymbols(response,attributes);

    });
};

$(document).ready(createMap);

The id's are being logged, but not being assigned to the individual markers


texts: Issues with ArcGIS Web AppBuilder creating custom widget by inheriting existing widget

texts:
I am trying to extend the inbuilt Select widget in the ArcGIS Web AppBuilder for Desktop using dojo. According to the documentation by ESRI we always inherit the 'dijit/_WidgetBase' module to create the custom widget.

define(['dojo/_base/declare', 'jimu/BaseWidget'],
function(declare, BaseWidget) {
  //To create a widget, you need to derive from BaseWidget.
  return declare([BaseWidget], {
    // Custom widget code goes here
    baseClass: 'jimu-widget-mywidget'
  });
});

Since I wanted the inherit he select widget I changed the code to the following:

define(['dojo/_base/declare', 'widgets/Select/Widget'],
function(declare, Select) {
  //To create a widget, you need to derive from BaseWidget.
  return declare([Select], {
    baseClass: 'jimu-widget-selectlsa',
    startup: function() {
      this.inherited(arguments);
      console.log("tesiting")
    }
  });
});

After copying all the HTML and CSS content from the default Select widget's to my custom widget's HTML and CSS the custom Select widget works fine.

Since I have inherited the default widget's Widget.js file I am able to call functions of the default widget form my custom widget.

However, in the default widget, there is one more file named SelectableLayerItem.js which has a method name _showActions. I tired adding the following code to access the _showActions method in the SelectableLayerItem.js file:

define(['dojo/_base/declare', 'widgets/Select/Widget','widgets/Select/SelectableLayerItem'],
function(declare, Select, SelectableLayerItem) {
  //To create a widget, you need to derive from BaseWidget.
  return declare([Select, SelectableLayerItem], {
    baseClass: 'jimu-widget-selectlsa',
    startup: function() {
      this.inherited(arguments);
      console.log("tesiting")
    },
    _showActions: function(){
        //overwrite
    }
  });
});

After running the application, the application runs without any error in the console. However, my custom widget is not able to overwrite the _showActions method in the SelectableLayerItem.js file.

Instead of calling the overwritten method it always calls the default SelectableLayerItem.js's _showActions method.

How can I inherit multiple modules and access their methods independently?


texts: What is the Earth Engine server side version of JavaScript's toFixed() function for rounding off to certain decimal places?

texts: I want to round off an ee.Number() object to two decimal places. Right now I am using .toFixed() of JavaScript, but how to do that using a server-side function? I tried toFloat() but that doesn't help.


texts: Accessing "old" ("previous date") Pansharpened Landsat data in ArcGIS

texts:
A few months ago I found an excellent pansharpened Landsat dataset to use as a basemap in ArcGIS by going to the Add Data icon, selecting ArcGIS Online and then selecting add pansharpened landsat. However, when I re-open the MXD now (a few months later), I see a new pansharpened landsat basemap – most likely from a more recent image. The new image is useless for this specific area as it contains many clouds.

Does anyone know if there is any way that I can access the old (previous date) pansharpened landsat basemap?

What would I do to access a certain date?

Or is there any way that this can be downloaded somewhere else and brought into ArcGIS?

I know where to download the regular landsat imagery but I am not sure about the pansharpened imagery.


texts: I'm working with 0.5m resolutions DSM'S and DTM's. and need to calculate a canopy height model. My question is: Is ".vrt" a possible output format for a gdal_calc operation?


texts:
I have this code:
driver = ogr.GetDriverByName('KML')
data = driver.Open(kmlFile)
layer = data.GetLayer()

geom = nextFeature.GetGeometryRef()

print(geom.ExportToJson())

I am trying to export any kind of output like GeoJSON or WKT but;
{ "type": "MultiPolygon", "coordinates": [ [ [ [ 33.228451, 39.725639, 0.0 ], [ 33.18737, 39.653671, 0.0 ],
[ 33.234299, 39.644131, 0.0 ], [ 33.280369, 39.671131, 0.0 ]

the output is something like this and looks like the geometry that comes from KML file is in 3D projection or something, or at least the program thinks it is.
Here is the spatial reference if you need to know:
{
  "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json",
  "type": "GeographicCRS",
  "name": "WGS 84",
  "datum": {
    "type": "GeodeticReferenceFrame",
    "name": "World Geodetic System 1984",
    "ellipsoid": {
      "name": "WGS 84",
      "semi_major_axis": 6378137,
      "inverse_flattening": 298.257223563
    }
  },
  "coordinate_system": {
    "subtype": "ellipsoidal",
    "axis": [
      {
        "name": "Latitude",
        "abbreviation": "lat",
        "direction": "north",
        "unit": "degree"
      },
      {
        "name": "Longitude",
        "abbreviation": "lon",
        "direction": "east",
        "unit": "degree"
      }
    ]
  },
  "id": {
    "authority": "EPSG",
    "code": 4326
  }
}


texts: Find Max Value and next highest two in Column using ArcPy?

texts:
I have a piece of code where I am able to find the max value. How can i edit this code to find the second largest, and then the third? I found this code here (https://geonet.esri.com/message/237880#237880). While there is another method that identifies a way to do this, it appears like it would be much slower (bulid array, sort, and pull [-1])

myField = "MY_FIELD"
minValue = arcpy.SearchCursor(inlayer, "", "", "", myField + " A").next().getValue(myField) #Get 1st row in ascending cursor sort
maxValue = arcpy.SearchCursor(inlayer, "", "", "", myField + " D").next().getValue(myField) #Get 1st row in descending cursor sort


texts:
I've tried to create a new Mapset/Location with Python for hours but I can't figure out how it works.

The problem is that I can't use GRASS.SCRIPT because it already needs an existing Mapset and Location when I want to set it up.

So I tried to use the Shell via the Subprocess Library.

import subprocess

gisbase = os.environ['GISBASE']
cmd = gisbase + "/bin/g.mapset.exe -c mapset=testmapset location=testlocation"

p = Popen(cmd, stdin=PIPE, stdout=PIPE, stderr=PIPE)

It runs fine but I don't get any output nor a location or a mapset. Nothing.

Do you have any ideas how I could make it work?


texts: Dissolve causes 'No Shapely geometry can be created from null value' in geopandas

texts:
When I do dissolve:

With this specific shapefile happens with other similar it is working fine.
Also I checked for null values in the geometry and there weren't any.

data2 = pst[pst['PROP_TYPE'].str.startswith(('01'))]
dslvd=data2.dissolve(by='PROP_TYPE')

TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point 223630.72975008655 4405823.482638143 at 223630.72975008655 4405823.482638143
TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point 223630.72975008655 4405823.482638143 at 223630.72975008655 4405823.482638143
TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point 223630.72975008655 4405823.482638143 at 223630.72975008655 4405823.482638143

 No Shapely geometry can be created from null value

During handling of the above exception, another exception occurred:

What could be wrong?


texts: Creating a VB.net Dataset in Visual Studio 2010 from an SDE Geodatabase

texts:
I am building an ArcMap Addin that contains several Windows Forms with controls that need to be bound to a database.  I am trying to decide whether or not to use a native SQL Server Database or an SDE Geodatabase.

I would like to link the forms to the features on the map, which suggests I should use the geodatabase, but I have not found anyway to bind geodatabase tables and fields to WinForms controls at design-time.

Has anyone successfully created a Dataset (xsd) in Visual Studio from an SDE Geodatabase?


texts: PostGIS ST_Intersects is slow and doesn't seem to use spatial index

texts:
There are several threads about ST_Intersects being slow (like this one). Main advise is to use spatial index. Well I know that and created spatial index beforehand, but the query is slow nevertheless (takes about 2 hours to complete). Result of explain seems to show no sign of spatial index being used.

I want to mark which LineStrings in one table (~3M records) intersects the target Polygon in other (~30k records). The polygon is huge, but not that huge in comparison to the extent of entire dataset and I expect to catch about 100k records from the highway_only_motor table.

The query:

UPDATE highway_only_motor
   SET has_business_value = true
 WHERE ST_Intersects(geom, (SELECT geom FROM admin_boundary WHERE id = 10236))

Explain result:

QUERY PLAN                                                                                |
---------------------------------------------------------------------------------------------------|
Update on highway_only_motor  (cost=8.30..1241137.48 rows=1431564 width=1782)                      |
  InitPlan 1 (returns $0)                                                                          |
    ->  Index Scan using admin_boundary_pkey on admin_boundary  (cost=0.29..8.30 rows=1 width=6768)|
          Index Cond: (id = 10236)                                                                 |
  ->  Seq Scan on highway_only_motor  (cost=0.00..1241129.18 rows=1431564 width=1782)              |
        Filter: st_intersects(geom, $0)                                                            |

Tables:

CREATE TABLE public.highway_only_motor (
    osm_id int8 NULL,
    geom geometry NULL,
    "name" text NULL,
    highway text NULL,
    "ref" text NULL,
    oneway text NULL,
    id serial NOT NULL,
    city varchar(255) NULL,
    state varchar(255) NULL,
    country varchar(255) NULL,
    has_business_value bool NOT NULL DEFAULT false,
    CONSTRAINT highway_only_motor_pkey PRIMARY KEY (id)
);
CREATE INDEX highway_only_motor_geom_idx ON public.highway_only_motor USING gist (geom);

CREATE TABLE public.admin_boundary (
    osm_id int8 NULL,
    admin_level int4 NULL,
    boundary text NULL,
    "name" text NULL,
    place text NULL,
    population text NULL,
    geom geometry NULL,
    id serial NOT NULL,
    CONSTRAINT admin_boundary_pkey PRIMARY KEY (id)
);
CREATE INDEX admin_boundary_geom_idx ON public.admin_boundary USING gist (geom);


texts: Distance calculation with v.net.path in QGIS 3.2

texts: I have 3 layers: origin and destination point shapefiles and network line shapefiles. I want to calculate the shortest path on the network between the origin (250 features) and destinations (20 features) by using the v.net.path command. I dont understand what the CENTRE POINT LAYER and NAME OF FILE CONTAINING START AND END POINTS is (as indicated on the photo) and I couldnt find a clear description about it. Could someone please explain?


texts: How to Make a .xml file for Esri shapefile?

texts:
Mostly for my work I use PostgreSQL or .tab files (MapInfo).
However, I am familiar with conversion from one format to another using ogr2ogr.
My problem is that how can I make the metdata file or .xml file for a shapefile.
I converted it using ogr2ogr and QGIS, but none of them makes the .xml file with it.

How can I make shapefile metatdata?

I can use any software including esri for that, just need to know how.


texts: Extract values from multiple raster to .csv

texts:
I have multiple raster .tif files. I'm using r and want to extract a cell's attribute value from all rasters based on an latitude and longitude input, and saved in a single .csv file. The output should be formatted as follows:

      Lat          Long       Value*       Value**
     20.15        77.12        12            20

*denotes attribute values from raster 1.
** denotes attribute values from raster 2.

Can anyone help me in this.


texts: The attached vector file contains thousands of individual polygons (a byproduct of the underlying raster image). The result of using polygonize function in qgis. My intent is to dissolve the polygons into few continuous polygons using clustering or some other prudent algorithm. Alternately not sure if i could make changes to raster file.


texts: Arcmap 10.1 Export Map Module Missing Save Button and Name Field

texts:
I'm trying to export a map to jpg and have lost the save buttons and name field in the module. The module points to the correct folder, but I've lost scroll bars and fields. Anyone have any ideas how to fix?  Copy map to clipboard works except for the loss of one layer.

Thanks.


texts: I know you can calculate using distance bands, but can you actually illustrate them on the map in ArcGIS? I have a presentation next week, and it would be a useful illustration for me to use.


texts: Creating a marker proportionally offset along a line in QGIS 2.4 when re-projecting

texts:
One of the nice new features in QGIS 2.4 is to be able to offset markers (e.g. arrow heads) along a line as part of the styling. I am trying to display a dataset of simple lines (not polylines) with an arrowhead 90% of the way along the line. My attempted method is as follows:

Layer Properties -> Style -> Symbol -> Change...

Select "Marker Line"
Marker placement: on first vertex only
Offset along line: Map Unit
Data Defined Properties...

Offset along line = $length * 0.9

This works fine when my layer projection was the same as for the project (all EPSG:27700 - OSGB 1936 / British National Grid). However, now that I have added an OpenStreetMap basemap layer, it seems that can't be re-projected so I have to switch my project CRS to that of OpenStreetMap (EPSG:3857 WGS 84 / Pseudo Mercator). Putting in a "fiddle factor of" 1.6 seems to do the trick ($length * 0.9 * 1.6), but that's a bit of an unpleasant hack. Better solutions?


texts: Approx distance between any 2 longitudes at a given latitude

texts:
Given a SINGLE longitude coordinate in decimal degrees, how do I calculate the distance in miles or km between latitude 0 and 1 at that longitude? ie what is the formula?

For example, given 51.43567 longitude, what is the distance between 0 latitude and 1 latitude at that longitude?

Note, that I do not want any radius around the point. Just a single number in miles or km.

EDIT: I posed the wrong question above. The correct one is below. I recognise the mkenndy's answer still points to the correct formula in one way or another, however, I need the point in bold below clearly set out.

Given a SINGLE latitude coordinate in decimal degrees, how do I calculate the distance in miles or km between any longitude at that latitude? ie what is the formula?

For example, without any concept of distance *, given 51.43567 latitude, what is the distance between 1 and 2 longitiude at that latitude?

.* known distances between longitude decimal degrees are only the equator and the poles.

The answer should just be a single number in miles or km.

Please use my example and show clear workings out in your answer to give the answer in miles or km.


texts: Plotting rectangle using minimum and maximum latitude and longitude in ArcMap?

texts:
Can we plot rectangle using minimum and maximum latitude and longitude in arcmap 10.2?

Here is a sample of the data i need to plot directly using latitude and longitude in polygon shapefile.


texts:
So I am having some issues getting my USGS DEM to display in 3d in ArcScene. I downloaded my DEM (.jpg) and brought it into ArcScene. I clicked properties and base heights and floating on a custom surface, and then after accepting, the DEM is still completely flat.

What am I doing wrong?

I thought I had already put it into a projected coordinate system.

What coordinate system is it in when I download it from USGS?

How would I change it correctly?

What I am trying to do ultimately is use the 3d elevations for a flood simulation.


texts: Free / Open source GIS Education sources

texts:
I'm looking for free / open source GIS education sources for advanced GIS users.  Specifically sources that cover programming, spatial analysis, modeling, spatial databases / data management.  I use mostly ArcGIS as a federal employee, but sources that go into conceptual concepts are appreciated too.

I'm a big fan of the Penn State Masters in GIS program because they offer many of their full classes under the MIT license.  I'm also reading through the lectures from the Dublin Institute of Technology.  I'm hoping that there are other high quality resources out there that I haven't come across yet.

While I would prefer actual courses that have some degree of interactiveness (like quizes or tasks associated with lessons), books should be considered if they are high-quality, free resources.


texts: Where are the ArcFM properties such as the 'Validation Rules' for each Feature Class stored in the database?  I would like to query directly using SQL (and xpath if necessary). I use ORacle 11G GeoDatabase with ArcFM 10.1


texts:
I'd like to save WPS result (ie. image file) in Geoserver's file system.
I found out Geoserver put the result in temp/wps folder.
Can I change  it to specific folder with specific file name?


texts:
I was clipping vector layers yesterday with no problems, and now all of a sudden it doesn't want to work and keeps giving me an error stating:

  Feature (2) has invalid geometry. Please fix the geometry or change the Processing setting to the "Ignore invalid input features" option.
  Execution failed after 0.05 seconds

I even created a new project and tried to re-clip a vector layer that I did yesterday, and now it doesn't work. Also, the error message says to change the Processing setting to the "Ignore invalid input features", but I can't see any options in the clipping window.

Setting up the clip:

Fails:


texts: Selection by attribute string with wildcard as part of name in ArcGIS Desktop?

texts:
I have an attribute field in my file geodatabase (GDB) where I'm trying to select all records that contain an underbar in the string field.  The issue is that the underbar is considered a SQL wildcard on its own.

To illustrate: I need all data similar to "ID" = '42W_42_42E' or "ID" = '31A_31'

Using "ID" LIKE '%_%' returns all records, which isn't helpful.  The perfect solution would simply return a count on the number of records with the underbar in their ID attribute.


texts:
I want to convert this project for other use from OL2 to OL4 but I can't convert the code and get back a map,

http://www.trafficorders.uk/

here the OL2 code:

  layer = new OpenLayers.Layer.WMS("UKMap", "http://www.trafficorders.uk/cgi-bin/wmsmaps", {
    layers: "bgmap"
  }, {
    isBaseLayer: true,
    transitionEffect: 'resize',
    singleTile: false,
    ratio: 1.1
  });

  var areaStyMap = new OpenLayers.StyleMap({
    "default": new OpenLayers.Style({
      fill: true,
      fillOpacity: 1,
      strokeColor: "#000000",
      fillColor: "${MapColour}",
      strokeWidth: 1

    }),
    "select": new OpenLayers.Style({
      strokeColor: "#FFFFFF",
      fillOpacity: 1,
      strokeWidth: 2,
      fillColor: "#FF9900",
      label: "${Name}",
      labelOutlineWidth: 3,
      labelOutlineColor: '#ffffff',
      fontSize: "12px",
      fontFamily: "Arial, Courier New",
      fontWeight: "bold"
    })
  });

  var strat = new OpenLayers.Strategy.Fixed();
  var pmproto = new OpenLayers.Protocol.WFS({
    url: "http://www.trafficorders.uk/cgi-bin/mapserv?map=/var/www/vhosts/trafficorders.uk/httpdocs/maps/wfsareas.map",
    featureType: "la_areas",
    featureNS: "http://mapserver.gis.umn.edu/mapserver",
    featurePrefix: "ms"
  });
  arealayer = new OpenLayers.Layer.Vector("la_areas", {
    strategies: [strat],
    protocol: pmproto,
    styleMap: areaStyMap
  });

and here my new code:

var vector = new ol.layer.Vector({
  source: new ol.source.Vector({
    url: 'http://www.traffwebdev.uk/cgi-bin/mapserv?map=/var/www/vhosts/traffwebdev.uk/httpdocs/hosted/maps/wfsareas.map',
    format: new ol.format.WFS({
      featureType: ["la_areas"],
      featureNS: "http://mapserver.gis.umn.edu/mapserver",
      featurePrefix: "ms"
    })
  })
});

the featureNS seems to do not works, I don't have nothing back while in the old code I have a collection of features.

What I want to achieve is something like this but with WFS instead of geoJson

https://openlayers.org/en/latest/examples/box-selection.html


texts: Is there any statistical method to proof if the method is suitable or not

texts:
I have recently done a phenological study using threshold method from MODIS NDVI images. But I'm suffered that I don't have any situ observation (ground truth) data to do validation. Is there any statistical method to validate that this method is suitable for remote sensing based phenological study without ground truth data.

Phenology is the study of periodic plant life cycle events (onset of green-up, maturity, offset of green-up etc) and how these are influenced by seasonal and interannual variations in climate, as well as habitat factors.

Thanks in advance


texts:
I tried gdal translate to png:

gdal_translate image.bsq image.jpg -of JPEG -outsize 10% 10% -scale

and jpeg:

gdal_translate image.bsq image.png -of PNG -outsize 10% 10% -scale

jpeg looks better and file size is smaller. In general should I expect poor result with png?

Using the -r cubic and bilinear are almost the same to my eyes:


texts: Error in Zonal Analysis Stand Alone Script using PyQGIS

texts:
I want to do zonal analysis from QGIS standalone script.For this, I am using following code to run through command line

Following code is save in .bat file.

import sys
from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
 # Initialize QGIS Application
QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis", True)
app = QgsApplication([], True)
QgsApplication.initQgis()
print 'Hello QGIS1!'
# Add the path to Processing framework
sys.path.append('C:\Users\USERNAME\.qgis2\python\plugins')

# Import and initialize Processing framework
from processing.core.Processing import Processing
Processing.initialize()
import processing

#zonal Analysis

 from qgis.analysis import QgsZonalStatistics

 #specify polygon shapefile vector
 polygonLayer = QgsVectorLayer('F:/test/road_15715.shp','zonepolygons', "ogr")

# specify raster filename
 rasterFilePath = 'F:/test/abcd.tif'

# usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString  &rasterFile, const QString &attributePrefix="", int rasterBand=1)

  zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath,'vi_', 1)
  zoneStat.calculateStatistics(None)

  print 'Zonal Analysis Done!'

While using code

    zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath,'vi_', 1)
    zoneStat.calculateStatistics(None)

it is working with QGIS Python console but while using in a standalone script, zonal analysis is not done.


texts:
I am trying to use Hyperion Data for my project, unfortunately even loading the Image Collection of Hyperion exceeds my maximum 2.5 GB Data Storage.
Hyperion has a revisitrate of 200 Days, so a one Year Daterange should be sufficient to yield only one picture.
I think there is something wrong with my .filterBounds
is there some kind of way to prefilter the Collection more precisely?

Link to GEE script


texts: the export to dxf seems like it worked, but Sketchup will not import it. I don't think this is a Sketchup issue, as Sketchup imports other dxf files just fine.


texts: How to add an ID based on attributes in ArcMap 10.1

texts:
Given are:

a polygon shapefile/feature class (FC) containing municipal districts
a point shapefile/FC containing locations of, lets say, trees.

My goal is to append an unique value (1, 2, 3, ...) to the point shapefile/FC which is:

sorted from north to south (or vice versa)
and sorted/grouped by the municipal districts (the priority of sorting lies on the municipal districts)

I managed to solve 1. by generating y-coordinates for each point and than running this script which calculates and sorts unique values by means of the the y-coordinate column:

However, I do not know how to solve 2. Finally my map should look like this (the numbers in red are the ones I'd like to calculate instead):


texts: modify LocationDisplayManger class drawing coordinates

texts: I noticed that LocationDisplayManger will show a marker on ur screen the current device location and am glad for that,but the problem is that i have to do a little bit adjustment before b/c a direct conversion of GPS coordinate to my Mapview coordinate has a little gap difference on my X and Y values.so i need to subtract the X and Y coordinates by some fixed distance before LocationDisplayManger put a marker on that coordinate.can any one tell me if that is possible and if so how?


texts: Using Field Calculator to fill out a list of unit numbers

texts: I am adding address points for individual units in trailer/RV parks and I would like to be able to use the Field Calculator to assign the unit numbers.  Is it possible to tell the Field Calculator a starting number and then have it run an equation on a selection like (n)+1 where n=the previous number?  If so will it assign the numbers based on the ObjectID?


texts: Query only geometry with ArcGIS REST service in Python

texts:
I'd like to query the geometry from a Mapserver: https://map.dfg.ca.gov/arcgis/rest/services/Project_BIOS_Public

The ulimtate goal is to query all the geometry for each layer and perform an intersect against a reference shapefile. I want to determine which layers intersect with a given watershed.

This is my first time working with the ArcGIS REST service. I'm a little confused on how I can extract only the geometry and the layer name... for the purposes of this exercise I just want a list of which layers intersect with a watershed I have as a shapefile.

I've written the following code (just specific URL for example)

import requests
import json
import urllib
from collections import defaultdict

queryURL = 'https://map.dfg.ca.gov/arcgis/rest/services/Project_BIOS_Public/q_BIOS_Public_pointslines00/MapServer/0/query'
params = urllib.parse.urlencode({'f': 'geojson', 'geometry':'Envelope', 'returnGeometry': 'true', 'OutFields':'*'}).encode("utf-8")

#requests query
req = urllib.request.Request(queryURL,params)

#opens the response
response = urllib.request.urlopen(req)

#loads response as json
jsonResult = json.load(response)

jsonResult

This returns the geometry but also all of the fields, and not the layer name. How do I get just the geometry and the layer name?


texts: GDB maintenance of muliple editors without ArcSDE

texts:
I'm working with ArcGIS 10.1 SP1. I have finished 18 months of collecting data and digitization and built an impressive GDB with a lot of tables, relationship classes and 100 feature classes. Now i need to deliver it to the costumer. The GDB need to be updated twice a year by me and the customer can make some changes at tables but not at feature classes.
What are the best practices for GDB maintenance ? I thought about editor tracking and using complex relationship classes. More ideas are most welcome.
Thanks a lot, Geog


texts: Merging feature class with same name from several geodatabases into single geodatabase using Arcpy?

texts:
I have 10 geodatabases, each one has 20 feature class on average, and every one feature class has the same name. What I want is a merge from every single feature class with the same name into another geodatabase.

I made a script where first all list all the feature class from the geodatabase, but I am stuck I don't know how to realize the merge from every several geodatabases into one geodatabase.

#-*- coding:utf-8 -*-
import arcpy
from arcpy import env

arcpy.env.workspace = r'C\...\GDBESTATALES_PRUEBA' #This is my folder where I have all my geodatabases

OutputMerge = r'C\...\Resultado\GDB_Nacional\GDB_Nacional.gdb'#This is the gdb where I want to save all the merges

List_workspaces = arcpy.ListWorkspaces()
List_featureclass = arcpy.ListFeatureClasses()

for ws in List_workspaces:
    print ws
    env.workspace = ws
    fcs = arcpy.ListFeatureClasses()
    for fc in fcs:
        print fc
    #I don't know what to do after this point

How do I merge the feature class from several geodatabases into another geodatabase?


texts: Create link for open specific object in the project map in QGIS

texts:
I am working to link Excel and QGIS. Is it pssible to create a link that can be put in Excel and when cklicked open the map/project zoomed to the specific object?
Greets


texts: How to place labels under a vector layer?

texts:
I want to place the labeling under a layer. In QGIS 2.2 apparently the label appears automatically at the top. Can someone help me. (see attachment)

example:


texts: Load National Elevation Data into a WMS server?

texts:
We have an application that can pull WMS data from the National Elevation Data set (NED) and use it to perform a number of run-time calculations.  We have been pulling it from NASA JPL WMS servers here (http://onearth.jpl.nasa.gov/) but we are pulling too much data for their servers to support.  We tried using the WMS tiling extension they support, but since it is non-standard (ie. it is not TMS) we can't find any proxy or caching tools that support it (eg. TileCache, MapProxy, etc).

Does anyone know of a WMS cache that does support the NASA JPL servers?  If not, can anyone give me some guidance on where to get the raw data and how to host it up ourselves locally?  (a tutorial link would be very helpful as we don't host any other WMS data locally and would basically be starting from scratch).


texts: How to create a continuous scale with distinct/custom color and value breaks with ggplot2 map

texts:
As stated in the title, I'm trying to create a continuous scale with distinct color and value breaks within the ggplot2 package in R.

I'm using some dummy code to illustrate my problem:

library(raster)
library(ggplot2)

#load shapefile
ken <- getData("GADM", country = "KEN", level = 1)

#fortify for ggplot mapping
ken@data$id <- rownames(ken@data)
ken.f <- fortify(ken, region = "id")
ken.df <- join(ken.f, ken@data, by = "id")

#plot
ggplot() +
  geom_polygon(data = ken.df,
               aes(long, lat, group = group, fill = CCN_1)) +
  coord_equal()

In this case, we've pulled in a shapefile of Kenya's regions that houses some internal data corresponding to each region; let's pretend that the attribute CCN_1 is what we want to map. It contains continuous integers. When mapped with this basic code, the output looks like this:

Instead of having the continuous scale, I'd like to define it with specific cuts/breaks. This is achieved rather easily through QGIS, as shown here:

How can I create a scale like the one shown above with ggplot2 in R? I find the R guides for fills particularly unhelpful when I search for help with this.


texts: Anything Python related causes crash to ArcGIS for Desktop?

texts:
Problem: ArcMap 10.3 crashes to desktop anytime I use a Python related function (Toolboxes, Field Calulator, Etc)

How I got to this point: I was trying to gain some familiarity with Python programming and opted to use the guide of Learn Python The Hard Way.
This involved downloading, installing Python and setting your environment using Powershell.
I believe this caused some issues since there was already an installation of Python in effect from having ArcGIS Installed. There was already a "Python27" folder on my main HDD.

Question: How do I tell my ArcMap to look in the "Python27\ArcGIS10.3" folder for Python? I'd prefer to avoid re-installing ArcGIS entirely.

Screenshots for thought:

^ From initial research, I found that I should have some sort of Python directory entered into these fields, but I'm not sure what specifically needs to go in there.


texts: Quantifying unevenness of polygon boundaries?

texts:
I have two polygons: Polygon 1 and Polygon 2.

Using two metrics, area and perimeter length, I want to express quantitatively that Polygon 1 has a more uneven/jagged/irregular perimeter than Polygon 2.

Each polygon has the same perimeter length but each covers quite different areas. To quantify the unevenness/jaggedness/irregularity of each polygon, should the calculation be:

area/perimeter

or

perimeter/area

I thought perimeter/area, but then I found this blog post which uses area/perimeter : http://www.r-bloggers.com/measuring-the-gerrymander-with-spatstat/


texts: v.surf.idw raster file conversion into a vector with same colours style

texts: Is it possible to convert from raster to vector the raster file resulting from grass interpolation algorithm "v.surf.idw" coloured with the quantile method and maintain such colours in the vector file in QGIS.


texts: Filtering above ground LiDAR points using a catalog in lidR?

texts:
I am analyzing LiDAR data collected two separate years- 2007 and 2015. My analysis is centered around changes in vegetation during this time period. The ultimate goal is to thin the data from the 2015 data set so it is consistent with the 2007 set, normalize the data, then calculate a few metrics for each hectare (average height, crown size, etc).

The 2007 data has a significantly lower point density than the 2015 data, so my question is about thinning the 2015 data to be more consistent with the 2007 point density.

For each area, the 2007 data has been split into two files, one for ground points and one for non-ground points. The 2015 data is all contained in one file covering the same area.

I am planning to normalize the data using a DTM created from the higher density 2015 data (assuming that the terrain has not changed dramatically in the 8 years between data collection events), and so I would like to save all of the 2015 "not ground" points and then thin the point cloud.

I have been using the lidR package in R for my analysis with a catalog of the LAS files.

My idea was this:

#Get the list of LAS files to work with
files <- list.files(path = "/...2015TestGroup", full.names = TRUE,  recursive=FALSE)

#Set up the catalog properties
ctg <- catalog(files)
opt_chunk_size(ctg) <- 500
opt_chunk_buffer(ctg) <- 20
opt_output_files(ctg) <- ".../Outputs/2015nonground/{XLEFT}_{YBOTTOM}_2015"

#Select only those points not classified as "ground"
aboveground <- lasfilter(ctg, Classification != 2L)

#Get the list of 2015 "non-ground" points
files <- list.files(path = ".../2015nonground", full.names = TRUE,  recursive=FALSE)
#Set up the catalog properties
ctg <- catalog(files)
opt_output_files(ctg) <- ".../2016thinned/{XLEFT}_{YBOTTOM}_2015"

#Thin the points so they are consistent with the 2007 data
thinned <- lasfilterdecimate(ctg, homogenize(1, 1))

However, I realized that the lasfilter function doesn't work on a catalog after receiving an error message.

Is there a way to filter to keep everything but the ground points? Or is there a different approach I should take that gets me to normalized data sets of similar point densities?


texts: Understanding QgsLayoutItemLegend of PyQGIS?

texts:
As I get more into QGIS scripting, I refer to the documentation often and usually without success. I have trouble deciphering how to interpret the classes, etc.

I found this thread which is really the same question and it helps a bit but I am still confused.

Let's take the QgsLayoutItemLegend class for example from the documentation:

Now, I think I understand that in order to create a legend item and add it to my Print Layout, I would do something like this...

QgsLayoutItemLegend.create(legend)

In the QGIS python editor I am prompted to include arguments for the function, which is helpful. It looks like this:

I look in the documentation for the QgsLayoutItemLegend.create() method and see this below, which I cannot make sense of:

I see nothing about the arguments this method takes in the documentation.

How do I create a Legend using QgsLayoutItemLegend?


texts:
I previously asked this question in mapinfo-L (google groups). Since, I think that GIS & Telecommunication is broad subject, I think GIS SE is the correct platform to ask these questions:

What are the uses of GIS in telecom?
How to use GIS for telecom?
Is there any free study and test material? There are a lot of them but none are free :( .


texts: How to detect burned areas from Landsat or MODIS datasets using ILWIS?

texts:
I'm trying to find a simple tutorial on how to detect and map burned areas from/with ILWIS . I use Landsat 5/7 or MODIS for detection.
Any idea?


texts: RGB color space conversion to HSV color space in Google Earth Engine using Sentinel 2 images

texts:
I am trying to convert RGB colour space to HSV colour space in GEE. I am using function rgbToHsv() for the same. I am able to get results for Landsat 8 images but not for Sentinel 2 images. The study area is the same for both of them.

Using Landsat 8 images:

var col = landsat.filter(ee.Filter.lte('CLOUD_COVER',15)).
             filterDate('2018-01-01','2018-12-31').filterBounds(pn);
print(col);
var img = col.first();
Map.addLayer(img,
imageVisParam, 'RGB'
);
Map.addLayer(img,
imageVisParam2, 'FCC'
);
var hsv = img.select(['B4','B3','B2']).rgbToHsv();
print(hsv,'HSV');
Map.addLayer(hsv,{},'HSV');

Link to code: https://code.earthengine.google.com/aea87858ab48b64f5b3981b3d7446a83

Using Sentinel 2 images:

print(col);
var img = col.median();
Map.addLayer(img,imageVisParam,'RGB');
Map.addLayer(img,imageVisParam2,'FCC');
var hsv = img.select(['B4','B3','B2']).rgbToHsv();
print(hsv);
Map.addLayer(hsv,{},'HSV');

Link to code: https://code.earthengine.google.com/75ba3f6db7ecc497f25111547816df67


texts:
I am trying to convert 1400 tiff files to GeoJson..tif file size ranges from 20 mb to 1Gb.
gdal_polygonize.py  is taking 24 hours for Big files 5 hours for small files.
Any way to improve performance?


texts:
I am hoping to run a simple Iterator operation in ModelBuilder, but everything I attempt is not producing the right results.

I have a Geodatabase with many point files in it. I want to run and Inverse Distance Weighted interpolation on each one.

The only iterator that will allow me to connect the to as an input feature is Feature Class, however this results in only one raster being output.

Other Iterators, such as Datasets and Files, will not give me the "Input Feature" option when I try to connect it.

Please let me know if you can identify what I'm doing wrong.


texts: Losing attribute table when working remotely

texts:
Because of the pandemic my work is strictly remote. I remote into my work computer via laptop. I have two monitors connected to that laptop. Sometimes when I move the ArcMap desktop program from one screen to another my attribute table is lost beyond the visual monitors and it takes a while to get it back. Online I found that unplugging and plugging my desktop from my monitors is the only thing that works. However it is only a stroke of luck that I can click on the attribute table before it disappears again. Every time this happens I loose a good 20 minutes.
Does anyone have a better solution that doesn't involve randomness/luck or me going to the office to change my remote computer settings (we are not allowed in the office until May).


texts: How to get Network Analyst to find nearest facility on globe?

texts: Currently, the nearest facility matrix results stay within the confines of the GCS84 projection, so the route from Los Angles to Japan travels EAST not WEST. How can I change the settings to permit it to solve on a orthographic world?


texts:
I have an MXD that someone else has created. They have put some sort of lock on the extent that stops me from zooming and panning. I can't for the life of me work out how to turn it off.

Does anyone know?


texts:
I am trying to add a Leaflet control measure my map, but I get an error in the console of browser. The error is in console is:
Uncaught TypeError: L.control.measure is not a function
I copied the folder of plugin to my public html, and added the url to css and js file of plugin in my index.html
L.control.measure().addTo(map);
and this is my index.html:
<!DOCTYPE html>

    <meta charset="utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.7.1/dist/leaflet.css" />
    <link rel="stylesheet" type="text/css" href="style.css">
    <link rel="stylesheet" type="text/css" href="./lib/leaflet-measure.css">
    <title>WEB GIS</title>

    <button onclick= fullScreenView()>Full Screen</button>

    <button onclick="window.print()" >Print Map</button>


texts: How to interpret the CDIP 9band energy table

texts:
I just stumbled over
http://cdip.ucsd.edu/themes/?zoom=auto&numcolorbands=10&tz=UTC&palette=cdip_classic&ll_fmt=dm&high=6.096&r=999&un=1&pb=1&d2=p70&u2=s:100:st:1:v:de_download_table:dt:202010
http://cdip.ucsd.edu/themes/data/download/?station=100&stream=p1&sub_stream=01&public=1&start=202010&end=202010&file_type=de&download_mode=browser
and I am curious how to interpret the values in the tables?
It lists cm^2 as unit but how is that converted into Joule ?


texts:
I am following this video tutorial on how to create a leaflet legend.

here is my full JS code

            function getCountyColor(popEst){
                if(popEst>100000000){
                    return 'red';
                }else if(popEst > 50000000){
                    return 'blue';
                }else{
                    return 'green';
                }
            };
            function countriesStyle(feature){
                return {
                    fillColor:getCountyColor(feature.properties.pop_est),
                    weight: 2,
                    opacity: 1,
                    color: 'white',
                    dashArray: 3,
                    fillOpacity: 0.7
                }
            };
            var map=L.map('map').setView([40.876,-74.9999],8);
            var countries_layer=L.geoJSON(countries,{
                style: countriesStyle
            }).addTo(map);
            map.fitBounds(countries_layer.getBounds());
            var legend = L.control({position: 'botomright'});
            legend.onAdd=function(map){
                var div=L.DomUtil.create('div','legend');
                var labels=["Population greater than 100000000","Population greater than 50000000",
                            "Population equal or less than 50000000"];
                var grades = [100000001,50000001,50000000];
                div.innerHTML='
Legend</div';
                for(var i=0; i <grades.length; i++){
                    div.innerHTML+='&nbsp;&nbsp;&nbsp;'
                    +labels[i]+'
';
                }
                return div;
            }
        legend.addTo(map);

before I tried adding a legend it looked like this

now since I added the legend code no map shows up and in the console I see this error

I am following this tutorial https://www.youtube.com/watch?v=WXSIU05It4g&index=3&list=PLNCPalajQvg55_lI2bkO2mvESqxrTwyJJ


texts: Getting unique value records from a field

texts:
Building a desktop application using C#.net and ArcObjects. I'm trying to retrieve records from a shapefile using IQueryFilter. Having trouble writing the where clause in my scenario and the applicable logic.

Here is the scenario:

I have a 'Category' and "Completed" field, (both are of type string).  What I need to do is, for a given category to see whether all the values in the completed field are unique. If values are not unique I need to display those records out.

Sample data

In the above example Category B and C does not have unique values. So all these records for B and C should be displayed.

I suppose a using an IQueryFilter interface should work. But puzzled on how I can achieve this by using a correct where clause. (I just need to get where the combination of Category and Completed fields are not unique)

Please help...


texts: Control visibility of specific fields of a PBF file

texts:
I downloaded a PBF file from GeoFabrik (OpenstreetMap data) and loaded it into QGIS 3.14.
It gives me each a layer for lines, multilinestrings, multipolygons, points and another one called "other relations". Inside the multipolygon layer under the fields tab I see a differentiation of buildings, lakes and such. Since I will not be using all of the data inside each layer but only selected fields of it, how can I set the visibility for each field?
To be clear: I do not want to hide the fields in a way so I will not see them in the attributes table but I want some fields not to appear in the map canvas.
I tried this but it did not work for me as the data of all fields was still visible (in both the attributes table and the canvas)


texts: Postgis query for nearest neighbor returns distant coordinates

texts:
I have schema for a table called "places" that looks like this with about 6000 rows:

name | county| long | lat | the_geom | id

Tuthilltown | Ulster | -74.1789 | 41.6843271513 | 0101000020E61000004850FC18738B52C0974C370898D74440 |  1

the_geom was created following a tutorial from the lat/long of each places like so:

UPDATE places SET the_geom = ST_GeomFromText('POINT(' || long || ' ' || lat || ')',4326);

These are places scattered across the state of NY.

I am currently in the far north, but when I attempt to search for my nearest neighbor, it returns points in the far east e.g. Montauk.

The basic query for nearest neighbor is supposed to look like this:

SELECT name, lat, long
FROM places
ORDER BY the_geom <-> st_setsrid(st_makepoint( "my coordinates" ),4326)
LIMIT 1;

My coordinates are not stored but rather given on the fly from my phone.

To fix this, I have thus far attempted to create an index on this table, but I am not sure how to search against the index or if it would give me the proper results.

Is this an issue with psql or is it my query (more likely)?


texts: From point to Grid of polygones or polylines

texts:
Is it possible to reconstruct a grid of polygones or polylines from the centroids?

I am looking for a arcGIS tool or script that will allow me to reconstruct a grid of polygons  or polylines from their centroids.


texts: Rookie Question! I have a shapefile in QGIS 1.8.0 of a country with 20000 divisions. I need the area of each division, but this is not a given attribute. When i click on a district, the "Identify Results" panel shows a "Derived" value for the area of each district along with the feature id. Can somebody tell me how I could extract these area values to a data file? Thanks!!


texts: Populating a table that uses a different geometry from another table of a different geometry

texts:
Does anyone know how I can populate a table that uses a different geometry from another table of a different geometry? (example geometry(points) -> geometry(Linestring) or geometry(polygon)
I'm new to SQL.
This is what the source table would look like: (note that this table has type, because in column type it it is either road or district)
CREATE TABLE points(
name text,
type text,
geom geometry(POINT));

This is what I'm putting data from that table into a table that looks like this:
CREATE table road(
id serial PRIMARY KEY,
name text,
geom geometry(LINESTRING));

I'm having trouble with this as I am required to input from points into road or district and then union then afterwards in another query. so far I've tried several things, but since I'm very new to SQL I'm not entirely sure what I am actually doing This is my code:
INSERT INTO district(name, geom) SELECT (name, geom) FROM points WHERE type =  'district';
INSERT INTO district SELECT * FROM points WHERE type = 'district';
INSERT INTO district SELECT name, points.geom FROM points where points.type = 'road';

so far I've gotten an assortment of errors, mainly saying:
ERROR: column "geom" does not exist

or
ERROR duplicate key value violates unique constraint "road_pkey"
DETAIL: key(id) - (Snowflake Lane) already exists.
SQL state: 23505

currently, the only way I've gotten it to work is if I manually put in all the roads and districts but I'm not supposed to do it that way since I need to UNION combine: List the name and geometry for place and etc.
I've been stuck on this part for days and still, haven't figured it out yet...


texts: How to incorporate stdmetrics in gridmetrics function?

texts:
I am attempting to feed the lidR gridmetrics function a list of metrics. I can generate the list of metrics within a function (following this example):

myMetrics = function(z, i)
{
  metrics = list(
     zmean   = mean(z),
     imean   = mean(i),
     zwimean = sum(z*i)/sum(i), # Mean elevation weighted by intensities
     zimean  = mean(z*i),       # Mean products of z by intensity
     zsqmean = sqrt(mean(z^2))  # Quadratic mean
   )

   return(metrics)
}

And then feed the list into the gridmetrics function:

metrics = grid_metrics(lidar, myMetrics(Z, Intensity), 20)

However, I cannot find a way to incorporate individual stdmetrics into the above myMetrics function.

How, for example, would I add the zqx or zkurt metrics (from this source) into the myMetrics function, which can then be fed into the gridmetrics function?


texts:
Is it possible to enable time aware on a layer programmatically? I can add a time-enabled layer to a non-time aware data frame as shown here:

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00s300000023000000.htm

However, how can I get the layer to be time enabled in the first place (programmatically)?

I see that it is possible in arcobjects but I would rather not get into that.


texts: Effiently incorporatinig metadata (.xml) into rastered landcover data (.tiff/.toff.ovr/.tiff.vat.dbf)

texts:
So I'm using this VGIN data to understand land cover surrounding points in Virginia in QGIS. The VGIN data came with ~100 tiles (and for each there is a .tif, .tif.ovr, and .tif.vat.dbf file) and two .xml files, one of which contains the land use classifications at specific values, along with other metadata.
I am 100% self-taught on QGIS and pretty much GIS entirely. Is there a way to be using these files in tandem so that when I pull a tile into my project it will have those classifications per-categorized, or do I have to do that with each tile?
The VGIN data would be SUPER useful, if only I could wrangle it! I'd be happy to have this be a Python learning opportunity, but don't want to overcomplicate if I don't have to!
Thanks


texts: Get attribute from touching line layer (QGIS)

texts:
I have two line layers. The features of both are touching each. Now I would like to label one of the line layers with the attributes of the touching layer.

The case: In the one layer are drain's (green) and in the other one are cables (red). The cables are inside the drain (digitized with snapping to the drains). Now I need the drain-id (green labels) labeled on every cable.

I have to use shape files! So I thought to use the aggregate function of QGIS...


texts: Using IDL poly_2d.pro with known p,q matrices for larger image

texts: I determine p and q 2x2 linear image transformation matrices (2D shift, 2D scale, rotate) needed in IDL poly_2d for the center part of an image (using SolarSoft auto_align_images.pro).  I want to apply this center-determined transformation then to the full image using poly_2D on that.  How must I modify p and q for this?


texts: How can I snap lines together, but only when there is nothing in between them?

texts:
I am mapping arable field boundaries, to see where hedges might go. I already have lines drawn round the edges of every arable field. However this creates too many lines as it gives me a boundary around every arable field when in practice two fields mostly share a boundary and have only one hedge between them. I managed to get rid of all the overlapping geometries by using Union to break them into separate features and then running Delete Overlapping Geometries to remove the duplicates. However I now need to remove the arable field boundaries that are additional but not quite on top of each other. I can do this using 'Snap geometry to layer' then again Deleting Overlapping Geometries. But I want to keep double boundaries where they fall either side of a road or track, as in this case only it is possible to have a hedge either side of the road or track. I have the road and track data - Is there a way I can snap and then delete only those hedge boundaries that do not fall either side of a road or track?
Here is an image showing the hedges (in orange) with no road or path between them. The intention is to snap the geometries, and then remove one of the two overlapping features.

Here on the other hand are hedges (in orange) with roads or paths (in yellow) between them. I need to find a way to keep both of these hedges so there is still a hedge either side of the road or path.


texts: Is there a function in QGIS similar to the 'Workspace Packager' in MapInfo that allows me to save a project and all the applicable layers together in a folder. I need to send a QGIS project to a client and I'm wondering if there's an automated way of gathering all the necessary files together.


texts: Connect to postgresql using ARCSDESQLExecute

texts:
I'm trying to connect to a postgresql database using ARCSDESQLExecute from ArcPy.
When running my script as stand alone (debug from pyscripter) it works fine.
When running my script as a tool from ArcMap it works fine.

I am able to publish the tool as a GPTool and to view the test-form from the rest interface.
When running the tool from the rest interface, it fails with:

sde_conn = arcpy.ArcSDESQLExecute("localsrv", "5432", "", "myusername", "mypassword")

File "c:\program files\arcgis\server\arcpy\arcpy\arcobjects\mixins.py", line 893, in
__init__ _BaseArcObject.__init__(self, *gp_fixargs((server, instance, database, user, password), True))
File "c:\program files\arcgis\server\arcpy\arcpy\arcobjects\_base.py", line 47, in __init__ for arg in args))
RuntimeError: Object: Connection error Operation Failed Failed to execute (MyGPTool). Failed to execute

I have also tried this:

sde_conn = arcpy.ArcSDESQLExecute("localsrv", "5432", "mydatabase", "myusername", "mypassword")

When I try to use an .sde file for arcpy.ArcSDESQLExecute i get the same error.
When trying to connect to a MSSQL server, it does not crash, only connecting to postgresql server fails.

I have tried to publish it on two different arcgis servers, same behavior.
Arcgis server versions tested: 10.2 and 10.1.1


texts:
I have downloaded the following two tiff files from here:

pnv_biome.type_biome00k_c_1km_s0..0cm_2000..2017_v0.1
pnv_biome.type_biome00k_sse_1km_s0..0cm_2000..2017_v0.1

Using r.report in QGIS 3.6 I have calculated area. Results are the following:

I could not find the legend. Hence, I do not know which land cover types correspond to numbers. Can anybody help me?


texts:
I got a problem with the rasterization of a shapefile. Basically my aim is to rasterize a shapefile 'Building' by defining boundaries (Output extent from another simple polygon layer) and output raster characteristics as in the image 1:

I'd like to obtain the features of the output raster image from a reference DEM - 'Layer property' >> 'Metadata' as shown in the following picture

.

(as I would like, later on, to merge the given 'building' raster and a reference DEM)

By the time I run this command, it takes literally 2 secs and the result is a black/White rasterized image with two available values (1.79769e+308 and -1.79769e+308), all the buildings not shown and, probably, not even rasterized ( I furthermore specify I tried many times to change the symbology/style of the output raster to see if something would have appeared but nothing...).


texts:
I have created a list of dates from which I would like to get Sentinel3 images, only to find that I have multiple images per day.
How can I reduce it to 1 image per day? I intend to do use median for the merging of the photos.

var region = ee.FeatureCollection([ee.Feature(Map.getBounds(true))])

//Dates of Interest
var start = ee.Date("2018-04-01");
var finish = ee.Date("2018-05-01");

//------------------- Data bases ---------------------//
var sen3 = ee.ImageCollection('COPERNICUS/S3/OLCI')
  .filterDate(start, finish)
  .filterBounds(region)
  .select(['Oa17_radiance', 'Oa08_radiance'])
  .map(function(image){
     return image
     .clip(region)
   })

//reference  dates
var sen2_for_dates = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterDate(start, finish)
  .filterBounds(region)
  .map(function(image){
    return image
    .clipToCollection(region)
    .updateMask(image.select('QA60').not())
    .addBands(image.normalizedDifference(['B8','B4']).rename('ndvi'))
  })
  .select('ndvi')

//------------------ Unique dates ----------------------//
var sen2_dates = sen2_for_dates.map(function(image){
  return image.set('simpleDateMillis', ee.Date(ee.Date(image.date().format('YYYY-MM-dd'))))
});
var listMillisSimple = ee.List(sen2_dates.aggregate_array('simpleDateMillis'));
var uniqueDatesSimple = listMillisSimple.distinct();
var SEN2_day0 = uniqueDatesSimple

//THIS IS THE lIST OF DATES I WOULD LIKE TO GET IT's IMAGES
var SEN3_day0 = uniqueDatesSimple.map(function(date){
  return ee.Date(date).advance(0, 'day');
});

var sen3_dates = sen3.map(function(image){
  return image.set('simpleDateMillis', ee.Date(ee.Date(image.date().format('YYYY-MM-dd'))))
});

//The ImageCollection that should be containing the images from the specific dates
var allImagesInRangeS30 = sen3_dates.filter(ee.Filter.inList("simpleDateMillis", SEN3_day0))

I get these structures:
SEN3_day0 - List of 6 elements, dates (as needed).
sen3_dates - ImageCollection of 409 images (to be filtered by the list of dates in SEN3_day0).
allImagesInRangeS30 - ImageCollection of 75 images (needs to be 6, as the list).

How can I merge images from the same day?

Thank you


texts:
I have a huge Digital Elevation Model in AutoCAD Civil 3D. I used the NNI smoothing method on the DEM (with 0.2 meters raster). I want to export the surface data as a point cloud in the LAS format.
How is this possible? Maybe there is a tool for this task.
Civil 3D is able to export the DEM as LandXML. Maybe LandXML can be converted to LAS.


texts: Concatenating Attributes of Points in Same Location using ArcGIS for Desktop?

texts:
I have a single feature class that contains multiple points with the same location, but different names and text attributes, and I would like to concatenate the attributes together and merge the points down to a single point.

So a dataset like this:

+----+-----+-----+------+-----------+------------+
| ID | Lat | Lon | Name | Atribute1 | Attribute2 |
+----+-----+-----+------+-----------+------------+
|  1 |  60 |  60 | a    | a         | q          |
|  2 |  60 |  60 | b    | a         | r          |
|  3 |  61 |  61 | c    | a         | s          |
|  4 |  62 |  62 | d    | b         | t          |
|  5 |  62 |  62 | e    | b         | u          |
|  6 |  62 |  62 | f    | b         | v          |
|  7 |  63 |  63 | g    | c         | w          |
|  8 |  64 |  64 | h    | c         | x          |
|  9 |  65 |  65 | i    | c         | y          |
| 10 |  65 |  65 | j    | d         | z          |
+----+-----+-----+------+-----------+------------+

I want it to end up something like this:

+----+-----+-----+---------+-----------+------------+
| ID | Lat | Lon |  Name   | Atribute1 | Attribute2 |
+----+-----+-----+---------+-----------+------------+
|  1 |  60 |  60 | a, b    | a         | q, r       |
|  3 |  61 |  61 | c       | a         | s          |
|  4 |  62 |  62 | d, e, f | b         | t, u, v    |
|  7 |  63 |  63 | g       | c         | w          |
|  8 |  64 |  64 | h       | c         | x          |
|  9 |  65 |  65 | I, j    | c, d      | y, z       |
+----+-----+-----+---------+-----------+------------+

I have tried a bunch of different methods to do things like run summary statistics, but that only does numbers or first/last not a concatenate.

I am currently working on ideas as to how to script it out in python, but I was hoping there is a tool which I am missing that will do this already. I am also not sure on what the best route to script it is either. I was thinking about looping through the original dataset, select by location, concatenate all the fields, and put them into the first point while deleting the rest. But not sure that would be the right method.


texts: Problem with projection of Slovenian Shapefile in QGIS

texts:
I hope someone can help me for this:
I'm doing my masters thesis in Slovenia (in the field of wildlife management), for this I got a Shapefile of slovenian data by the University. As a main layer, I use CORINELandCover and World borders (TM_WORLD_BORDERS-0.3) in WGS84. When I import the slovenian layer (Data layer strictly confidential, but exactly the same problem when I import this layer: Grid 100m x 100m.zip), it absolutely doesn't match the location where it should be projected on the CORINE- and WorldBorder-layer. I'don't know the coordinate system of the slovenian shapes, but when importet into QGIS, it tells me that it's WGS84, which cannot be true, because when I zoom to the slovenian layer, the coordinates shown in the "Current map coordinate" (right, down) are definitely not WGS84 (it gives me the following extends: 342637,26823 : 656728,197123).

I have no idea what the problem is, maybe QGIS doesn't know this coordinate system, but in this case, shouldn't it report an error message regarding this?


texts: Dictionary.toArray: Unable to convert dictionary to array, key has a null value

texts:
In Google Earth Engine I'm running the Zonal Time Series Plotter.  It runs fine with some coordinates, but when I run certain other coordinates, I receive the error: "Dictionary.toArray: Unable to convert dictionary to array, key '1985' has a null value."
For example, the coordinates in this script produce the intended outcome.
While the coordinates in this script produce the error.  With nothing else in the script changing.
Could someone explain what may be producing the error?  Changing the scale of Analysis from 1000 to a smaller value will make the error go away.


texts: Join/append multiple, identical polygon layers into a single layer whilst retaining all attributes?

texts:
I have run some analysis (zonalstats) on different raster layers using the same analysis polygons. I want to concatenate all these different layers together, i.e. retain a single polygon and have a "wide" table of attributes for each one (needs to be a batch process rather than one layer at a time).

The attributes have the same name (e.g. max, min, average) so there would need to be a way to identify which polygon layer they came from when creating the merged table. I could just go into the dbf file and copy and paste into a spreadsheet as the analysis polygons are identical but this doesn't seem very elegant and I am hoping to do this with larger datasets.


texts:
I am running into problems configuring the jms clustering of Geoserver.

Here is the background:

I have 4 virtualboxes, each running a Geoserver 2.7.1 with the
jms-cluster extension. I would like to keep the embedded broker in a
matter of simpler installation. Also I would like to preserve this kind
of peer to peer network because it supports easily the add or removal of
node and because clients, do not have to know the broker network config.

However, my network admin announced me multicast will not be allowed in
production environment. So as I know a priori the ip addresses of all
the geoserver instances, I would like to change the multicast to a set
of static ip.

But, for now, I did not succeed... jms brokers do not see each others.

I have even tried with a specific broker.xml:

<?xml version="1.0" encoding="UTF-8"?>
<beans
    xmlns="http://www.springframework.org/schema/beans"
    xmlns:amq="http://activemq.apache.org/schema/core"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
    xsi:schemaLocation="http://www.springframework.org/schema/beans http://www.springframework.org/schema/beans/spring-beans-2.0.xsd
                        http://activemq.apache.org/schema/core http://activemq.apache.org/schema/core/activemq-core.xsd">
  <bean class="org.springframework.beans.factory.config.PropertyPlaceholderConfigurer"/>
  <broker useJmx="true" xmlns="http://activemq.apache.org/schema/core">

    <amq:persistenceAdapter>
      <kahaDB directory="${activemq.base}/kahadb" lockKeepAlivePeriod="0" />
    </amq:persistenceAdapter>

    <networkConnectors xmlns="http://activemq.apache.org/schema/core">
      <networkConnector uri="static:(tcp://host1:61616,tcp://host2:61616,tcp://host3:61616,tcp://localhost:61616)" />
    </networkConnectors>

    <transportConnectors>
      <transportConnector name="openwire" uri="${activemq.transportConnectors.server.uri}"
                <!-- I HAVE TRIED WITH AND WITHOUT THIS NEXT LINE -->
                     discoveryUri="static:(tcp://host1:61616,tcp://host2:61616,tcp://host3:61616,tcp://localhost:61616)"
/>
    </transportConnectors>
  </broker>
</beans>

Is there a way to configure the activemq discovery with a fixed set
of ip?

I have already checked the activemq documentation
(http://activemq.apache.org/networks-of-brokers.html) but I can have it
working so perhaps it is linked to the fact the broker is embedded.


texts: Unable to Add ArcGIS Server to ArcGIS Desktop

texts:
I am using ArcMap 10.3.1 and trying to make a connection with federated ArcGIS Server 10.6.1.
I am Getting this error message:

We were unable to connect to: http://arcgisserver:6443/arcgis. Error :
Bad syntax in request. (status code 400)

What could the problem be?


texts: Downloading vector latlong boundary data of India's Telecom circles from OpenStreetMap?

texts:
I recently came across a couple of really useful sites for downloading vector latlong boundary data for free

converts entered address into lat/long coordinates (has a batch mode):
http://stevemorse.org/jcal/latlonbatch.html?direction=forward

Tool that lets you enter a lat/long coordinate, and then get vector data for the administrative regions that point is in....has an API
http://global.mapit.mysociety.org/

The source of the data in the last link is OpenStreetMap, and these tools together have been great for getting lat/long boundary data of the districts within India's states.  However, I'm wondering if the same boundary data exists for India's telecom circles, similar to a map of US area codes.


texts:
I am using the OpenGeo suite Snapshot 2.2 I have rules in GeoServer that display layers to certain users, and only allow editing to certain layers by certain users. This works fine when using GeoExplorer, the user must login to see the layerrs / edit, but when using the SDK I get no login box, meaning the user stays anonymous which prevents editing.

I read HERE, OpenGeo SDK insturctions that in the gxp.viewer I can add authenticate? Is this what I need to add to bring the login dialogue up? if so how?


texts:
I have calculated the NDVI for my study area. However, I am having challenges to export the NDVI images.

I want to export a single image with several NDVI layers. I have tried the mosaic function but it's giving me an image with a single layer. Can someone assist me on this one?

var l8collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
.filterBounds(Border) // study area
.filterDate('2013-04-01', '2018-08-31'); // period

print(l8collection);

// Mask pixels with clouds and cloud shadows
// Surface reflectance products come with a 'pixel_qa' band based on CFMask

// create function to mask clouds, cloud shadows, snow
var maskClouds = function(image){
  // make a new single band image from the pixel qa band
  var pixel_qa = image.select('pixel_qa');
  // keep clear (0) and water (1) pixels
  return image.updateMask(pixel_qa.eq(322));
};

// use "map" to apply the function to each image in the collection
var l8masked = l8collection.map(maskClouds);

print (l8masked);

// visualize the first image in the collection, pre- and post- mask
var visParams = {bands: ['B4','B3','B2'], min: 150, max: 2000}
 Map.centerObject(Border, 9);
Map.addLayer(ee.Image(l8masked.first()).clip(Border), visParams, 'clouds masked', false)
Map.addLayer(ee.Image(l8collection.first()).clip(Border), visParams, 'original', false)

// Compute the Normalized Difference Vegetation Index (NDVI).
var getNDVI = function(img){
  return img.addBands(img.normalizedDifference(['B5','B4']).rename('NDVI'));
};

// map over image collection
var l8ndvi = l8masked.map(getNDVI);

print(l8ndvi);

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

print(ndvi);

// Create a task that you can launch from the Tasks tab.
Export.image.toDrive({
  image: ndvi.mosaic().clip(Border),
  description: 'NDVI',
  scale: 30,
  maxPixels:1e13,
  skipEmptyTiles: true
  });


texts: How does ArcGIS calculate distance of 2 points in custom projection?

texts:
I used a software called WRF to export (lat,lon) and want to draw them in ArcGIS.
but in this software, all the coordinates are caculated under this settings: the earth is as a sphere which the radius is 6370km.
Details:
-lambert projection:
cen_lat:24
cen_lon:115
true_lat1:4
true_lat2:44
false-easting:45000
false_northing:45000

How to set these numbers in ArcGIS if these coordinates are projected correctly and distance of these two points is the same as in that software?


texts:
I've just started working as a GIS technician in a land agency firm. I'm handling large datasets such as historic house sales (~500,000 rows of data per year) and executing heavy duty operations, sub-sequentially putting pressure on a laptop and PC which I've been given to work on. Currently I'm concentrating on 2D visualisation of data but anticipate 3D will come into play further down the line. I'm using QGIS running at version 2.18.5.

Laptop:

Windows 10, 64 bit
Intel(R) Core(TM) i7-6500U CPU @ 2.50GHz 2601Mhz, 2 cores, 4 logical processors
Installed Physical Memory 8.00 GB
Total Available Physical Memory 7.89 GB
Available Physical Memory 3.74 GB
Total Virtual Memory 17.8 GB
Available Virtual Memory 12.7 GB GB
AMD Radeon R7 M370 Graphics Card

PC:

Windows 7, 64 bit
Intel(R) Xeon(R) CPU E3-1225 v5 @ 3.30GHz, 3301 Mhz, 4 cores, 4 logical processors
Installed Physical Memory 8.00 GB
Total Available Physical Memory 7.89 GB
Available Physical Memory 3.85 GB
Total Virtual Memory 15.8 GB
Available Virtual Memory 11.7 GB
(On-board) Intel HD Graphics P530 Graphics Card

I've created a stress test with the following shapefile, (EPSG '27700') operating an 'inverse' distance buffer on it (conduct at your own choosing): https://www.dropbox.com/s/0e4wyp13mwbb85d/Benchmark%20Test.zip?dl=0
The shapefile shows all the planning local authorities in the UK, the boundaries of which are made up with many vertices. The test with 'shrink' the polygon features due to a negative buffer setting. The algorithm is the most demanding I've come across hence I've set it up to be my benchmark test.

Both PC and Laptop choke at 95% if distance is set at '-1000' and segments left at '5'. At 2 hours into the algorithium, there was no progress beyond 95% although it took <5 mins to get to 95%.

If the buffer is increased to '-100' the laptop and PC both manage it in good time (min:sec).

Laptop: Time: 5:50

PC:  Time: 4:30

The company is planning to get some more powerful hardware for me. The test is purely meant as a comparison and to gain some perspective.

If you are running shorter times, or have a good perception for what PC component might make QGIS run faster what is it?


texts: Adding "0" to front of all values in column using ArcGIS Field Calculator?

texts:
I have a table for each state with Zip Code values. The 5-digit Zip Codes that start with "0" have been reduced to 4-digits. I cannot figure out how to add a "0" to the front of all the values in the table. Tried a number of field-calculations but just can't seem to get it right.

Using ArcInfo 10.1


texts: Polyline disappearing intermittently in layout view

texts:
I'm using ArcGIS Desktop 10.4.1 and I'm running into an issue with a polyline disappearing. This exhibit has 4 data frames on it each with a combination of polygons and polylines.  The data frame where the issue occurs is 5.5 inch by 7 inch, with a scale of 1 inch = 820 feet.  This specific layer is symbolized into 2 fields, I'll call them 1 and 2 for simplicity.  Line 1 line type is HiddenX2, line 2 is Hidden, both width of 1.0.  They are straight lines, only 2 vertices for each.

The issue occurs when I switch from data view to layout view.  In data view, both lines show up exactly how I want.  When I switch to layout view, line 1 shows up fine, however, line 2 messes up with parts of the line disappearing. For example, one of the lines, about half is missing while another line is missing the first and last third.  Then there will be some that show up entirely.  One thing I noticed is when I switch the line type to continuous, it shows up.  Or if I keep the line type at Hidden but change the scale to 1 inch = 1000 feet, it'll show up as well.

Has anyone seen this?

Is there a fix for it?

Is it a graphics deal with the line type/scale?

What I ended up doing to cover it up is create a line in layout view to fake it.

Data View on Top, Layout View on Bottom


texts: Using Geocode Service Published to ArcGIS Server Within Silverlight Web Map

texts:
I'm trying to replicate this sample on my web map.  I have been able to build it successfully as-is within VS 2010.  However, I run into some issues when I try to get it to work with my services.  I've successfully been able to display my map service, but I’m getting snagged up on the geocoding part.  My geocode service  is published and working fine (tested it within REST). Also, a ClientAccessPolicy.XML is on our server so I don't believe there are access issues.

On the xaml side of things, I've removed all but the Address stackpanel. I think my problems are arising on the C# side.  This is the bit of code I believe is giving me troubles:

   private void FindAddressButton_Click(object sender, RoutedEventArgs e)
    {
        Locator locatorTask = new Locator("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/" +
            "Locators/ESRI_Geocode_USA/GeocodeServer");
        locatorTask.AddressToLocationsCompleted += LocatorTask_AddressToLocationsCompleted;
        locatorTask.Failed += LocatorTask_Failed;
        AddressToLocationsParameters addressParams = new AddressToLocationsParameters();
        Dictionary<string, string> address = addressParams.Address;

        if (!string.IsNullOrEmpty(State.Text))
            address.Add("Address", Address.Text);
        if (!string.IsNullOrEmpty(City.Text))
            address.Add("City", City.Text);
        if (!string.IsNullOrEmpty(State.Text))
            address.Add("State", State.Text);
        if (!string.IsNullOrEmpty(Zip.Text))
            address.Add("Zip", Zip.Text);

        if (addressParams.Address.Count == 0)
        {
          GraphicsLayer graphicsLayer = MyMap.Layers["MyGraphicsLayer"] as GraphicsLayer;
          if (graphicsLayer !=null)
            graphicsLayer.ClearGraphics();
          CandidateScrollViewer.Content = null;
          CandidatePanelGrid.Visibility = Visibility.Collapsed;
          MessageBox.Show("No address field entered.");
          return;
        }

        locatorTask.AddressToLocationsAsync(addressParams);
    }

I've replaced the GeocodeServer URL with my service. A little further down are IF statements that deal with parameters and I’m not sure how to handle. I initially tried to remove all but the first IF statements, but I get errors when I go to interact with the map (it builds successfully).  I'm not a programmer so I'm not sure how involved it would be make the manipulations necessary to get this function to work. In essence, all I need is a single input parameter, Address. I have successfully manipulated the xaml side to display a single input parameter, know I just need to manipulate C# code to make this work. Any help would be awesome.


texts: How to use an existing authentication database in source manager

texts:
I have an existing authentication db that the layers use when loading the project but when I try and add a layer using the source manager it does not seem to know about it and insists that I create a new one.  Put another way there was nothing listed under the Authentication configurations, when I filled out the basic tab it asked me if I wanted to save the configuration.

I assume I am missing something or does the source manager not use the new auth db?


texts: Source for high-res subscription imagery service

texts: I'm developing a diagramming tool using aerial imagery as a reference layer. This tool should allow the user to zoom in quite close to a roadway or intersection. Unfortunately none of the public map services I have tested supports zooming in close enough. Are there subscription map services out there that have imagery with resolution high enough to zoom in quite close to the roadway?


texts: gdal for reprojection of VSCMO VIIRS data

texts:
I downloaded h5 VIIRS file by CLASS NOOA website.

In particular I am interested in Snow Cover Binary Maps EDR and I need to reproject them in lat/lon or UTM.

I tried to use gdal warp to this aim but, although I get correct latitude, I get completely wrong longitude (for example in my AOI I got 170 degrees instead of 10 degrees E), when I open tif file using qgis or envi

Using gdalinfo I get the names of the subdatasets. In particular one of them is the subdataset:

SUBDATASET_16_NAME=HDF5:"GITCO-VSCMO_npp_d20150102_t1051499_e1057303_b16487_c20160728125230331648_noaa_ops.h5"://All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap
SUBDATASET_16_DESC=[6144x6400]//All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap (8-bit unsigned character)

consequently the commmand line I tried is the following:

gdalwarp -t_srs 'epsg:4326' HDF5:"GITCO-VSCMO_npp_d20141212_t1043482_e1049286_b16189_c20160530182604182662_noaa_ops.h5"://All_Data/VIIRS-SCD-BINARY-SNOW-MAP-EDR_All/SnowCoverBinaryMap rip.tif

I also tried to insert -s_srs 'epsg:4326' but the result is the same.

Do you have any suggestion, please?


texts:
Is there any way or shortcut to reorder QGIS layers without the mouse?

I'm having problems with a remote machine where the drag and drop mouse method does not work.


texts: Error with updating a legend's style item in pycharm

texts:
I am trying to update an item in the legend with a style item, this is the code:
#Add Legend
legend = arcpy.mapping.ListLayoutElements(mxd, "LEGEND_ELEMENT", "Legend") [0]
legendLayers = legend.listLegendItemLayers()

#Using Style Manager to edit Legend
styleItem = arcpy.mapping.ListStyleItems("C:\\Users\\Windows\\Desktop\\test\\style manager\\noise(dBA).style", "Legend Items", "noise(dBA)")
legend.update(insertLayer, styleItem)

but I keep receiving the following error message:
legend.update(insertLayer, styleItem)
AttributeError: 'LegendElement' object has no attribute 'update'

I am a complete newbie with PyCharm and coding!


texts: Can I create blended images with GDAL command line utilities?

texts:
I've written a program in Python with GDAL bindings that generates a single band georeferenced GeoTIFF as output.  From there I use the GDAL commandline tool gdaldem to colorize that image to highlight certain classes of pixels.  I'd like to also lay it on top of a basemap DEM as final output but I need to do this procedurally in my program as opposed to load both layers into a GIS tool and blend from there.

Are there any GDAL command line utilities I can use to blend two GeoTIFFs together and output a pretty png?

Here's a specific example with some random datasets:

  DEM                                  Pixels

  Pixels + DEM Blended (This is what I want to do with command line tools)

Update:

I ended up using the blend function from PIL like this:

basemap_image = Image.open('basemap.png').convert('RGBA')
overlay_image = Image.open('overlay.png').convert('RGBA')
combined = Image.blend(basemap_image, overlay_image, 0.8)
combined.save("out.png")


texts: ArcGIS server 9.3.1 Javascript API via REST : Transparency issue

texts:
I have two map services , one is tiled service returns a satellite image , and the other is dynamic displays the roads ,POIs and Parcels.

I'm trying to make the dynamic one shows the roads only above the satellite image , but i'm always getting white background with the roads and it covers the satellite image.

Here's my code :

map = new esri.Map("mapDiv");
var basemap = new esri.layers.ArcGISTiledMapServiceLayer("http://URL_TO_SATELLITE/MapServer");
map.addLayer(basemap);

var imageParameters = new esri.layers.ImageParameters();
imageParameters.format = "png24";
imageParameters.height=100;
imageParameters.width=100;
imageParameters.layerIds=[7]; // Only Roads !
imageParameters.layerOption = esri.layers.ImageParameters.LAYER_OPTION_SHOW;
imageParameters.transparent=true;

var entrmap = new esri.layers.ArcGISDynamicMapServiceLayer("http://URL_TO_DYNAMIC/MapServer",
              {"imageParameters":imageParameters});

map.addLayer(entrmap);

Thank you for your help :)


texts: SDO_Inside equivalent in PostGIS for 3D objects


texts: What is the best way to deal with inconsistent gSSURGO data across state boundaries?

texts:
I am calculating available water content (awc) using the values aws0_999 and tk0_999a from the gSSURGO Value1 table for South Dakota and Wyoming to use in a soil water balance model (SWB) to estimate recharge. The problem is that the tk0_999a (soil thickness) values from Wyoming are orders of magnitude smaller than the South Dakota values for the same soil across the border. This will definitely affect the estimated recharge from the SWB.
Are there any suggestions on the best way to deal with this problem?


texts: I have a Cartodb api for a separate map, I need to display a linestring for each time the sql is run, but not store the linestring in the database.


texts:
var bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'NDVI', 'BG', 'BR', 'SLOPE', 'SRTM', 'NDWI'];
var trainingTesting = withLayers.select(bands).sampleRegions({
  collection: new_merge,
  properties: ['landcover'],
  scale: 30
}).randomColumn('random');

print (trainingTesting);

// Partition the training
var training = trainingTesting.filter(ee.Filter.lte('random', 0.7));
print('Training n =', training.aggregate_count('.all'));

var testing = trainingTesting.filter(ee.Filter.gt('random', 0.7));
print('Testing n =', testing.aggregate_count('.all'));

Export.table.toDrive({
  collection: ee.FeatureCollection(training),
  description: 'Perlis_training',
  fileFormat: 'CSV'
});

I tried to export the whole samples before doing any splitting and the samples were all okay..

However, I want to export on the 70% of the samples which I've labeled it as training but the information on the coordinate system was not included in the export. How can I extract the coordinate of the training samples after being split?


texts: Leaflet circle has two radius properties: _mRadius, which is the radius in meter, and _radius. Is _radius in radians or what is the unit of it?