texts: Restoring a PostgreSQL/PostGIS backup into a new database with fresh extensions

texts:
I'm working on a migration from a very old server (Ubuntu 10.04/Postgres 8.x, etc) to a new instance (18.04) with up-to-date versions of everything in our web stack.

Our PostgreSQL database uses a few extensions—notably PostGIS and fuzzystrmatch—as well as several custom functions. Probably pretty typical.

I expect a pg_dump like this..

pg_dump -U username -h 127.0.0.1 dbname > dbname.backup.sql

..will include everything: extensions, custom functions, and data. But, I don't want the old extensions.

If I initially create a same-named, but empty db, like this..

createdb -U postgres -T template0 --owner=dbowner dbname

..then manually add the extensions like this..

psql -p 5432
\connect dbname;
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION ogr_fdw;
CREATE EXTENSION pgrouting;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
\q

..so that everything is recent and up-to-date..

Will a db restore from the .sql file created from a vanilla pg_dump, like above, actually overwrite (a'la CREATE OR REPLACE) any same-named functions in the new extensions?

I'm starting to assume so, and I'm re-thinking this.

tl;dr
Is the best way to do this to first invoke the --data-only option in pg_dump for our db backup, then build our new, same-named schema manually, individually add our custom functions, then import the data-only .sql file?

FWIW I'm hoping to avoid that ^ approach, because I don't want to overlook any of our custom functions by moving them manually. While many of them have a prefix signalling our organization, unfortunately, some don't. Furthermore, I'm concerned if I do a whole db restore, then try to upgrade the extensions, the upgrade process will be painful because the source db/extension versions are so out-of-date to begin with.


texts:
I'm working with data that constantly changes, and would like to write a script to automate my tasks.

My tasks are as follows:

create 3 polygons->select earthquakes within those polygons->edit earthquakes attribute table to include a field that labels earthquakes by region.

Is there a way to edit attribute tables and use the editor functions with ArcPy? I would prefer to use ModelBuilder, but I don't believe this offers the ability to edit features/tables. I'm still new to ArcPy and will post solutions if I find any.

EDIT: I'm not asking for the code I just want to know if this workflow is possible


texts:
Which professional certificates are most desired by employers seeking GIS professionals?

Are there any online courses available online that provides certification?


texts: Exporting map.graphics to image using ArcGIS API for JavaScript?

texts:
I'm developing a GIS Map Application, I'm using Esri JS API .

I have a facility to highlight some parts in the map, however I'm wondering if I can - in any way - export these graphics into an image for printing purposes.

I don't have the printing service enabled.


texts: Creating and writing GeoPackage Lookup Tables

texts:
I am writing Python using OGR to create and write a single vector layer in a GeoPackage. Some of the feature attributes need to be foreign keys to several lookup tables. My needs are completely standard SQLite capability.
Is it required or recommended that I use the same OGR driver that I use for the feature layer to create and write to the lookup tables or can I use a native Python SQLite3 connection?


texts: Defining features geometry in new created layer using PyQGIS?

texts:
I'm trying to create points in a new layer with a given distance along all lines of another layer.

Here my function so far:

from PyQt4.QtGui import *
from PyQt4.QtCore import *
from PyQt4.QtCore import QVariant
from qgis.core import *
from qgis.core import (QgsFeature, QgsGeometry,
                   QgsVectorLayer, QgsMapLayerRegistry,
                   QgsField, QGis)
from qgis.gui import *
from qgis.utils import iface
import math

def pointsAlongLine(dist):
  try:
    inputlayer = QgsMapLayerRegistry.instance().mapLayersByName('route')[0]
    print ('route-layer found')
    sum = 0
    sumdistances = []
    epsg = inputlayer.crs().postgisSrid()
    uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
    memlayer = QgsVectorLayer(uri, 'itpoint', 'memory')
    QgsMapLayerRegistry.instance().addMapLayer(memlayer)
    iface.setActiveLayer(inputlayer)
    inputlayer.selectAll()
    feat = inputlayer.selectedFeatures()
    for f in feat:
      length = f.geometry().length()
      print(length)
      for sum in range (0, (math.ceil(length) - dist), dist):
        sumdistances.append(sum)
      geompoints = [f.geometry().interpolate(distance).exportToWkt()
                       for distance in range (0, (math.ceil(length) - dist), dist)]
      memlayer.startEditing()
      prov = memlayer.dataProvider()
      feats = [ QgsFeature() for i in range(len(sumdistances)) ]
      for i, ft in enumerate(feats):
        ft.setAttributes([i])
        if i>(len(geompoints) -1):
            break
        ft.setGeometry(QgsGeometry.fromWkt(geompoints[i]))
      prov.addFeatures(feats)
      memlayer.commitChanges()
      memlayer.updateExtents()
  except IndexError:
    print ('route does not exist')

It seems that all works fine (layer is created and all points calculated and set correctly) but now I realized, that the new created points don't have their geometry defined.

How do I manage that and give them point geometries?


texts: Getting multiband raster's min and max values with Python in QGIS?

texts:
Is there anyway to have access to minimum and maximum pixel values from multi-band raster using Python in QGIS? This is QGIS and I loaded 3 different bands as rasters.

If there's a way to have access to them. How to do it?

GetDefaultHistogram gives different values for minimum and maximum. See below:

Is there a way to get the scale numbers from the histogram? For all bands?


texts: Reducing raster legend classes to the ones shown in my map in QGIS

texts:
I have a landcover raster layer. I use the styles of a QML-file. Is there an option to reduce my style classes shown in the legend to only the relevant classes shown in the map?
The QML has 44 classes, but on the map, there are only 26 classes shown. In ArcGIS, there was an option to automatically reduce the legend to the relevant classes. Is that also possible in QGIS?


texts:
I am trying to open a txt file (the textfile contains coordinates. It has a x and y column) within ArcMap 10.0. But everytime I want to open it (Add data --> Add XY Data) the software does not show me the X and Y field. Even when I try to open the txt via Data -> Add Data --> Add Data I get an error message "General function failure".
Even when I try to create a txt file like done in http://resources.arcgis.com/en/help/main/10.1%20/index.html#//005s00000010000000 I am experiencing the same problems. The same happens also if I convert the txt to csv.

Does anybody have an idea why I can't open any txt file in Arcmap?

Edit: I can open the file in ArcMap on my computer at home. Only on my works computer I am experiencing the problem that I can't open any text file. Any ideas?


texts: Compute flow accumulation only from flow direction?

texts:
I am looking to calculate flow accumulation from a flow direction raster. In ArcGIS there is Flow Accumulation (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/flow-accumulation.htm), that seems to do the job. However, I am looking something in opensource software, but all that I have found require using a DEM.

Any function in GRASS or SAGA?


texts: How to transfer table from localhost server to another server in PostgreSQL

texts:
I have a table with more than 1000000 rows (point,line,polygons) in localhost server and I want to transfer this table to another database in another server with PostgreSQL.

What is the best and fastest way to do this?


texts: Does 'r.patch' act only on "loaded" rasters?

texts: Can I run the GRASS command 'r.patch' on raster files which are not "loaded" into the current GRASS session? That is, rasters which are simply in my folder but not in a GRASS location & mapset?


texts: Adding cartodb voyager basemap to QGIS and Mapinfo

texts:
How do I add cartodb_voyager as a basemap to QGIS and Mapinfo? The map can be found here: https://leaflet-extras.github.io/leaflet-providers/preview/#filter=CartoDB.Voyager

I am using qgis v3.2 so some of the addons which would presumably add it have not been updated yet. I am using mapinfo v16 as well.


texts: Calculating field in which null values may be present?

texts:
I would like to calculate a field (short integer) using other fields (short integer) present in the table. In some cases these fields may have null values for one or more of the fields. For example, the new field name is 'total' and total is based upon adding field1 + field2 + field3. Some records (rows) have values for field1, field2, and field3, but some rows only have values for field1 and field2, then field3 is null. Do I need to go and replace all null values with a numeric value? Or, is there some workaround within the VB or python script?

I'm using ArcMap 10.1

I've used a python example from below and have been unable to get it to work. Here is the code I'm using:

#calculate new field 'summed' in feature class 'test'
inFeatures='test'
fieldName01='summed'
exSummed = "stack(!field01!+!field02!+!field03!)"
codeBlock = "def stack(*args):
    return sum(filter(None, args))"

arcpy.CalculateField_management(inFeatures, fieldName01, exSummed, "PYTHON", codeBlock)

where 'summed' is the field (short) I would like to populate and field01, field02, field03 are all short integers which may or may not have NULL's present.

I receive an error when I run this...


texts: GeoWebCache OpenLayers WMS request structure

texts:
I'm developing a GIS app that makes use of GeoServer and GeoWebCache for serving requests, and OpenLayers for frontend. WMS requests to the GeoServer are served with no issues. But when I make a request through GeoWebCache, nothing is returned to the client. I've enabled Direct WMS Integration in GeoWebCache configuration. I'm using GeoWebCache integrated in GeoServer. My request is as shown below:

var wmslayer = new OpenLayers.Layer.WMS(
   "Kenya Administrative Boundaries",
   "http://localhost:8081/geoserver/gwc/service/wms?",
   {layers: 'kenya_admin', format:'image/png', transparent: true},
   {opacity: 1.0, isBaseLayer: false, visibility: true}
);
map.addLayer(wmsLayer);

But the layer is not loaded. Any advice on what I may be doing wrong. Thanks in advance.


texts: How does one populate a column in Carto with calculated values that have rounded decimal places?

texts:
I am now able to populate a new column with a calculated percentage with this SQL:

UPDATE <table>
SET <newColumn>=subtotal/total

Could you assist me with the required SQL query to set the populated result to 2 decimal places only?

I have tried:

SELECT round(perc_gr10_11,2) FROM ctn_youth_education

which gives me a syntax error, and I have also tried:

SELECT round( CAST(perc_gr10_11 as numeric),2) FROM ctn_youth_education

which gives me (obviously) a selection of the table with the rounding applied.

However, I would like the values in the column to reflect only 2 decimal places, not just a selection of the table.  The purpose of this is for the pop-up labeling which looks weird with too many decimal places.


texts: How to improve the performance of an OpenLayers based API containing a large number of layers?

texts:
I have a GIS SDK with WMS layers and Vector layers (with information from Geoserver).

I have methods for setting a style map, but I have a performance problem.

When I set same style for 30 layers, for example, the UI gets very slow because the 30 layers have to refresh.

Is it possible to reduce the refresh time for the layers?

I have been thinking about some this, but I need some advice:

a) Is it possible to refresh vector layers without downloading the whole WFS data again? I just want to set a style with some rules.

b) Is it possible to reduce the refreshing time of WMS layers, maybe juggling with the number/size of the tiles (I am using single tile now)?

c) I tried creating a layer group in Geoserver, but I don't know how I can apply a style to each of its layers. I don't even know if it is possible to manipulate the layers individually when they are in a layer group. Is it possible?

When I change zooms or when I pan the map, it is very slow too. I think it is related to this question, so I leave it as a P.S.


texts: Select a file under two subfolders and save it in the subfolders

texts:
I want to select a file under two subfolders and save it in the subfolders I have set in the ModelBuilder in ArcGIS.

How can I do it?


texts: Is there a way to retrieve remote data through ECWP protocol in QGIS? I've read that .ecw files can be opened with QGIS, but haven't found anything concerning this particular protocol and QGIS.


texts: How add a variable from buffer selections?

texts:
New to QGIS; don't know how to use commands. Here is my question!

I use a polygon layer of forest fires in Canada (N=50 000+)
I have a point layer of centroids of specific municipalities in Canada relevant for the research (N=145)
I buffered and selected all the fires located whitin a 75km radius from these municipalities, and created a new forest fire layer (Fires_Buffer_75km; N=4905)

Now, for each of these fire polygons selected, I want to add a new variable: the NAME of the municipality from which they got buffered and selected.


texts: QGIS 3 symbols like bridleway and construction are missing. When I have updated from other versions there has never been a problem. Where are they now?


texts: Is there any possibility to draw an editable map using OSM data?

texts:
I am a newbie on openstreetmap related developments, I need to build a web application using PHP or Python frameworks, that could draw map using OSM data. It should also allows to edit attributes of nodes (ex: highway, signal lights, foot path, parking etc), edited data should not need to be updated to OSM servers.

I couldn't find any proper guide. Let me know if there any guides or tutorials to do so.

Where should I start?


texts: Arcobjects ConstructBuffers produces Error -2147220889?

texts:
I use the Esri-function ConstructBuffers in my arcobjects program. The call of this method produces a Error with code -2147220889.

Looking up the error constants it says:

  FDO_E_WORKSPACE_INVALID_KEYSETID -2147220889 Returned Keyset id is
  invalid.

what does that mean?

The geometriy list that I provide for the buffer creation is created by myself in my Software. I did some operations to make sure the geometries are correct like: Simplify, SnapToSpatialReference.
The buffer creation is executed in a loop and strangely it works with most geometry lists that I provide as input but only crashes with a few of them. I can not see why that is so. Who has an idea what could be the cause?

            IGeometryBag subcurvesBag = new GeometryBagClass();
            IGeometryCollection subcurvesCollection = (IGeometryCollection) subcurvesBag;
            GeometryBag buffersGeoBag = new GeometryBagClass();
            IGeometryCollection buffersGeoColl = (IGeometryCollection) buffersGeoBag;

            foreach (IPolyline subCurve in subCurveList)
            {
                if (subCurve.IsEmpty)
                {
                    subCurveList.Remove(subCurve);
                }
                subCurve.SpatialReference = map.SpatialReference;
                subCurve.PrepareForTopologicalOperation();//simplify
                subcurvesCollection.AddGeometry(subCurve);
            }

            IBufferConstructionProperties bufferConstProp = (IBufferConstructionProperties) bufferConstruction;
            bufferConstProp.EndOption = esriBufferConstructionEndEnum.esriBufferFlat;
            bufferConstruction.ConstructBuffers((IEnumGeometry) subcurvesCollection, bufferSize, buffersGeoColl);


texts: Split polygon by MultiLineString - shapely


texts:
I'm trying to create a geometric network in a FILE GEODATABASE, but I get:

  Error: Cannot create geometric network. The selected feature dataset
  does not contain any feature classes which can participate in a new
  geometric network.

However, when I do this on a PERSONAL GEODATABASE, it works fine.

Could this be a VERSIONING issue?

Some posts I've seen mention the VERSION toolbar, but mine is greyed out.

I'm on ARCINFO(desktop), WIN10.


texts:
I was given two DEM files covering the same area of Northern California (lets call them File_A and File_B). They're derived from LiDAR that was flown separately in the last 10 years. My goal ultimately is to compare these two datasets to see change over time. Both files have undefined vertical datums.
File_A has decent metadata and I was able to identify the vertical datum as NAVD88.  This also checks out when I compare actual elevation data points from a known good elevation file to File_A's values at that same location. Points are <<1m different. Everything is groovy so far.
Unfortunately, File_B doesn't have metadata to confirm the vertical datum.  I'm confident it isn't NAVD88 because pixels across the entire dataset are 9 to 33 meters lower than File_A with no points in File_B sitting higher than File_A (average ~23 m difference).
So my question: how would you go about identifying the original vertical datum of File_B so that it can be transformed to match File_A to enable comparisons?
If location is important, here's the center of the files: 120.6°W 40.1°N and they extent ~6 miles in all directions from that point.


texts: How to merge two (or more) features in OpenLayers?

texts:
Context: I'm building a WebApp with OpenLayers that allows users to represent on a map their commuting trips by selecting a bus stop, choosing a bus line and finally selecting where they get off the bus.

My data consists of a Shapefile with Points (bus stops) and another Shapefile with MultiLineStrings representing every allowed segment between two bus stops. They contain information about bus lines, starting and ending stops, direction, etc. All this data is served from a GeoServer and retrieved through WFS to OpenLayers.

I've made an algorithm to get all the MultiLineStrings needed to represent user's commute and it works fine. So, if user travels from stop 0 to stop 4, my trip data consists in 4 MultiLineString Features.

Problem: I'd like to merge the n MultiLineStrings of a trip into a new single one Feature that represents the whole trip from end to end (because later I need the user to rate it and stuff...)

I've looked for some feature merging method in OpenLayers' API but I didn't find anything that helps me...

So, in short, is it possible to merge two or more features in a new one?


texts: Creating Maximum-Value Composite of MODIS images using QGIS

texts:
I am trying to create a Maximum-Value composite of 4 MODIS images that are in .hdf format and I tried with 'Raster Composite' of ArcGIS but I think that's not the appropriate to do it. Can I do it using QGIS?

Data I Have: 05_2005, 06_2005, 07_2005, 08_2005 (MODIS Land Surface Temperature Data in HDF format)

I looked at this question
Creating Maximum-Value Composite using QGIS?
but this is for raster in .TIF format.


texts: Error: Exported bands must have compatible data types; found inconsistent types: Float32 and Float64

texts:
Here is my link

https://code.earthengine.google.com/ef1f791ab0223983a1177ddb7bd10699

It has 3 bands, all of them are float type.
Why this error to exporting images?
What do I need to modify?


texts: Does Avenza MAPublisher update a shapefile if it has been edited in another application (e.g. ArcMap)?

texts:
I have imported a shapefile of trails into Adobe Illustrator using the MAPublisher plugin. Then I edited that shapefile in ArcMap.

Is there a way to re-link to that file so that the Illustrator layer reflects the changes made?


texts: Satellite data to watch atmosphere and ionosphere

texts: Which bandwidths are best suited to watch atmosphere and ionosphere of  a region? And what are the sources (satellite sensors/names etc)? Any suggestion would definitely help me to narrow down my search window. I want to monitor temporal changes of atmospheric and ionospheric constituents. In case I am on a wrong channel, kindly redirect me.


texts: TerraClimate: Error generating chart: Data column(s) for axis #0 cannot be of type string

texts:
I am using "TerraClimate: Monthly Climate and Climatic Water Balance for Global Terrestrial Surfaces, University of Idaho" dataset...
When I try my code, I have this error message Error generating chart: Data column(s) for axis #0 cannot be of type string.
I'm interesting to vapor pressure deficit. Here is my code:
var dataset = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE')
.filter(ee.Filter.date('1958-01-01', '2020-12-01'))
.filterBounds(geometry);
print (dataset);

var Vapor_pressure_deficit = dataset.select('vpd');
print (Vapor_pressure_deficit);

var VPD = Vapor_pressure_deficit.map(function(img){
var date = img.get('system:time_start');
return img.multiply(0.01).set('system_time_start', date);
});

var createTS = function(img){
var date = img.get('system_time_start');
var value = img.reduceRegion(ee.Reducer.mean(), geometry).get('vpd');
var ft = ee.Feature(null, {'system:time_start': date,
                             'date': ee.Date(date).format('Y/M/d'),
                             'value': value});
  return ft;
};
var vpd_final = VPD.map(createTS);

var graph = ui.Chart.feature.byFeature(vpd_final, 'system:time_start', 'value');

print(graph.setChartType("ColumnChart")
           .setOptions({vAxis: {title: 'VPD [kPa]'},
                        hAxis: {title: 'Date'}}));
Export.table.toDrive({collection: vpd_final, selectors: 'date, value'});


texts: Calculating the percentage of overlapping areas from different layers in QGIS

texts:
I work with QGIS 3.16 Hannover and I have three polygon layers 'A', 'B', 'C' with multiple geometries in each.
In the image, layer 'A' has the bold polygons and, for each of them (meaning within each of them), I need to calculate the percentage of the area of all the blue polygons (layer 'B') that is covered by the squared polygons (layer 'C')

Following the guidelines given by @MrXsquared in this post, I built the following code in the field calculator for layer 'A', but it returns null or wrong values
area(
     intersection(
                  (buffer(collect_geometries(overlay_intersects('B', $geometry)),0)),
                  (buffer(collect_geometries(overlay_intersects('C', $geometry)),0))))
*100/area(buffer(collect_geometries(overlay_intersects('B', $geometry)),0))

Any idea of what's wrong?


texts: Publishing multiple rasters under single map service?

texts:
I have a new project in design where I need to publish a web service that will allow access for clients (like ArcGIS Server Map Services).

The only thing is, I have a big raster database (~300+ Rasters) and I need to allow users to turn each raster on/off (like in a group layer).

I experimented a little in ArcGIS 10 and stumbled upon the Group Layer which gives me a great solution to my problem:

I created a Mosaic Dataset for each raster.
I created a layer (*.lyr) for each Mosaic DS.
I gathered the layers under 1 group layer.
I published the group layer to the ArcGIS Server using an MSD file.

The service works fine and I can turn each sub layer on/off from the client. Unfortunately, I can't really test it with a large amount of rasters (I tried it with 4 medium sized rasters).

Are there any limitations for the group layer?

It seems kind of strange putting all rasters under one layer but I'm really excited by the fact that I can publish a single service for all my rasters.


texts: Extract address data / Query point in administrative boundary

texts:
I try to create a DB in postgreSQL with all the postal adress data for a specific country.

Actually I have this table full of data:

  id integer NOT NULL DEFAULT nextval('plz_id_seq'::regclass),
    country_code character varying(20),
    postal_code character varying(20),
    city_name character varying(180),
    state_name character varying(180),
    state_code character varying(20),
    community_name character varying(180),
    community_code character varying(20),
    latitude double precision,
    longitude double precision,
    CONSTRAINT plz_pkey PRIMARY KEY (id)

The next step should be, to create a new table with the streets and housenumbers and link it to the right postal code area. For this I imported the specific OSM data into postgis with osm2pgsql and the hstore option.

Where Im now stuck, is the query to find the OSM postal_code area/polygon that fits to my plz data.

I tried this query, but get nothing back:

  SELECT * FROM planet_osm_polygon
  WHERE boundary='postal_code' AND ST_CONTAINS(way, ST_PointFromText('POINT(48 7.85)', 900913));

Where ist my error?


texts: Excel co-ordinate file to polygon shapefile

texts:
So I am making a map of fishing intensity. Currently, I have an excel table with the mid-latitude and mid-longitude for each of my intensity data points. I also have the c squared value for each point on my map. As the long/lat data is the centre of each point, I want to create square polygons around this point. I am using ArcGIS Desktop 10.6.1 with an Advanced licence.

I am able to input this data as a point shapefile no problem (via add xy data), but I want to be able to show this map as a relative intensity scale using graduated colours. I can kind of do this with points, but they slightly overlap each other and it doesn't look too great (see image).

How can I either put this type of data in as a polygon from excel or convert my point shapefile to a polygon shapefile?

I have tried aggregate points, but it doesn't work as I don't know the distance between my points in their long/lat form. I remember being able to do this a few months ago in class, to get the figure below, but I've stupidly lost the polygon shapefile for it!


texts: Free website host for embedding tilemill/mapbox web map

texts:
I am hoping to design a simple interactive web map using TileMill. I am then planning on uploading that map to Mapbox.com, where I understand I will be given an embed code to put on a website. My question is: will this mapbox map seamlessly embed onto a Wordpress or Wix website? Or should I use a different free website platform? Nota bene: I have little HTML editing experience.

Thanks!


texts: How to join a table to a shapefile with non-matching IDs and names (similar strings)?

texts:
I'm having an annoying problem that I'm trying to find an automated solution for. The shorthand version is that I have a shapefile and a table of created data for regions within countries. The created data table does NOT have any sort of standardized GIDs/admin codes to match to shapefiles, and the region names are not exact matches either. Let's take a closer look; here's my dummy data frame + shapefile.

library(rgdal)

#load in shapefile
arm <- readOGR("D:/Country-Shapefiles/ARM_adm_shp", layer = "ARM_adm1")

#create dummy data frame
id <- c(100:110)
name <- c("Aragatsotn", "Ararat", "Armavir", "Gaghark'unik'", "Kotayk", "Lorri",
          "Shirak", "Syunik'", "Tavush", "Vayots' Dzor", "Yerevan City")
value <- runif(11, 0.0, 1.0)
df <- data.frame(id, name, value)

So what I have is a table with seemingly random IDs, region names, and a value to be plotted with a choropleth map. Looks like this:

> df
    id          name     value
1  100    Aragatsotn 0.6923852
2  101        Ararat 0.5762024
3  102       Armavir 0.4688358
4  103 Gaghark'unik' 0.4702253
5  104        Kotayk 0.9347992
6  105         Lorri 0.1937813
7  106        Shirak 0.5162604
8  107       Syunik' 0.4332389
9  108        Tavush 0.9889513
10 109  Vayots' Dzor 0.2182024
11 110  Yerevan City 0.5791886

Looking at the shapefile attributes of interest, we've got this:

> arm@data[c("ID_1", "NAME_1")]

       ID_1      NAME_1
    0     1  Aragatsotn
    1     2      Ararat
    2     3     Armavir
    3     4      Erevan
    4     5 Gegharkunik
    5     6      Kotayk
    6     7        Lori
    7     8      Shirak
    8     9      Syunik
    9    10      Tavush
    10   11 Vayots Dzor

Ideally, df would include some sort of matching admin IDs to join to the shapefile. Whoever created the data that I'm using did not follow these conventions, unfortunately. Alternatively, it'd be great to match the region names themselves...but as you can see, there are slight variations in each name.

Matching by hand is always a backup solution, but who wants to take the time to do that? ;) But really, even barring laziness, the project that I'm working on will be mapping dozens and dozens of different countries, so I'm looking for an automated solution that can do everything without having to do anything by hand. Is this possible? Can I somehow match these -almost- region names to the shapefiles?

Sidenote: I'm looking into grepl for partial string matches per this post, but I'm not sure if this is a potential solution because I'll need to draw from the column names rather than inputting each region name by hand.

EDIT: When I match the IDs by hand, what I've done is create a new column in my data frame and adding the exact matching terms from the shapefile. Unfortunately, because of the peculiarities of the data, the order of the names don't match up either, so this still requires some manual input. I'm hoping for some sort of completely automated solution (if it's even possible).


texts: osm2pgsql Ways importing very slowly in .pbf, much slower than .osm

texts:
I'm currently importing the planet-latest.osm.pbf file I downloaded from ftp://ftp.spline.de/pub/openstreetmap/pbf/planet-latest.osm.pbf some time last week. I'm importing it using the following command:

nohup osm2pgsql -l -U gis -c -v -r pbf --unlogged --drop -s data/planet-latest.osm.pbf &

My problem is that the ways are importing very slowly and based on current speed (0.04k/s) and the OSM stats it will take 86 days to import them all. By contrast the nodes imported at 151.3k/s. Even if I waited 86 days for the ways to import (I won't) I would then need to import the relations and sit through all the post-processing.

I recently imported another planet file in .osm format using the same command line arguments and, although it ultimately failed due to a parsing error it was moving along at a much faster speed.

From what I can see in the docs I've done everything I can to improve the import speed.

Can anyone tell me why the .pbf format import is moving so much slower over ways than the .osm format?

I'm using PostGIS 9.4 and osm2pgsql 0.88.1. Disks are on a fast SAN, the machine has a single Intel Xeon @ CPU 2.67GHz and 4GB memory. OS is Ubuntu 14.04 Server.


texts: DeleteFields without receiving error message when no fields are deleted

texts:
The code below deletes fields I deem unnecessary (those not defined in 'keepFields') as I would like it to.

def EnvCompliance_Del(Workspace):
    env.workspace = Workspace
    Fields = arcpy.ListFields('aFeatureClass')
    keepFields = ['Field_A', 'Field_B', 'Field_C')
    in_table = 'aFeatureClass' + '.dbf'
    drop_fields = [x.name for x in Fields if x.name not in keepFields]
    arcpy.DeleteField_management(in_table, drop_fields)

EnvCompliance_Del(r'C:database.gdb')

When there are no extra fields, however, I receive an error message that states:

  raise e ExecuteError: Failed to execute. Parameters are not valid.
  ERROR 000735: Drop Field: Value is required Failed to execute
  (DeleteField).

What can I add so that code prints some message stating something similar to "no unnecessary fields within 'aFeatureClass' exist" rather than the error message. My initial thought is an if-else statement, where the 'else' portion prints the message, but I am unable to successfully incorporate it.


texts: Select by location and copy to new shape and name based on attribute

texts:
I have a point features and polygon features.  I need to:

Select points that intersect the polygons
Copy the points to a new shapefile
Name this new point shapefile using an attribute of the polygon it intersects

Some of these polygons overlap each other.  If one point intersects 2 polygons, this point would need to be included in two new shapefiles.

I'm very new to Python, and therefore I'm trying to achieve this by modifying the code below.  I believe most of the basic pieces are there.

#import libraries
import arcpy, os, sys

#set current workspace
arcpy.env.workspace = r"C:\data"
arcpy.env.overwriteOutput = True

#assign variables
rawSurvey = "rawSurvey" #rawSurvey is the point file from the survey company
surface= "surface"       # surface is the approach surface polygon
output = r"C:\Users\data\airspace_Working2.gdb"

#create cursor to cycle through each "Type" of surface
cursor = arcpy.da.UpdateCursor ('surface', "Type")
for row in cursor:
    # selects the survey features intersecting (under) the surfaces
    # Make a layer and select rawSurvey points which overlap each polygon
    arcpy.MakeFeatureLayer_management ('rawSurvey', 'Tes_lyr')
    arcpy.SelectLayerByLocation_management ('Tes_lyr', 'intersect', 'surface')

# If features matched criteria write them to new feature classes appSurvey_"type" based on which type of surface they ovelrap
matchcount = int(arcpy.GetCount_management ('Tes_lyr').getOutput(0))
if matchcount == 0:
    print ('no features matched spatial and attribute criteria')
else:
    arcpy.CopyFeatures_management ('Tes_Lyr', 'appSurvey')
#print the number of features matched
print (str(matchcount) + ' points were under the surfaces and written as appSurvey')


texts:
I have a div element

And I have a map

 var map = new ol.Map({
    layers: [map_layer],
    target: 'map',
    view: Map_1_view
});

I try this to set target

 $("#map").data("my_map",  map);

but it didn't work. How can I solve this?


texts: OpenLayers: Combine multiple OGC filters

texts:
I have two OGC filters client side:

<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <ogc:PropertyIsEqualTo matchCase="true">
        <ogc:PropertyName>code</ogc:PropertyName>
        <ogc:Literal>secret</ogc:Literal>
    </ogc:PropertyIsEqualTo>
</ogc:Filter>

<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    <ogc:Intersects>
        <ogc:PropertyName>geom</ogc:PropertyName>
        <gml:Box xmlns:gml="http://www.opengis.net/gml">
            <gml:coordinates decimal="." cs="," ts=" ">1298543.0367787,7107400.531374 1299312.4070614,7203464.9016567</gml:coordinates>
        </gml:Box>
    </ogc:Intersects>
</ogc:Filter>

I have both filters represented as strings and I want to combine them into one filter. I can see in the documentation that OpenLayers can create a Filter Object from a string, but it doesn't seem to work for me.

Am I doing something wrong and/or how should I do it?


texts: Correct Projection for calculating centroids of a group of points?

texts:
I have a long list of points covering all of North America for which I would like to calculate the centroid (by grouping points based on the "Name" column).

My approach is to take the WGS84 co-ordinates and project them into North America Equidistant Conic; and then use gCentroid from rgeos to calculate the centres, like so:

#Using: http://www.spatialreference.org/ref/esri/102010/ we get the Proj4js format
no_am_eq_co <- "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
wgs84 <- "+proj=longlat +datum=WGS84"

# FROM: Coordinates are geographic latitude/longitudes
coordinates(in_data) <- c("lon", "lat")
proj4string(in_data) <- CRS(wgs84)

# TO: Project into North America Equidistant Conic
df <- spTransform(in_data, CRS(no_am_eq_co))

# Get centroids
ctrs <- lapply(unique(df$Name), function(x) gCentroid(df[df$Name==x,]))
ctrsout <- setNames( ctrs , unique(df$Name ) )

# Create data frame
df <- do.call(rbind, lapply(ctrsout, data.frame, stringsAsFactors=FALSE))
coordinates(df) <- c("x", "y")
proj4string(df) <- CRS(no_am_eq_co)
df <- as.data.frame(spTransform(df, CRS(wgs84)))
names(df) <- c("longitude", "latitude")

However, I wonder if using for example the lambert projection would be better (i.e. more accurate centroids)?

It is my understanding that the equidistant conic minimises distance distortions(which is what I think I need to get centroids), compared to the lambert conformal conic which minimises shape distortions.


texts: What is the process to calculate Helmert parameters?

texts:
I have some double points (the same points in 2 different datums).
I'd like to calculate the Helmert parameters but I don't understand how.

NOTE: I don't need software.  I need an algorithm in order to write a small Delphi application.


texts: Import SQL Server Spatial Data into PostGIS

texts:
Is there a free tool to import SQL Server spatial data into PostGIS?

Specifically I want to be able to provide connection strings to my SQL Server and PostgreSQL instances, specify the table(s) I want to import, and click a button to do the import.


texts:
I'm trying to post a new SLD style to geoserver and associate it to an existing layer, but keep getting errors. I'm doing the following:

curl -v -u XXXX:xxxxx -X POST "http://localhost/geoserver/rest/layers/landusemap_t11/styles" -H "content-type: application/vnd.ogc.sld+xml" -H "Content-Type: landusemap/collection/t11/defaultstyle/t11.sld" -d "<?xml version=\"1.0\" encoding=\"UTF-8\"?><StyleInfoPost>\t<name>landusemap_t11</name>\t<filename>t11.sld</filename></StyleInfoPost>"

I'm getting:

  Could not determine format. Try setting the Content-type header.

What am I doing wrong?


texts: Measuring distance between points in ArcMap - what coordinate system to use?

texts:
This is a basic -newbie - question. I have some handwritten data points taken off a GPS - so I am assuming they are in WGS 84, and I want to load them in ArcMap 9.3, via excel, and measure the distance between them (they are transects). I would like to know how best to maintain accuracy, and what format the data the coordinates should be when I load them into ArcMap.

The coordinates look like: 39. 51 550, 148.16 841 but the way they have been written they could be 39. 51 .50, 148.16.841 or 39 51.550, 148 16.841. Do these look like WGS 84 coordinates and what is the best way to input them into excel? Do I have to set or reset the projection in ArcMap before I take the measurements?

As you can see I am not confident with the basics of the coordinate system and projection system so any help and tips would be appreciated.


texts: Problem with service area gap with Network Analyst

texts:
I am working on running Network Analyst to create a Service Area from a single location to all around it. Esri's network dataset was missing a couple crucial street lines, and so I uploading their street shapefile and added in three streets. Upon running the Service Area, I get a couple expected gaps, due to zero streets in that area (parks and forests etc). However, in the image below, I have a neighborhood where I have added in two streets that do not appear to be getting picked up, and thus creating a gap in this area.

These added lines have been calculated for Length_Mi, Speed Limit, and Minutes for time for travel (all the things required to create a network dataset)

I have checked that these lines are connected to endpoints (zoomed all the way in), but still no luck.

What else could I be missing?


texts: error: 'Cannot read property 'divideBy' of undefined', L.Proj with esri.leaflet

texts:
I am trying my first attempt at adding an (in-house published) ArcGIS basemap to a Leaflet based web application, using Proj4Leaflet and esri-leaflet.

Here is my js:

var crs = new L.Proj.CRS('EPSG: 2263', '+proj=lcc  +lat_1=40 +lat_2=41 +lat_0=40 +lon_0=-74 +x_0=984250 +y_0=0 +datum=NAD83 +units=us-ft +no_defs', {
 origin: [-120039300, 145506800],
 resolutions: [
  295.138888888889,
  217.013888888889,
  108.506944444444,
  55.5555555555556,
  27.7777777777778,
  13.8888888888889,
  6.94444444444444,
  3.47222222222222,
  1.73611111111111,
  0.8671875
]
});

var map = L.map('map',{
  crs: crs
  }).setView([40.721, -73.994], 5);

L.esri.tiledMapLayer('http://devserver/ArcGIS/rest/services/basemap/MapServer',{
    maxZoom: 9,
    minZoom: 0

}).addTo(map);

All of the projections parameters came from the properties of the source mxd of the tiled map service.  I tried integer-only values only because all the examples I found used integer values as well.
I was not able to find definitive documentation on how to define the properties for different projections in an L.Proj.CRDS object, so I tried following examples I found online.

This code results in a blank page, and I receive the error 'Cannot read property 'divideBy' of undefined'.  I am assuming something is missing or wrong in my reprojection parameters, but I haven't had any luck yet.


texts: How to convert GML Surface to SHP Using OGR?

texts:
I have a very large .gml file with multiple layers in it.  I would like to export one of these layers to .shp.  It is a polygon layer (of lakes).  However, I am getting the error:

  ERROR 1: Unrecognised geometry type <Surface>.

The ogr command I was using:
ogr2ogr -f "ESRI Shapefile" output.shp input.gml layername

Any thoughts?  Is it possible the polygon layer is somehow incorrectly defined as a surface in the GML?  Any way I can force it to be a polygon?  Or is it something else entirely?  I don't often work with .gml, so may be missing something obvious...

EDIT 20101202:
Here is a snip of (what appears to be) the relevant info from an example feature in the example file linked below:

<gml:featureMember>
    <gml:surfaceProperty>
        <gml:Surface srsName="EPSG:4140" srsDimension="2">
            <gml:patches>
                <gml:PolygonPatch>
                    <gml:exterior>
                        <gml:LinearRing>
                           <gml:posList>
                               52.2023177 -77.596532 52.2024007 -77.5963559 52.2028141 -77.5962899 52.2028362 -77.5964181 52.2026119 -77.5967049 52.2024106 -77.5967968 52.2023273 -77.5966966 52.2023177 -77.596532
                           </gml:posList>
                        </gml:LinearRing>
                    </gml:exterior>
             </gml:PolygonPatch>
        </gml:patches>
    </gml:Surface>
</gml:surfaceProperty>

So, the file is identifying it as a 'surface' as well as a 'PolygonPatch'?  Perhaps the question should be:  what is a 'surface' in a .gml file?


texts: Preserving raster .tif.aux.xml with gdalwarp

texts:
I created a raster in R and wrote it to a .tif file which created an associated .tif.aux.xml file with the raster attribute table. Later I wanted to call gdalwarp to downsample and/or crop it. However, it does not create a new .tif.aux.xml with the attribute table, so the raster's attributes are lost. How can I preserve the attributes?
reproducible example
R code to generate raster files
library(raster)

m <- matrix(sample(1:5, size=100, replace=TRUE), ncol = 10)
r <- raster(m)
r <- ratify(r)
levels(r) <- c('city','mountain','river','lake','ocean')
colortable(r) <- c('blue','green','red','purple','yellow')

r <- raster(nrow=10, ncol=10)
values(r) = 1
r[51:100] = 2
r[3:6, 1:5] = 3
r <- ratify(r)

rat <- levels(r)[[1]]
rat$landcover <- c("Pine", "Oak", "Meadow")
rat$code <- c(12,25,30)
levels(r) <- rat
r

colortable(r) <- c('black', 'darkgreen','darkred','goldenrod')
plot(r)

# Write the raster
r_tif <- writeRaster(r, 'testraster.tif')

This creates two files, testraster.tif and testraster.tif.aux.xml
Now I execute this system command:
gdalwarp -tr 1 1 testraster.tif testrasterwarped.tif

At this point if I go back to R and call:
raster('testrasterwarped.tif')

I get no attribute table:
class      : RasterLayer
dimensions : 171, 342, 58482  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : -180, 162, -81, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : C:/Users/qread/Documents/temp/testrasterwarped.tif
names      : testrasterwarped
values     : 0, 255  (min, max)

But if I manually copy testraster.tif.aux.xml to testrasterwarped.tif.aux.xml I get:
class      : RasterLayer
dimensions : 171, 342, 58482  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : -180, 162, -81, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : C:/Users/qread/Documents/temp/testrasterwarped.tif
names      : testrasterwarped
values     : 0, 255  (min, max)
attributes :
 ID category
  0
  1     Pine
  2      Oak
  3   Meadow

How to preserve the attribute table to begin with?


texts: Converting elevation from EU-DEM (ETRS89) to WGS84

texts:
I'm having a hard time trying to use elevation data from EU-DEM.

This data is in ETRS89 (EPSG:4258), according to their site.

Now, when I extract the elevation for a shore, I get 0 meters of altitude from those GeoTIFF files. This must be elevation over sea level, as it is a beach. A GPS using EPSG:4326, however, gives around 50 meters for the same spot (which matches the Geoid height for that location).

I need to be able to compare elevation values from GPS devices with data from this DEM files, I thought this would be maybe because the geoids are  different (one is GRS80 and the other is WGS84). It turns out they are practically equivalent, and I tried converting between the two and the z value is not changed:

echo "-8.790094 42.198734 0.0" | gdaltransform -s_srs EPSG:4258 -t_srs EPSG:4326
-8.790094 42.198734 0

So my questions are:

What kind of elevation is provided by those EU-DEM files? It doesn't seem to be height over the ellipsoid, but rather height over mean sea level. I can't believe this is the case, as it would difficult using GPS 3D positions along with elevation data, and it seemed to me more straighforward that elevation survey satellites referenced elevation data to a standard ellipsoid.
What is the correct way to convert between elevation values in these GeoTIFF files to WGS84 elevation coordinates?


texts: Converting CSV with LineString into shapefile using ArcGIS Desktop?


texts:
In the definition of surface plots in my book, (An introduction to data mining, chapter 3) it is said that:

  Surface plots are often used to describe mathematical functions or physical surfaces that vary in a relatively smooth manner.

I don't understand what exactly is meant by the above sentence. What special property do mathematical functions and physical surfaces have that is appropriate for the case of surface plot?


texts: Why do raster values appear in display but not in attribute table?

texts:
I received a raster dataset from a colleague consisting of integer values between 1 to 4. The raster displays the correct range of values, but the attribute value contains only a count for values of '0', and no information about the other values that I know are also present based on the display.

This is a problem because I need to perform a series of manipulations and analysis on this data, but these have proved impossible because of this disconnect between how the data is displayed and how it is recorded in the table.

The following images hopefully show clearly the issue I am encountering. The display issues are not my main problem, but seem to symptoms of something happening with this file that I do not understand.

Symbology:

Attribute table:

Display at large spatial scale:

Display after zooming in:


texts: Plotting large shapefiles with matplotlib

texts:
I'm writing a little GIS tool in Python and I want to display shapefiles. I've successfully gotten matplotlib working using this answer, but it is very slow and heavy to load any shapefile that happens to be over 5mb. The code follows:

if len(shpFilePath) > 0:
            try:
                shp = shapefile.Reader(shpFilePath)
                plt.figure()
                try:
                    for shape in shp.shapeRecords():
                        x = [i[0] for i in shape.shape.points[:]]
                        y = [i[1] for i in shape.shape.points[:]]
                        plt.plot(x, y)
                    plt.show(1)
                except AssertionError:
                    self.error_popup('Error', 'Shapefile does not contain points.',
                                     'Check feature type and try again')

            except shapefile.ShapefileException:
                self.error_popup('Error', 'Shapefile not found.', 'Ensure that the path is correct and try again.')

Is there anything I can do to improve the speed and the memory that this consumes? I'm sure the problem lies in the creation of the x and y lists.

I loaded a 10mb shapefile yesterday and it ended up using nearly 1gb of RAM just to display it. I've been looking at other libraries but it seems like matplotlib is the most frequently recommended one.


texts: Getting query for OSMquery-Tool using dynamic codeblock

texts:
The OSMQuery tool (https://github.com/riccardoklinger/OSMquery/releases/tag/1.0.0) makes a request to Overpass API and returns a layer. The in_query makes the request to the Overpass API and must be written in Overpass QL. It is to be passed like a code block. So with """code""". I like to have the queries generated by algorithm in a later step. For this I need to be able to exchange elements of the query.
Does anyone have an idea how I can achieve this using ArcPy with ArcGIS Pro?
As the script stands it does not work. For it to work #2 must be selected. But this option is static and not dynamic.
import arcpy

arcpy.ImportToolbox(r"C:\Arbeitsordner\OSMquery-1.0.0\OSMquery-1.0.0\OSMQuery.pyt")

arcpy.env.workspace = "C:/Users/user/Documents/ArcGIS/Projects/BraunDaten/Default.gdb"
#tabelle = "Straßenverzeichnis_osmST_Braunschweig"

#Feld "shape_length" hinzufügen
#arcpy.management.AddField(tabelle, "shape_length", "DOUBLE")

ags_name = "Braunschweig"
str_name = "Celler Heerstraße"
plz = "38114"

semi= """"""
ags_f = f"(area[name={ags_name}];) ->.a;"
plz_f = f"(area[postal_code={plz}];) ->.b;"
str_f = f"way(area.a)(area.b)[name={str_name}];"
in_query = semi + ags_f + plz_f + str_f + semi
# in_query = "\'" + "\'" + "\'" + ags_f + plz_f + str_f + "\'" + "\'" +"\'"

#1: in_query = f"(area[name={ags_name}];) ->.a;(area[postal_code={plz}];) ->.b;way(area.a)(area.b)[name={str_name}];"
#2: in_query = """(area[name="Braunschweig"];) ->.a;(area[postal_code="38114"];) ->.b;way(area.a)(area.b)[name="Celler Heerstraße"];"""
#3: in_query = f"""(area[name={ags_name}];) ->.a;(area[postal_code={plz}];) ->.b;way(area.a)(area.b)[name={str_name}];"""

arcpy.OSMQueryToolbox.GetOSMDataExpert(in_query)


texts: Adding a new link to an existing road shapefile in QGIS

texts: I have an existing road file of Auckland. I wish to add a new road to the existing shapefile so I can calculate a revised distance isochrone. However, when I add the new link, it is ignored by the existing network.


texts: Using GeoTools to open Esri File Geodatabase (FileGDB)?

texts: How can I open an Esri File Geodatabase (FileGDB) in GeoTools?


texts: Creating multiple imageCollections in loop in Google Earth Engine

texts:
I am writing a large loop and for each iteration, I want to use / manipulate a different imageCollection in order to use as little memory as possible.

However, it seems that this is impossible due to how imageCollection() interprets what's a string and what's not. Does anyone know why or, if it's not impossible, what GEE specifications are?

Here is some code:

// Make a list of Features.
var features = [
  ee.Feature(ee.Geometry.Point(97.58, 24.36), {name: 'sample1'}),
  ee.Feature(ee.Geometry.Point(99.45, 20.26), {name: 'sample2'}),
  ee.Feature(ee.Geometry.Point(93.24, 21.42), {name: 'sample3'}),
  ee.Feature(ee.Geometry.Point(96.13, 24.12), {name: 'sample4'}),
  ee.Feature(ee.Geometry.Point(96.88, 22.48), {name: 'sample5'})
];

// Create a FeatureCollection from the list and print it.
var points = ee.FeatureCollection(features);

// list containing 3 satellite sources
var sat =
  ee.List(["LANDSAT/LC08/C01/T1_TOA",
    "MODIS/006/MOD09GQ",
  "MODIS/006/MOD09GA",
  ]);

//this doesn't work
print(ee.ImageCollection(ee.String(sat.get(0))));

//nor this
var test = sat.map(ee.ImageCollection().
                  .filterDate('2015-06-23', '2020-01-31');
print('test', test);

//nor this
var satindex = ee.List.sequence(0,2);

satindex.size().evaluate(function(nlist){
  for(var j = 0; j<nlist; j++){
    var satname = ee.String(sat.get(j))
    print(satname)

    var satimgs = ee.ImageCollection(satname)
                  .filterDate('2015-06-23', '2020-01-31')
    print(satimgs)
  }
});

//but this one does work
var mytest = function(satellite){
  var imgs = ee.ImageCollection(satellite)
    .filterDate('2015-06-23', '2020-01-31')
    .filterBounds(points)
  return imgs
};
print('mytest', mytest("LANDSAT/LC08/C01/T1_TOA"));


texts: Why does @data@values slot for RasterLayer only contain logical(0) and not the actual values?

texts:
Trying to get to the bottom of why, when I read in a raster of NDVI, the @data@values slot does not contain the actual values until I set them manually. For example:

    NDVI <- raster("./filename.tif", crs="+proj=longlat +datum=WGS84")
    NDVI@data@values
            ## returns: logical(0)

This did not happen with other rasters I had loaded in by the same method, so I'm confused. I wish I could be more specific but I do not remember doing anything differently before. It's easy enough to get the values manually, using:

    NDVI1@data@values <- getValues(NDVI19east)

But it's still a pain to have to do so for every file. So, two-part question:

Why did this happen in the first place? I understand that it could have something to do with how the raster file is being stored (i.e. whether it's in memory or not) but I can't really understand how that changes the methods I should use for accessing the data...
Is there a way to automate this process (perhaps using a method similar to lapply) for reading files as RasterLayers and accessing values for those files? My current project involves reading 6-10 files at a time for NDVI, Rainfall and other environmental variables, in order to combine them and perform some weighted overlays. It would be helpful to automate the process of importing the data.


texts: Find data: Number of households by income level per U.S. state?

texts:
Do you by any chance know an online source, where I can find the number of households by income level per U.S. state?

I have been going through all of the U.S. census data, but I only can find total US number of households by income level or median household by state.


texts: How to convert a raster from Integer format to Float

texts:
I have a classification image in integer format but need it in a float for weightings before running a least cost path analysis.

I have tried the raster calculator:

float("raster_1")

But this did not complete the conversion.

How should I complete this task?


texts: Code to return false if two preconditions are false in Modelbuilder

texts:
I'm trying to use Calculate Value in Modelbuilder's Arcmap to run a process if one of two booleans preconditions is true. So if both of precondintion are false the process shouldn't run, otherwise it should run.
First I used the following python code:

Expression: pergunta("%parameter1%","%parameter2%")
Code Block: def pergunta(par1,par2):
                import arcpy
                if par1 == "false" and par2=="false":
                               return "false"

This returned false every time, even when it was supposed to be true.
So I changed the code block for the following:

def pergunta(par1,par2):
                import arcpy
                if par1 == "false" and par2=="false":
                               return "false"
                else:
                               return "true"

In this case it returned true every time, even when it supposed to be false.

My last attempt for the code block was the following:

def pergunta(par1,par2):
                import arcpy
                if par1 == "false" and par2=="false":
                               return "false"
                elif par1 == "true" and par2=="false":
                               return "true"
                elif par1 == "false" and par2=="true":
                               return "true"
                elif par1 == "true" and par2=="true":
                               return "true"

Once again, it returned false every time, even when it was supposed to be true.

It seems like no matter what I do, it is not considering if the preconditions are false or true.
What am I doing wrong?


texts: How to spatially join attributes of a polygon closest/nearest to a point in QGIS?

texts: this seems like it should be a simple task but so far I've not had any luck. I have a points shapefile (showing where manholes are on public roads), and a polygon shapefile showing where building foot prints are. I want to take the attributes of the closest building footprint polygon to a manhole point - using a nearest distance algorithm. Is there any simple way to achieve this in QGIS?


texts: How to disable advanced label expressions for all layers using ArcObjects?

texts:
I have some MXDs with a large number of labeled layers. Recently I discovered that the Tiled Labels to Annotation tool didn't work with advanced Python label expressions (with FindLabel functions), so I had to simplify the expressions. But when changing the label expressions using Arcpy, the "Advanced" label expression checkbox remained checked.

I know I can't uncheck the checkboxes using Arcpy, but how can I disable the advanced expressions with ArcObjects? I'm completely new to ArcObjects, so can anyone give me some pointers how to start with this?


texts:
I am trying do some analysis of two feature layers in ArcGIS online. Basically, I want to extract and show only the points that are within the boundaries of the second layer. Trying to perform "Select by location" type of analysis but in AGOL. Filtering through the definition query won't work since these are two different feature layers.

Is there a way to do that?


texts: openlayers - export map to canvas problem

texts:
I am trying to make a export control to export the map as image file.
I found this example:
http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/examples/exportMapCanvas.html
so i added the control to the map and the function

function exportMap() {
    var canvas = OpenLayers.Util.getElement("exportedImage");
    exportMapControl.trigger(canvas);
}

I also added the ExportMap.js file:

http://dev.openlayers.org/sandbox/camptocamp/canvas/openlayers/lib/OpenLayers/Control/ExportMap.js

but, i am getting this error, when i am trying to run exportMap():

Uncaught TypeError: Failed to execute 'drawImage' on 'CanvasRenderingContext2D': No function was found that matched the signature provided. ExportMap.js:128

What that mean? how can i fix it?


texts: Input a geometry spline through a negative buffer

texts:
Update My question is a duplicate of Extracting centerline of a Complex-Polygon in PostGIS/Python (though if anyone has a programmatic solution - that is, not ArcGIS but rather Py, C, Java, etc. - I would be interested).

Let's say I have shape such as the below:

From this shape, I want to perform a negative buffer that produces the following:

The rough idea was to perform a negative buffer that exceeded the diameter of the geometry at any given point. This does not work, though, and leads to an empty geometry.

How might I perform a geometric operation on a shape so as to impute its centermost points? Is some form of negative buffering a reasonable approach? Or, is there a more well-word path to achieving the spline of such geometries?


texts: Converting Gtiff to ASCII-Grid using QGIS gdal_translate?

texts:
Gdal_translate implemented in QGIS is often recommended for converting to Ascii-Grid. My problem is that every output which is generated from my geotiff-raster layer is binary instead of plain-text in my Editor.

gdal_translate -of GTiff D:/.../lu_krs_attr.tif D:/.../grid/test_grid.asc

How can I generate an Ascii-Grid-File in plain-text with header etc?

The GTiff was originally created in an Arc*-Software, it has 1 channel and a raster attribute table implemented (which is not assigned in QGIS and which I don't want to convert either).


texts:
Purpose:To take an input raster image and clip its shape to a parcel block. In later stages this model will be applied to a loop.

Logic: The Model receives an input raster image and calculates a centroid as a feature class. The centroid is then used in a spatial selection from a grid (PlatMap in this case). The selection of a single plat is then used by the clip tool.

Area of Trouble: the Python script that generates a point feature class that allows the selection of the specific plat I need.

Question: From the script presented below, what needs to change to get this to work? I have some experience coding in Python, but this is my first script using ArcPy. I could use more clarification on are how the input and output variables are called into and out of the script, as well as the write function for the created feature class.

First Attempt:

Currently my parameters (in properties for the script) are set to:

RasterInput  =  Raster Layer
RasterOutpur  = Feature Class

Codeblock is largely inspired by Blord-Castillo's response to a question here.

import arcpy, os

# Calculates the center of the raster image.

def RasterCenter(raster):

  #raster: string reference to raster

    raster = arcpy.Raster(raster)
    fcname = "{}_center".format(os.path.basename(str(raster)))
    x = raster.extent.XMin + (raster.extent.XMax - raster.extent.XMin)/2
    y = raster.extent.YMin + (raster.extent.YMax - raster.extent.YMin)/2

  #Combine the two to create a feature class containing the coordinates.

    featureclass = arcpy.CreateFeatureclass_management("in_memory", fcname, "POINT",spatial_reference = raster.spatialReference)
    with arcpy.da.InsertCursor(featureclass, ['SHAPE@XY']) as cursor:
        cursor.insertRow(((x, y),))
    RasterOutput = arcpy.setParameter(featureclass)

RasterInput = arcpy.GetParameterAsText(0)

RasterCenter(RasterInput)

Second Attempt:

In response to ahmadhanb's revisions.

I changed some other portions of the text to be clearer to me, but the biggest error I still receive when running the script is as follows. It is the same error I receive when running the revision from ahmadhanb. My suspicion is the the spatial_reference field, which I changed to the Authority code because it should be the same for every time I run this.

My parameters (in properties for the script) are set to:

rasterFile = Raster
Output_fc = Feature Class

import arcpy, os

# Calculates the center of the raster image.

def RasterCenter(rasterInput):

  #raster: string reference to raster
    raster = arcpy.Raster(rasterInput)
    item = os.path.split(rasterInput)
    name,ext = os.path.splitext(item[-1])
    fcname = "{0}_center".format(name)

  # Set the template file path.
    template = "C:/Users/preed/Desktop/STAGING/1799_Leak_Detection/LeakDetection.gdb/CentroidTemplate"

  # Setting Spatial reference to NAD 1983 StatePlane Washington North FIPS 4601 feet
  # Authority Code = 32048
    spatial_reference = arcpy.SpatialReference(32048)

  #Fetches the approximate X and Y coordinates we'll be using.
    x = raster.extent.XMin + (raster.extent.XMax - raster.extent.XMin)/2
    y = raster.extent.YMin + (raster.extent.YMax - raster.extent.YMin)/2

  #Combine the two to create a feature class containing the coordinates.
    featureclass = arcpy.CreateFeatureclass_management("in_memory", fcname, "POINT", template ,"DISABLED", "DISABLED", spatial_reference)
    with arcpy.da.InsertCursor(featureclass, ['SHAPE@XY']) as cursor:
        cursor.insertRow(((x, y),))

    return featureclass

rasterFile = arcpy.GetParameterAsText(0)

output_fc = RasterCenter(rasterFile)

save = arcpy.CopyFeatures_management(output_fc,arcpy.GetParameterAsText(1))


texts: Is there a way to create pixel centroids inside polygons with PostGIS ?

texts:
I am trying to find a solution to do the function of the QGIS-tool 'Generate points(pixel centroids) inside polygons' in PostGIS.

In the picture below you can see the output. Is there a way to do it in PostGIS?


texts: How to convert a CSV from DMS to WGS84 in QGIS?

texts:
I would like to automatically convert from a Degree - minute - second system to WGS84 EPSG:32632. I have a CSV with these kind of values:

LATITUDINE  LONGITUDINE
46°01’24,7” 11°13’47,5”
45°42’07,5” 10°55’11,3”
46°01’37,6” 11°06’41,7”
46°15’03,7” 11°11’00,1”

and I would like to automatically convert them.
There is a way? They are a lot, if not I would do with a normal converter...
I try to import them in QGIS (but only a few of them works) and export them in some way, but I sincerely do not understand how to. Any suggestion?
Thanks!


texts: Cluster Points on spatial network in R/QGIS

texts:
I am searching for a method for clustering spatial points along spatial lines.
The points are connected to the lines.

I tried nearest neighbour clustering and gProject/gInterpolate from the rgeos package. As you can see in the image, both results have the same problem: the clustering along lines doesn't make sense. In the image below, the point colours visualize the different clusters.

Is there a solution in R or in QGIS to cluster the points by their location on a line/network?


texts: Using Maximum Likelihood Classification and raster border in ArcGIS Desktop?

texts:
I have tried multiple times to see if this was an issue with my supervised classification, but I don't think it is. I created my training samples and created a signature file. Per my course, I need to run Maximum Likelihood. However, every time I do, my first class fills in extra border(s) around my results.

I am running ArcMap 10.3. I have attached an image of my Landsat5 imagery before I ran the Maximum Likelihood classification, as well as after (which is the issue I am having). If I had to guess, this is something very simple I missed.


texts:
I have a dataset with sensors measuring e.g. Nitrogen. Every sensor has an unique ID, but for every sensor there are multiple points (for each year one). However, when adding the new values for this year not all fields are filled resulting in the following schema:
  ID        Year       Location description
SensorA     2020
SensorA     2019       UndertheBridge

What I want to establish is that location description for 2020 is filled with the one of 2019, if 2019 is filled. How could I best establish that? My attempts so far in the field calculator have not the desired result.


texts: Deciding when to use projected coordinate system (PCS) or geographic coordinate system (GCS)?

texts:
While I understand the technical differences between them both, I'd like an expert's breakdown of criteria to look out for that would determine when I should use a PCS or GCS.

I understand that a PCS is an additional algorithm placed on a GCS to transform it from a 3D placement to a 2D placement on a map. For that reason, PCS units are metrics of length (e.g. meters or kilometers) to determine a file's location and shape. GCS on the other hand is its 3D placement on a spheroid, and its units are in degrees, minutes, and seconds.


texts: Geocortex Essentials 4.7 Issues with HTTP ArcGIS Services

texts:
Just upgraded to Geocortex Essentials 4.7 and am having some issues getting ArcGIS Services into it from an ArcGIS Server using http as oppose to https.

It allows me to connect and add the service, but in the view nothing will come up. I can't seem to get a good answer from Geocortex (I don't think they know) and I can't find the Admin or User Guide for 4.7.

I do know that http and https services are not allowed in the same Geocortex web app, but all I am adding are http services.
Https services work fine, no problems
Geocortex Essentials is https.


texts: Tobler's hiking function in PathDistance ArcGIS 10.3.1

texts:
I am working on least cost paths for an archaeological project and I am trying to implement Tobler’s hiking function on ArcGIS 10.3.1 for Desktop.

I am using SRTM 1 arc second  global downloaded on earth explorer.

My DEM characteristics are: GRID format, 16 bit signed integer, band 1, WGS 84 UTM zone 37N (Syria and Lebanon). I am working on a large zone, as I am trying to caclculate the least cost past between two cities around 200 km away from each other.

I am using Path Distance. My inputs are: feature source data: shp (point indicating the city); cost surface: slope in percent rise reclassified 8 classes; surface raster: my SRTM DEM. Vertical factors parameters: input vertical raster: my SRTM DEM; vertical factor: Tobler's table (this table http://kaitlinyanchar.com/wp-content/uploads/2015/04/ToblerAway.txt).

The least cost path I obtain using Tobler’s function coincides with the least cost path I obtain just using the slope as cost surface (i.e. slope in percent rise reclassified 8 classes and then I use cost distance and cost path to obtain a least cost path) without implementing Tobler’s function.

I have tried to produce a least cost path following Nicolas Tripcevich’s tutorial (i.e. without a cost raster) but I only get fairly evenly spaced lines that look like a spider web. I tried all the solutions I could find in different posts on gistackexchange or geonet (i.e. tried to use “same as display”, or change the DEM to 32 bit floating point, or convert my input point to a raster, or changed the maximum number of unique values to render in ArcMap options) but I invariably gets the spider web if I do not use a cost surface or if I use a cost surface with all values set to 1.

I tried to change the values in the table, just using random values in the negative slope and again I get the same least cost path that I get using only the slope or with the original Tobler's table.

There’s clearly a problem with the table because – as far as I can understand – the PathDistance does not take into account Tobler's table (the vertical factor) at all during the process.

Could it be a problem with the format of the table as described in Cost-Distance Tool not producing logical answer (see second answer)?

I do not have ArcGIS Pro so I cannot try to use it.

How can I try to change the table so PathDistance takes it into account?
Or how can I try to implement it using the command line in ArcGIS 10.3.1?

Please note that I am not very familiar with command lines and Python.


texts:
I'm using Mapbox for a web project and would like to have a drawing tool for closed polygons where the drawn line is always snapped to a grid.
Specifically, I'm using a 100m global grid (MGRS) and want for all the drawn shapes to stored by my server as an array of the MGRS grid squares.
Any suggestions either for a drawing tool that is only on the grid lines or for a different way to let users draw a polygon and to turn it into an array of grid squares?


texts: Unable to make "IDataStatistics:UniqueValues" work in .Net

texts:
I would like to use IDataStatistics:UniqueValues to iterate through a set of unique values.  This worked flawlessly in VB6 and VBA but I am having a difficult time making it work in VS 2013, ArcGIS 10.3, targeting .Net 3.5.

I am following the example shown here

Debugging InvalidCastException was unhandled by user code error from ArcObjects DataStat? is a similar post with no resolution. I've already tried the suggestions provided there to no avail.

The error I get always occurs on the  pEnumVar = CType(pData.UniqueValues, IEnumVariantSimple) line with Err.Description: "Unable to cast object of type 'System.Runtime.InteropServices.CustomMarshalers.EnumeratorViewOfEnumVariant' to type 'ESRI.ArcGIS.esriSystem.IEnumVariantSimple'."

I would really like to make this work rather than having to nest my loop to get the unique values the long way...

Here my code:

'pField As IField, pSelectionSet As ISelectionSet2 passed in

Dim pData As IDataStatistics = New DataStatistics
Dim pStatResults As IStatisticsResults
Dim pEnumVar As IEnumVariantSimple
Dim value As Object

... ...

pSelectionSet.Search(Nothing, False, pCursor)
pData = New DataStatistics
pData.Field = pField.Name
pData.Cursor = pCursor
EnumVar = CType(pData.UniqueValues, IEnumVariantSimple)
value = pEnumVar.Next
    Do Until IsDBNull(value)    ...     ...


texts: Spatial join between two FeatureCollections in Earth Engine

texts:
Given two FeatureCollections in Earth Engine, one being a regular grid, the other a number of country shapefiles, I would like to produce a new FeatureCollection with the geometry of the regular grid, but, in addition, the underlying country name, which is stored in a property of the country shapefile FeatureCollections.

In particular, as I am used in other GIS softwares, would it be possible to specify intersecting/overlapping as the condition for the spatial join?

Here is the EE demo link: https://code.earthengine.google.com/e50d27dd851e07f98cdabc69efc33da8


texts:
Where can I find or enable the processing toolbox in QGIS 2.18? Doesn't seem to be an option in the "View" drop menu, and the help file doesn't have anything about how to get it to display

I attached a screenshot of what QGIS looks like to me, maybe it's there somewhere and I'm missing it?


texts: ArcGIS Raster Calculator reference xy coordinates

texts:
I'm using ArcGIS Raster Calculator (Spatial Analyst Tools > Map Algebra > Raster Calculator). I've got the correct extension licences and have run basic Con(,,) successfully.

I would like to reference XY coordinates in my calculations but don't know how. Any suggestions?

My use case is to check 2 partially overlapping rasters and create a new one, with interpolation in the overlap region based on X coordinate.

Pseudo code below showing what I want to achieve.

[x] is where I want to reference the X coordinate of the cell being calculated
the numbers 7000 and 7500 are X coordinates relating to my raster overlap, leading to an interpolation function being defined

con ( IsNull("raster1") , "raster2" , con (IsNull("raster2"),"raster1", "raster1" + "raster2"*([x]-7000)/7500-7000 ))

I've got zero skill with Python, but am familiar with ModelBuilder and the old VB calculators etc. I'd prefer a solution that uses the GUI. If Python is the only way then I'll need a bit of hand holding to work out how to apply it and make it work.


texts: GeoServer WMS GetFeatureInfo only returning one Feature?

texts:
This is how I get feature info from a group layer when clicking on the map (using WMS GetFeatureInfo):

let view = this.map.getView()
let viewResolution = view.getResolution();
var source = this.groupLayer.getSource();

var url = source.getGetFeatureInfoUrl(
  this.map.getEventCoordinate(event),
  viewResolution,
  view.getProjection(),
  {'INFO_FORMAT': 'application/json'}
);

this.expedientesService.getExpsClick(url).subscribe(exps => {
  console.log(exps)
})

The problem: This request is returning only one feature even if there are more than one on that clicked point.

The URL generated by getGetFeatureInfoUrl has the following format:

http://foo?
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetFeatureInfo&
FORMAT=image%2Fpng&
TRANSPARENT=true&
QUERY_LAYERS=groupLayer&
LAYERS=groupLayer&
INFO_FORMAT=application%2Fjson&
I=50&
J=50&
CRS=EPSG%3Axyz&
STYLES=&
WIDTH=101&
HEIGHT=101&
BBOX=570384.8749920629%2C4722144.92875535%2C570386.7576818952%2C4722146.811445182

(I think I and J are always 50, don't know why)

  As seen here If there is more than one feature at the point then
  it may return information about all the features (up to
  featureCount) that are there.

Any idea of what is wrong here?

I have an auxiliary layer containing all the features, and when making a WFS GetFeature request, multiple features are returned.


texts: Dashed line with different colour fill and stroke in QGIS

texts:
How can I create a dashed line with a different fill and stroke colour?

I have a white dash line on my map "former range" which doesn't show up in the legend, which has a white background. I want to keep them both white. The only way to make the line show in the legend is to assign a darker stroke colour to border the dashes.

I've tried various formatting options, e.g. converting the lines to arrows and setting the heads at size 0 but I can only make the stroke dashed, not the fill. I've also tried with the layer as a line and a polygon.


texts: Write GeoPandas df to PostGIS - Linestring Z geometry

texts:
I am trying to write a geodataframe with Linestring Z geometry to a PostGIS table.

This rise the error:

ProgrammingError: (psycopg2.ProgrammingError) can't adapt type 'LineString'
[SQL: INSERT INTO test (geom) VALUES (ST_GeomFromEWKT(%(geom)s))]

Command

gdf.to_sql(table_name, engine, if_exists='append', index=False,
       dtype={'geom': Geometry(geometry_type='LINESTRINGZ', srid='32633',dimension='3')})

The geometry was created with

df['geom'] = df['geom'].apply(wkt.loads)
gdf = geopandas.GeoDataFrame(df, geometry='geom')

and this is the head

                                                geom
0  LINESTRING Z (272353.39 7040567.82 98.48, 2723...

I dont know what is failing. I though it was because my features has Z value. But from postgis doc about ST_GeomFromEWKT" This function supports 3d and will not drop the z-index".


texts:
How do I embed a live traffic map from google?

When I click "share or embed map", it says that the map is built for me, so other users will not view the same map as me. Is there a way around this?


texts: Are there any good libraries for visually creating simple Google Maps?

texts:
I'm working on a simple webapp that lets the journalists I work with easily construct styled Google maps.

In short, I need an interface that lets users place markers, set zoomlevel and choose a center lat/long — all without writing any JavaScript.

I know Google Maps has a "My Maps" functionality that provides this already, but I'm needing that kind of functionality as a way of adding data to my app.

Something like this...

Anyone know of such a library? I'm currently using a series of text fields using LocationPicker for the markers, but I still need to be able to set center and zoom — and I can't shake the feeling there's probably already something out there that does this...

Thanks!

Edit: This seems close, if a little clunky... http://www.birdtheme.org/useful/v3tool.html


texts: Transferring data from one raster to another in QGIS

texts:
I have 2 raster layers. The first is a mobile signal coverage raster with 0 and 1 values. The second is a population raster with value from 0.89 to 236. I want to transfer the population layer to the mobile signal coverage raster. I am using the raster calculation function in QGIS. Also, the extent that I use is the mobile signal extent. The calculation I am using is:
(“mobilesignal@1” >= 0) * “population@1"

But the result is not correct and the population numbers on the mobile signal raster are way too low.
Can anyone help me with the correct calculation?


texts: Transfer subtype Descriptions form long integer value into a string at new diffrent field

texts: I need to convert entire geodatebase into shp files.One of the layers have 30 subtypes. I want and maintain the subtype data. Is there a way to Automatically populate the name of the subtype (Transfer subtype Descriptions)in to a new string field in order to not do it manually and not lussing the subtype data ?


texts:
I am trying to load a base map (streets) and have a feature layer from the arcGIS online system.  I can view the base map, but I am receiving a few errors when i try to combine the base map and feature layer.

Here is my code...

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1, maximum-scale=1, user-scalable=no">
  <title>ArcGIS DevLabs: JavaScript Starter App</title>

<link rel="stylesheet" href="https://js.arcgis.com/4.9/esri/css/main.css">

What am i missing??


texts:
Am I missing something or is it not possible to read terrain data (such as forests, roads, rivers, lakes, etc) through google?

I'm trying to create an app where you can walk in a scene dynamically created from terrain data from google - live/semi live.

Is this possible?


texts:
Is it possible to view either live or very recent (past few days) aerial imagery for UK?

Ideally images would be satellite imagery, high resolution and free.


texts: Nearest Neighbour Classification in eCognition

texts:
I have a few papers and youtube videos explaining how to conduct the NN classification in eCognition but it isnt giving me any output.

I think the confusion is arising when trying to define the statistics in eCognition Developer.

Could someone like a step by step guide for dummies or explain it to me in very simple terms?


texts: Linking Google Sheet to Hosted Feature Layer in ArcGIS Developer Online (free version)

texts:
I am trying to create a ArcGIS Online Dashboard that is linked to Google Sheets. I have done this successfully in a regular ArcGIS Online map by publishing the Google Sheet as a CSV file and uploading accordingly; however, this type of layer cannot be read by the Dashboard since it is not a "Hosted Feature Layer".
Therefore, I tried to upload the Google Sheet using the following steps:

Navigated to "Add Item from Cloud Drive"
Selected the correct Google Sheet
Reviewed the data (coordinates, etc.)
Uploaded the Google Sheet

Despite following these steps, I constantly get a "An Error Occurred" message and the new hosted feature layer is blank. I have tried changing the "Share Settings" on the Google Drive, but to no avail.


texts: Binary cloud Mask from MODIS (MOD35_L2) using MRTSWATH

texts:
When the MRTSwath tool is fed with MOD35_L2 cloud mask HDF file, the tool provides six tif images (b0 till b5). Each tif file has 8 bit information.
Are these 6 images somehow indicating the 48 bits (https://modis-images.gsfc.nasa.gov/_docs/CMUSERSGUIDE.pdf) in MODIS cloud mask?

I am intersted only in the first 3 bits of the 48 bit mask. I need only the pixels which are 'cloudy' and 'uncertain clear'.

I have extracted the binary cloud mask using LDOPE 1.7 (by parsing bits 0-2==001 OR 0-2==101 ,Cloud_Mask SDS). The resulting HDF  file was fed to MRTSwath. But the mask derived like this looks too cloudy for my purpose.

syntax:

create_mask -of=output.hdf -on=255 -off=0
 -mask="input.hdf,Cloud_Mask,0-2==001,OR,input.hdf,Cloud_Mask,0-2==001"

Is there any alternate way to create a mask from the 6 tif images (MRTSwath)?


texts: Is there a way to use windowing functions (lag, lead) within the QGIS Field Calculator?

texts:
I am keen to compute a windowing function on values in QGIS. I was not able to find a lag/lead function. I think this could be emulated using @rownumber and $currentfeature, but I am not expert in using these.

I need to find a value for an attribute for a current feature. This seems to be possible using attribute($currentfeature,"fieldname"). I am unsure how to specify something like the previous feature, equivalent to treating the list of values for an attribute in a feature collection (pseudocode):

fc[@rownumber]["fieldname"]

I could then do fc[@rownumber-offset]["fieldname"], assuming it evaluates to NULL for the first entry. Is there any way to do this? I am sure people computing say acceleration on trajecotry data would deal with this.

kudos: if the answer hints at treatment of the sorting in formats that do not support this well (shapefile).


texts: Open source implementation of "feature density" algorithm

texts:
I'm working on a dataset of house searches - 1 polygon for each area in which people want to find a house from an online search tool.

This produces many overlapping polygons. It would be great to get a basic understanding of the places that are most frequently searched.

You can do this in QGIS and other programs using transparency:

But ideally, you would create a new layer with the search density for each cell in a grid. This is exactly what the "feature density" algorithm does in ArcView 3.2 and I've seen it work to great effect. It creates a raster map from a mass of overlaid polygons.

So the question is: is there an implementation of this algorithm for my polygon dataset for open source programs. I'm an R and QGIS user so one of those would be preferable. Many thanks!


texts:
What I am trying to do is, compare attributes in a shapefile to another updated version of that same shapefile just in a different folder. Right now I am just using the Compare Feature tool, which works okay, its just a rather long process when you have 3000 shapefiles to compare.  Instead of one file at a time, let it run through each of their respective folders and compare features that way.  I currently have code that will iterate through a folder, and I have the code that compares. Im just not sure how to put the two codes together. Is there a better way?
ArcGIS Desktop 10.2.1, Python 2.7.5

First code snippet

import arcpy
try:
 # Set local variables

    base_features = r"C:\Users\usr\Desktop\basepoints 2013\01001_points.shp"

    test_features = r"C:\Users\usr\Desktop\Testpoints apr2014\01001_points.shp"
    sort_field = "APN"
    compare_type = "GEOMETRY_ONLY"
    ignore_option = ""
    xy_tolerance = ".00000001 METERS"
    m_tolerance=".001"
    z_tolerance=".001"
    attribute_tolerance = "#"
    omit_field = "FID;Shape;PARCEL_ID;STATE_CODE;CNTY_CODE;APN2;OWNER;ADDR;CITY;STATE;ZIP;PLUS;STD_ADDR;STD_CITY;STD_STATE;STD_ZIP;STD_PLUS;TYPE_CODE;LONGITUDE;LATITUDE;CENSUS_TR;BLOCK_NBR;LOT_NBR;RANGE;TOWNSHIP;SECTION;QRT_SECT;FZ_COM_PNL;LAND_USE;CNTY_USE1;CNTY_USE2;M_HOME_IND;ZONING;PROP_IND;MUNI_NAME;MUNI_CODE;SUB_TR_NUM;SUB_PLT_BK;SUB_PLT_PG;SUB_NAME;FRONT_FT;DEPTH_FT;LAND_ACRES;LAND_SQ_FT;LOT_AREA;UBLD_SQ_FT;BLD_SQ_IND;FRPLCE_TYP;FRPLCE_NBR;FRPLCE_IND;PKG_SQ_FT;BSMT_SQ_FT;ADJ_SQ_FT;GR_SQ_FT;GF_SQ_FT;LIV_SQ_FT;BLD_SQ_FT;YR_BLT;EFF_YR_BLT;BEDROOMS;ROOMS;BATHS_CAL;BATHS;FULL_BATHS;HALF_BATHS;QTR1_BATHS;QTR3_BATHS;BATH_FIX;AIR_COND;BSMT_FNSH;BLD_CODE;BLD_IMP_CD;CONDITION;CONSTR_TYP;EXT_WALLS;FOUNDATION;FLOOR;FRAME;GARAGE;HEATING;PKG_SPACES;PKG_TYP;POOL;POOL_CD;QUALITY;ROOF_COVER;ROOF_TYP;STORY_CD;STORY_NBR;BLD_VIEW;LOC_INF;BLD_UNITS;STYLE;UNITS_NBR;ENERGY_USE;FUEL;SEWER;WATER"
    continue_compare = "CONTINUE_COMPARE"
    compare_file = r"C:\Users\usr\Desktop\compare_parcel_baja.txt"

    compare_result=arcpy.FeatureCompare_management(base_features, test_features, sort_field, compare_type, ignore_option, xy_tolerance, m_tolerance, z_tolerance, attribute_tolerance, omit_field, continue_compare, compare_file)
    print  compare_result.getOutput(1)
    print arcpy.GetMessages()
except:
    print arcpy.GetMessages()

Second code snippet

##======== iterating code ==============
import os.path

topdir= r"X:\CoreData_Jan15_2013\ParcelPoint\Data\points"

exten='.shp'

def step(ext,dirname,names):
    ext=ext.lower()

    for name in names:
        if name.lower().endswith(ext):
            print(os.path.join(dirname,name))

os.path.walk(topdir,step,exten)


texts: How to calculate the Zone of Impact using OpenWind

texts:
After downloading OpenWind, I am trying to calculate the ZVI for a single Turbine. There is no documention regarding this from the developers.

At this point I just want to compare the output of the visibilty analysis i.e the Viewshed.

I have loaded my DGM and defined it as my Elevation layer with Properties>Interpretation>Terrain Elevation.

I have a point shapefile which contains a single location for a turbine. And this is where I get lost. How do I define the turbine specifications and let the software know that the point is the turbine location?

Thanks for any pointers,

yours,

Robert


texts:
I have downloaded Scotland 50cm (phase 3) DTMs for a location I am studying using Digimap. I don't know what has happened to the tiles. When I added them from ArcCatalog I had to select create pyramids for the tiles.
If that is the problem, what can I do?
I am using ArcMap 10.8.1.


texts: Creating Evenly Spaced Nodes/Lines or Set Spacing Rule in ArcGIS Desktop?

texts:
I have written a Python script for ArcMap 10 that uses the x/y coords from a shapefile to create perpendicular transects to the existing line.  Since these shapefiles are actually rivers the lines are very curved.  In order to accurately represent the river's shape there are more points in the curves.  My problem is that the program I created takes all of these points to create multiple transects and in the curves it looks like a horrible mess because there are so many points.

Is there any way to create some kind of spacing rule where my program chooses evenly spaced nodes so I can avoid this messy overlap?


texts: Adding points on a map parcel from an attribute in the attributes table

texts:
I am using QGIS 3.4.2 and a have a map of a subdivision with roads, road names, and lot outlines. I also have added an attribute, by using a join to a CSV file, of whether there is a structure (house) on a parcel. What I would like to do next is add a point (not sure if this is the right word) or an icon, on a parcel if there is a house there.

Can anyone point me in the correct direction?


texts: How to do change detection between two Landsat scenes?

texts:
I want to do change detection and calculate the changes step by step using two images: the first image Landsat 1986, and the second Landsat 2000.

Actually I have used image differences in ArcMap and I get raster values for RGB ranging from 213 to -69.  What I should do to make this result more meaningful?


texts: Building GDAL With Read/Write Support of FileGDB

texts:
Has anyone been successful in building gdal with read/write support of file geodatabases?  I am able to build successfully, just never see FileGDB as in the supported vector formats.  Here is what I am doing:

Download latest gdal trunk (example: gdalsvn-trunk-20110822.zip)
Open the makegdal10.sln file in Visual Studio 2010
Edit the nmake.opt file, the ONLY thing I modify at this point is the GDAL_HOME variable, the FGDB section is still commented out.  Everything else left as is.
Build Solution
------>SUCCESS
Uncomment the FGDB Section and set path to FGDB SDK
Open Visual Studio command line tools and cd to /ogr/ogrsf_frmts/filegdb
Execute nmake /f makefile.vc plugin
------>SUCCESS - I now have ogr_FileGDB.dll in that folder

Now, I copy the ogr_FileGDB.dll and the FileGDBAPI.dll to the bin folder containing ogr2ogr.exe.  When I run ogr2ogr --formats, I never get FileGDB as a format.

Am I missing something here?


texts: GeoNetwork: Exception when deployed to Tomcat 7

texts:
I deployed GeoNetwork 2.10 to a Tomcat 7 instance, replacing the ROOT application. GeoNetwork is accessible from the command line, returning the usual HTML message requiring JavaScript. But when it is accessed with a graphical web browser it sends back the exception below. The various log files reporting exactly the same without adding any further information.

Any clues on what may be causing this exception are welcome.

  type Exception report

  message

  description The server encountered an internal error () that prevented
  it from fulfilling this request.

  exception

  javax.servlet.ServletException: Exception raised
    jeeves.server.JeevesEngine.init(JeevesEngine.java:236)
    jeeves.server.sources.http.JeevesServlet.init(JeevesServlet.java:91)
    javax.servlet.GenericServlet.init(GenericServlet.java:160)
    org.apache.catalina.authenticator.AuthenticatorBase.invoke(AuthenticatorBase.java:472)
    org.apache.catalina.valves.ErrorReportValve.invoke(ErrorReportValve.java:98)
    org.apache.catalina.valves.AccessLogValve.invoke(AccessLogValve.java:927)
    org.apache.catalina.connector.CoyoteAdapter.service(CoyoteAdapter.java:407)
    org.apache.coyote.http11.AbstractHttp11Processor.process(AbstractHttp11Processor.java:987)
    org.apache.coyote.AbstractProtocol$AbstractConnectionHandler.process(AbstractProtocol.java:579)
    org.apache.tomcat.util.net.JIoEndpoint$SocketProcessor.run(JIoEndpoint.java:309)
    java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1146)
    java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:615)
    java.lang.Thread.run(Thread.java:701)

  root cause

  java.lang.StringIndexOutOfBoundsException: String index out of range:
  -1    java.lang.String.substring(String.java:1949)    java.lang.String.substring(String.java:1916)
    jeeves.monitor.MonitorManager.(MonitorManager.java:76)
    jeeves.server.JeevesEngine.init(JeevesEngine.java:141)
    jeeves.server.sources.http.JeevesServlet.init(JeevesServlet.java:91)
    javax.servlet.GenericServlet.init(GenericServlet.java:160)
    org.apache.catalina.authenticator.AuthenticatorBase.invoke(AuthenticatorBase.java:472)
    org.apache.catalina.valves.ErrorReportValve.invoke(ErrorReportValve.java:98)
    org.apache.catalina.valves.AccessLogValve.invoke(AccessLogValve.java:927)
    org.apache.catalina.connector.CoyoteAdapter.service(CoyoteAdapter.java:407)
    org.apache.coyote.http11.AbstractHttp11Processor.process(AbstractHttp11Processor.java:987)
    org.apache.coyote.AbstractProtocol$AbstractConnectionHandler.process(AbstractProtocol.java:579)
    org.apache.tomcat.util.net.JIoEndpoint$SocketProcessor.run(JIoEndpoint.java:309)
    java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1146)
    java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:615)
    java.lang.Thread.run(Thread.java:701)

  note The full stack trace of
  the root cause is available in the Apache Tomcat/7.0.26 logs.

Update: Following Délawen's answer, this is now issue 443 at GitHub. For now I will leave this question marked as unanswered.


texts: Removing duplicate polygons having different attribute values using ArcPy?

texts:
I'm trying to create a national.shp to cover the entire countries land. I made index called "province" to help me update the information in my attribute table from the national shapefile.

The idea is to create each provincial polygon files first and then merge into the national file.

The problem is that after the data merging, I found that an overlap exists for my file, (same polygon, two values). Is there a way to get rid of one of the values?

I tried to use dissolve tool, but that way I will have to lose my "province" field, which will be needed if I want to update my national.shp... I tried to trace the province value back, just can't seem to find a method.


texts:
I can not get it to work. I started with complex sql and ended with 'hello word' example, but I get errors:
My sqls are:

SELECT t.id, t.user_id, t.coordinates, t.size, t.price FROM my_table t WHERE t.id=%id_param%
Error: Failed to create SQL view: ERROR: syntax error at or near ")" Position: 100
SELECT t.id, t.user_id, t.coordinates, t.size, t.price FROM my_table t WHERE t.id='%id_param%'
Error: Failed to create SQL view: ERROR: invalid input syntax for integer: "" Position: 100

I am not sure, what I am missing.
I am using approx 2 months old OpenGeo suite (4.0.2).

  About GeoServer General information about GeoServer Build Information Version
  2.4-SNAPSHOT Git Revision f295d37d70c3b1f0e51942e8b031253e4d75459 Build Date 20-Jan-2014 22:29 GeoTools Version 10-SNAPSHOT (rev
  3bb3df4f410750baf9ebc48c9c475a70642aadb6) GeoWebCache Version
  1.5.0 (rev 1.5.x/f981cafd349886bb2a1ef45067f26b1aa8e4aff4)


texts: Forcing color ramps to center around zero in QGIS

texts:
I have several columns of data that have percentage values. Some are positive and some are negative. In one particular column, the values range from -0.70 to +0.20. The mean value of the column is around -0.08.

I want to color these features using a divergent color ramp (for example, the Spectral ramp). But the symbology options in QGIS seem to not give me an option where the actual color divergence happens, so I get classifications like these:

Note how the yellow color is arbitrarily centered on -0.1, not 0.

How can I force the central color of the color ramp to be centered on 0?

For example, how can I tell QGIS to do "Pretty breaks", but on a given range instead of the data column itself? Do I have to do this manually?


texts: GeoServer + Leaflet display map wrongly in square zones of contrast

texts:
I have a GeoServer(2.16.2) publishing a raster layer through WMS(1.1.1). On the client side, I have a Leaflet(1.7.1) canvas display the raster layer.
As can be seen from the figures below, map is being displayed unbalanced in intensity/bright within the tiles; At each level the tiles are wrongly balanced.
This is the code I'm running on the client side:

Any ideas on how to fix this?


texts:
I'm using Arc GIS Toolbox Hydrology...everything works fine. Moreover, it is possible to derive the potential stream velocity just out of the DEM? Informations about rainfall, infiltration, hydraulic radius or surface roughness are not available and the project should be kept easy!

I need streamclasses like: fast velocity, slow velocity and standing water.

Does anybody have an idea how to derive such a classification just out of a DEM?


texts: Finding which point is furthest from something using ArcGIS for Desktop?

texts:
Here is an example of what I am looking to do:

Which town in South Dakota is farthest from a Hospital?

Which tool would I use to be able to do this?


texts: Summing number of intersections with separate buffer layer using QGIS?

texts:
I have two layers: 1) census blocks, and 2) point data.

My goal is to calculate the sum of the point data within the census blocks to eventually create a choropleth map and to run regression models with the sums of the points as a dependent variable; however, some of the points are directly between blocks and are not truly within them, so the points in polygons tool seems to miss a lot of points.

Because I am interested in the effect of these points on the proximal blocks, I had the idea of creating a 15m buffer around the points, and summing up the number of buffers that each census block intersected, but I am at a loss for how to do this in QGIS. So far, I have tried to join attributes by location, intersections, unions, and they all don't seem to get me where I need to be.

I attached a sample picture of my layout with the intersecting census blocks selected.