texts: Arcmap layout mode increasing label size?

texts:
I have a layer of rivers, but I only want to include some of them as including all would clutter the map. So i create a new layer by selection, and add that to the layout, add the labels, but they are too small. So i go into properties, symbol, and change the size to larger.. but when i change it to larger all the rivers appear instead of just the selection.

How can I get only the selection to have large labels?


texts: Field Calculator - Calculating Sum based on Different Fields

texts:
I'll explain this the best I can. I have 3 fields that I need to work with in Field Calculator:

Link ID
DesignLength
AsBuiltLength

DesignLength is length of each line feature calculated by ArcMap. Each LinkID instance has its own specific DesignLength. However, there are ID numbers that are used multiple times. What I need to do is sum the DesignLength for one number (that could have multiple rows) into the AsBuiltLength field using Field Calculator.

Here is a screenshot of my table.

So for example, 0102 would have the same number in both DesignLength and AsBuiltLength because there is only one instance of 0102. But for 0104, I need all 13 instances of 0104 to be added together, and then that same number to be put in AsBuiltLength in all 13 rows for 0104.

Any ideas? I'm struggling to get anywhere successful in Field Calculator.


texts: Exploring ArcObjects through C# or VB from Python background

texts:
While I understand that this may be off topic to an extent, I am in a bit of a bind and could use some help directing my thought process and research.

I have little programming background aside from:
1) Scripting in a language similar to C
2) Practical use of Python for use in geoprocessing with ArcGIS

This has worked just fine for me so far, but there appears to be a need for me to move into a .NET-based language so that I can start making use of the tools exposed in ArcObjects that ESRI has not yet made available through Python. A decision is being made by my management (with some direction from me) on whether the company should move forward and create a standard for coding in ArcObjects that is aimed at using VB or C#. The main considerations (since I am the only person who will be taking part in this) are future extendability, "safety" of code, and readability. Most of these concerns appear to be things that can be dealt with in the writing of the code as opposed to the language itself, but I don't know enough to be certain.

I've read an article or two comparing C# to VB in terms of .NET, but I would like to ask for further reasoning to consider when deciding on one language over another for a person of my background. What are the most important things to keep in mind when choosing which of the two languages to begin study in for use with ArcObjects?

Sorry if this is too verbose! I simply want my need to be understood.

Thanks.

EDIT: Responses have been great. I went with Dan's answer since it covered all the ground I asked about and went the extra step to point me at IPython.


texts: Georeferencing CAD drawing and setting ground control points in GRASS?

texts:
I'm fairly new to QGIS and Grass GIS, and learning the interface.

I have a CAD drawing that has no coordinate system assigned and needs to be georeferenced to NAD83 HARN Northern Colorado (feet).

I'm trying to use v.rectify to match known control points on the CAD drawing to their counterparts on a properly georeferenced vector map downloaded from the county website.

I keep getting the error:

  Poorly placed control points. Can not generate the transformation
  equation.

So I assume I'm not going about setting these points properly.

I'm using the Digitize New Points function in GRASS to get the source coordinates at road intersections on the CAD drawing and in the County data, then putting them in notepad and using that for my input file.


texts: Is Landsat 4-5-7 L1G geometric corrected?

texts:
I'm using Landsat L1G for our project and need L1G to be geometric corrected for later use. One thing I know for sure that L1G need to be employed a Digital Elevation Model (DEM) for topographic accuracy. I follow a guide to use ArcMap (from ArcGIS) to correct the uncorrected Landsat image. But the process required DEM and DOQ (Digital orthophoto quadrangle) image to run. DEM is available but DOQ is not available outside US.

As I read, ArcMap using DOQ to establish the GCPs for uncorrected image. Since the DOQ is not available, I wondering if only using DEM to correct the L1G image is good enough? In other hands, can anyone know that the Landsat L1G is corrected with GCPs?

Updated:

We decided to use data from L1T source instead. But I face another problem can't find any information about satellite observation zenith angle/azimuth of Landsat image. We need it to make sure two different image has the same condition for further calculation. Any one have idea where I can find that information?


texts:
I am interested in estimating road width (distance from road center line to the nearest building on both sides) in PostGIS. I am aware of the cases/exceptions where there can be buildings on only one side of the road (or no buildings at all). Thus in the following query, it is assumed that buildings exist on both sides of the road center line within 100 meters search radius. I have two tables i.e., road st_wi and building bu_wi in PostgreSQL 9.5 database. The building table contains gid, geom and height columns for each building. The sample scenario is given as follows:

The Approach:

The approach I am adopting to estimate street width is:

Buffer road to 1m to start
Intersect that buffer with building polygon
Count the number of rows
if count = 0,increase buffer by 1m
Stop if two buildings (on opposite sides) are selected and return that buffer distance

The Query:

So far, the following query works but have many limitations. For example, the buffer around road center line increments up to 15 meters only.

Do  Language plpgsql
$$
Declare
r integer;
Begin
    EXECUTE 'Drop table if exists t';
    EXECUTE 'Create table t (gid int, geom geometry)';

   /* buffer upto 15 meters */
    For r in 1..15
    LOOP
        EXECUTE '
        Insert into t (gid, geom)
        Select
            b.gid,
            b.geom
        From st_wi a
        -- joining buildings to road within 100m search radius
        Left join bu_wi b on st_dwithin(b.geom, a.geom, 100)
        Where st_intersects(st_buffer(b.geom,$1), a.geom)
        and a.gid not in (select gid from t)
        and b.geom && a.geom
        group by b.gid, b.geom
        having count(a.geom) = 2'
        using r;
    END LOOP;
END;
$$;

Then:

Select * from t;

The above approach might not be the most efficient one. I am open  for suggestions.


texts: How should students cite maps I create for them?

texts:
I teach QGIS within my high school's IB Geography course at our school here in China and I'm struggling to get accurate base maps.

I'm considering using a drone to take aerial photos of the village we study and tracing it in QGIS to create a shapefile for the students to use as a base map.

How would students cite that map/shapefile?

I've checked for sample citations and directions for maps, and found instructions here, but I'm struggling to build an appropriate MLA citation for the students to use. I'm guessing it should include author name, map title, date of creation, and location accessed (school's online learning platform). We're also supposed to anonymise the work, so we can't put our school's name in the citation.


texts:
I have a shapefile with polygons and another one with points. I would like to create a new shapefile with those polygons that don't contain any point inside. I tried the function sp::over in R and it returns a vector indicating which polygons have points, but I have been unable to create the new shapefile from there.

 library(sp)
 polygons <-over(basins,points_proj)


texts: ExtractSubDataset in ArcPy loop giving ERROR 999999?

texts:
In ArcMap, I've created a layer group that I would like to add extracted subdatasets from a file to.  It is a .grb file and I've successfully done this with the provided "Extract Subdataset" tool.  Is it possible to do this in python with a loop to name each out_raster in the group according to index?  For example,
for subdataset_index = "00", then
out_raster = "filename+"010"+" and so on for all indices in the original file?  I've seen this done manually with the batch method, but I have 44 subdatasets to extract in each file.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
idx = 0
while idx < 10:
  in_ras = "C:/Data/gribFile.grb"
  out_ras = "C:/ArcGIS/geoDB.gdb/Grb_+'FL'+idx+'10'"
  arcpy.ExtractSubDataset_management(in_ras, out_ras, idx)
  idx = idx+1

I guess I'm getting stuck on how to grab each index from each file as object.  Below is error I'm getting.

  Error message: Runtime error  Traceback (most recent call last):
  File "", line 8, in    File "c:\program files
  (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 14129, in
  ExtractSubDataset     raise e ExecuteError: ERROR 999999: Error
  executing function. The table name is invalid. No spatial reference
  exists. Invalid function arguments Invalid function arguments An
  invalid SQL statement was used. [SELECT ObjectID, UUID, Type FROM
  GDB_Items WHERE PhysicalName = 'Grb_+'FL'+IDX+'10'' AND Type IN
  ('{70737809-852C-4A03-9E22-2CECEA5B9BFA}')] Failed to execute
  (ExtractSubDataset).

Was able to get it working with this:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
#Levels array
lvs = {myArray}

for ftLevel in lvs:
level = ftLevel
FL = lvs[ftLevel]
gFile = "C:/filepath"
gFileDate = "G_2016"
gFileHour = "000HR"
gRasterFileloc = "C:/filepath/raster"
in_ras = gFile+gFileDate+gFileHour+".grb"
out_ras = gRasterFileloc+gFileDate+gFileHour+FL
arcpy.ExtractSubDataset_management(in_ras, out_ras, level)


texts: Unable to get the feature collection data while exporting as kml file

texts:
I am trying to get the kml file by giving simple polygon data but not able get exact kml Please suggest me if iam wrong.

code snippet:

GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory(null);

WKTReader reader = new WKTReader( geometryFactory );

Encoder lEncoder = new Encoder(new KMLConfiguration());

OutputStream lFileOutputStream = null;

lEncoder.setIndenting(true);

try
{
    lFileOutputStream = new FileOutputStream("my.kml");

    Polygon polygon = (Polygon) reader.read("POLYGON((20 10, 30 0, 40 10, 30 20, 20 10))");
    lEncoder.encode(polygon, KML.kml, lFileOutputStream);
    lFileOutputStream.close();
}
catch (FileNotFoundException e1)
{
    e1.printStackTrace();
}
catch (ParseException e)
{
    e.printStackTrace();
}
catch (IOException e)
{
    e.printStackTrace();
}

output getting in KML File

< ?xml version="1.0" encoding="UTF-8"?>
< kml:kml xmlns:atom="http://www.w3.org/2005/Atom"
    xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:xal="urn:oasis:names:tc:ciq:xsdschema:xAL:2.0"/>


texts: GRASS Error: WARNING: segment zero_fill(): Unable to write (No space left on device)

texts:
I was trying to define the drainage patterns and to delineate some catchments based on a 30m DEM using

r.watershed --overwrite elev=DEM accumulation=Accum drainage=Drainage -m

but grass would constantly crash and i would obtain the following error

WARNING: segment zero_fill(): Unable to write (No space left on device)
WARNING: seg_open(): could not write segment file

I did some searching about this error here and here, and this help me understanding that the program was running out of memory when processing the order, however it did not provide much more insight into how to address this issue.

Does anyone know what i can do to solve this and run r.watershed?


texts: Ungrouping French Guiana from France in Natural Earth data using QGIS?

texts:
I have a layer of countries, and I am trying to select all of Central and South America and export it to GeoJSON.

I am using this Admin 0 - Countries without boundary lakes from this link:
https://www.naturalearthdata.com/downloads/110m-cultural-vectors/

For some reason, I cannot select French Guiana alone. When I select it (with either an expression or the point-and-click), France gets selected as well. I searched for French Guiana in the attribute table and I couldn't find it, I'm afraid it is grouped in with France somehow.

How do I separate French Guiana from France so I can select all of Central and South America and not select European France?


texts: Portions of edited shapefile will not draw

texts:
I edited a shapefile of approximately 1500+ pipe segments yesterday in QGIS.  When opened today (first in ArcMAP by another user) only about 10% of the pipes were drawn.  The same happened when the shapefile was subsequently opened again in QGIS.  The attribute table still shows the 1500+ pipes.

Any help as to why they are not drawing?

Attached image shows original file (rose coloured lines) versus edited (blue lines).


texts: I want to delete a polygon and redraw but I want to save the attributes?

texts: I am trying to redo a polygon it is so messed up that it can't just be fixed it needs to be delete and redrawn but I don't know how to do this and save the attributes.


texts: How to run OSGeo4W shell for molodensky transform?

texts:
I'm a beginner with Proj, so I just downloaded proj, and I'm using the OSGeo4W Shell and when I try to run the example

proj=molodensky a=6378160 rf=298.25 da=-23 df=-8.120449e-8  dx=-134 dy=-48 dz=149 abridged that is listed on the Proj website I get the error:

projection initialization failure
cause: invalid projection system error (-9999)
program abnormally terminated

I've tried looking up this error message but I can't find anything about it. Do I need to download something else to get the command to work? I'm on a Windows 10 computer.


texts: How to split multiband image into image tiles using Rasterio?

texts:
I have been using this answer on GIS SE to split 3-band (RGB) imagery into 256x256 3-band image tiles:

import os, gdal

in_path = '/path/to/indata/'
input_filename = 'dtm_5.tif'

out_path = '/path/to/output_folder/'
output_filename = 'tile_'

tile_size_x = 256
tile_size_y = 256

ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
xsize = band.XSize
ysize = band.YSize

for i in range(0, xsize, tile_size_x):
    for j in range(0, ysize, tile_size_y):
        com_string = "gdal_translate -of GTIFF -srcwin " + str(i)+ ", " + str(j) + ", " + str(tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) + str(input_filename) + " " + str(out_path) + str(output_filename) + str(i) + "_" + str(j) + ".tif"
        os.system(com_string)

Is there a comparable Rasterio or numpy approach?


texts: Create one id variable for each dataframe in a list of dataframes with name of each dataframe

texts:
I have a list of dataframes. For each dataframe I want to create one new variable or column called uniqueid and I want this unique id to be the name of the dataframe from the list. I have tried the following but it does not work.
# List all shapefiles
dfs <- Filter(function(x) is(x, "data.frame"), mget(ls()))

for( i in seq_along(dfs)){
  dfs[[i]]$uniqueid<-print(paste(dfs[[i]]))
}

When I do the following each dataframe get a uniqueid equals to ID. SO I guess the part that is not working correctly on my code is print(paste(dfs[[i]])) however I can't figure out what the right fix is.
# List all shapefiles
dfs <- Filter(function(x) is(x, "data.frame"), mget(ls()))

for( i in seq_along(dfs)){
  dfs[[i]]$uniqueid<-'ID'
}


texts: ESRI ASCII. Latitude and longitude in a in raster

texts:
I've started a new project that involves a ESRI ASCII grid. As far as i know it is a map of Germany or the UK (I'm not sure which one).
This is how the header looks like:

NCOLS 8083
NROWS 12952
XLLCORNER -59900
YLLCORNER -58100
CELLSIZE 100
NODATA_value 0

I decided to use POSTGIS to import it into a raster table (I named the table testraster). The import worked smoothly, and I managed to play around with it, using some postgis functions.

For example:

SELECT x, y, val, ST_AsText(geom) FROM (SELECT (ST_PixelAsCentroids(rast, 1, true)).*
FROM testraster WHERE rid = 1) foo;

This select returns all the center points of my cells(cell = a value in the raster). Only the cells that have values are taken in consideration.
The very small part of the result looks like this:

x       y       val     st_astext
4955    353 -109    POINT(435550 -51050)
4956    353 -112    POINT(435650 -51050)
4957    353 -112    POINT(435750 -51050)
4958    353 -114    POINT(435850 -51050)
4959    353 -115    POINT(435950 -51050)

This takes me to my question.

How do i know the latitude and longitude of any of these points? What is the equivalent of 435940 and -51050?
Or, a more realistic scenario: I'm in Berlin at (52.52001, 13.40495), what is the cell that i should be looking for in my raster?


texts: Forcing to include Fields in WMS/WFS in pyQGIS (processing)?

texts:
I've got a *.qgs project I have to publish in QGIS server with 60 layers. I don't want to publish every fields so I made a processing script to exclude them using QgsVectorLayer.setExcludeAttributesWMS() and QgsVectorLayer.setExcludeAttributesWFS().

I would like to know how to do the opposite, a kind of setIncludeAttributesWFS.

##WMS=group

from qgis.core import *
from PyQt4.QtCore import *

myRegistry = QgsMapLayerRegistry.instance().mapLayers()

myList=[]

for layer in myRegistry.values():
    if (layer.type() == QgsMapLayer.VectorLayer):
        #default cleanning
        myList[:] = []

        fields = layer.pendingFields()

        for field in attributs:
                if field .name() != 'type':
                myList.append(field.name())

        #exclude WFS/WMS
        layer.setExcludeAttributesWMS(myList)
        layer.setExcludeAttributesWFS(myList)


texts: PyQGIS importing CSV changes projection of basemap

texts:
I am using PyQGIS to import a csv file with a lat and long, when doing this I am using the appropriate crs of EPSG:4326. I'm plotting this onto Google Maps.

I load my basemap, then import my CSV. The issue is that my basemap projection then changes to 4326 and I need it to remain on 3857.

I've tried importing the basemap after the CSV and moving it down in the layers, however this still changes the projections.

import requests
from PyQt5.QtGui import *
from PyQt5.QtCore import *
from qgis.core import *
from qgis.utils import iface
from qgis import core

#Use Google Street Map as QGIS basemap.

service_url = "mt1.google.com/vt/lyrs=m&x={x}&y={y}&z={z}"
service_uri = "type=xyz&zmin=0&zmax=21&url=https://"+requests.utils.quote(service_url)
tms_layer = iface.addRasterLayer(service_uri, "GoogleSat", "wms")

#Import CSV and plot.
uri = 'file:///home/user/fred.csv?type=csv&xField=%s&yField=%s&crs=%s' % ("Site Longitude", "Site Latitude", "EPSG:4326")
layer_csv = QgsVectorLayer(uri, 'fred', 'delimitedtext')
layer_csv.isValid()
QgsProject.instance().addMapLayer(layer_csv

I'll be the first to admit I'm a novice with QGIS.


texts: HTML Custom Attribute does not register as HTML code

texts:
I am working on creating a hyperlink by configuring the popup through Attribute Expressions. The purpose is to create a telephone call link for all of the Municipalites. I have all of the numbers for all of the police departments  and want to insert them into some HTML text which will be populated by creating a virtual field in Attribute Expressions.

When I have the full HTML in the field, the hyperlink is created, but when I write out the code by using Attribute Expressions below the HTML text just shows up with no link.

I would prefer to keep using a list of attributes instead of a custom display as I think this looks more professional, and would like to populate the HTML code as apposed to changing the HTML every time so that people more involved with the call center activities can make changes to the phone numbers.

var phone = $feature.Police_Num1

'' + phone + ''

Has anyone else had this issue or know how to solve for it?

UPDATED


texts: QGIS Problem loading or adding WFS Layer

texts:
I have a problem loading and adding data from the following a WFS to my QGIS project:
https://data.wien.gv.at/daten/geo?version=1.1.0
I can connect to the service and see the data available. (See picture below)

But after adding a layer with the data to my project I get the following error and no data is displayed:

Layer ogdwien:FMZKGEBOGD: Download of features for layer
ogdwien:FMZKGEBOGD failed or partially failed: Server generated an
exception in GetFeature response:
org.opengis.referencing.NoSuchAuthorityCodeException: No authority was
defined for code "". Did you forget "AUTHORITY:NUMBER"? No authority
was defined for code "". Did you forget "AUTHORITY:NUMBER"?. You may
attempt reloading the layer with F5

As far as I know you don’t need an authority number to use the service since it is an open data service from the city.
I can export the data to a shapefile. But then some of the data is missing. In my case not all the buildings are exported, just some of them.
I am using Windows 10 and the QGIS Version before was 3.4.11. I now upgraded to the new version 3.14.1 but I still have the same problem.


texts: Calculating feature count of attribute type using arcade in ArcGIS Online

texts:
I have a simple table saved as a feature in ArcGIS Online. The table structure is very simple:

ID,Date,Type,Location,Total,Sum_typ,cumul

In the last three fields I would like to calculate statistics using arcade:

FID,Date,Type,Location,Total,Sum_typ,sum_loc
0,19.03.2020,confirmed,Hamburg,8,5,3
2,19.03.2020,confirmed,Hannover,8,5,3
3,21.03.2020,not confirmed,Hamburg,8,3,3
4,21.03.2020,not confirmed,Hannover,8,3,3
5,21.03.2020,confirmed,Bremen,8,5,2
6,24.03.2020,not confirmed,Hannover,8,3,3
7,25.03.2020,confirmed,Bremen,8,5,2
8,25.03.2020,confirmed,Hamburg,8,5,3

Total= would be the total number of records in the table which the following arcade calculation would work.

 var tot = sum($layer)

Here I need help

Sum_typ = should contain the total number or records grouped by type

-any record which is confirmed should have a 5 as this is the total amount of "confirmed" cases in the data.

-any record which is not confirmed should have a 3 as this is the total amount of "not confirmed" cases in the data.

sum_loc = total number of records grouped by location
The same as above bur groupd according to location:

therefore any data from Hamburg = 3
therefore any data from Hannover = 3
therefore any data from Bremen =2

I can't do this in any other way as far as I can see apart from with arcade. I need to keep the data dynamic so using a summarize function does not work as this creates a static view which is not updated when new records are added.


texts: Is this correct behavior for ST_ConcaveHull (PostGIS)?

texts:
I'm getting this error on the ST_ConcaveHull function call and after reading the docs, I'm unclear as to whether this Error is correct behavior or not.

Here is the business end of the query (vx are points):
select id , ST_ConcaveHull(ST_Collect(vx), 0.80)

And here's the error:

********** Error **********
ERROR: ExteriorRing: geom is not a polygon
SQL state: XX000
Context: PL/pgSQL function st_concavehull(geometry,double precision,boolean) line 67 at assignment

Since I'm not passing polygons of any sort, I assume the ST_ConcaveHull function is creating them internally.  I assume that there must be 1 or more grouping of points that don't produce an expected polygon.  And this wouldn't surprise me, but I'd rather just get back a null geometry rather than an error that borks the query.

Is the behaviour as-designed, and if so, what is the best way to isolate the data that is causing the error (I assume I would need to pre-filter the data somehow)?

I can also tell you that
select id , ST_ConcaveHull(ST_Collect(vx), 0.99)
 completes successfully.

PostGIS Version:

  POSTGIS="2.0.4" GEOS="3.4.2-CAPI-1.8.2 r0" PROJ="Rel. 4.8.0, 6 March
  2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.7.8"
  LIBJSON="UNKNOWN" (core procs from "2.0.4" need upgrade) RASTER
  (raster procs from "2.0.4" need upgrade)

Postgresql Version:

  PostgreSQL 9.2.2, compiled by Visual C++ build 1600, 64-bit


texts: Creating Land Use Change Map in ArcGIS for Desktop?

texts:
I need to create a map with multiple dataframes showing the percent change in land use - one each for agriculture, developed, and forest. Percent change has to be labeled on each county and color coded for increase, decrease, and no change.  I have the data figured out and ready to map.

I can't get started with creating the map and was wondering if you could give me a boost.

How can I separate the fields and just map the data picked from the table shown?

I included a screen capture of the table and basic county map.

I'm using ArcGIS 10.2 for Desktop.


texts: How can I export elevation layer data and OpenStreetMap layer data as shapefile in QGIS?

texts:
I am a newbie in QGIS looking to obtain elevation data with respect to a map. As I could not find any free elevation data, therefore I adopted the following approach.

Elevation data: I got the elevation data by selecting an area of interest on the ESRI satellite map and used the SRTM Downloader to extract the elevation data. The elevation data file was downloaded as .hgt format. I opened this file using qgis to see the elevation of the area of interest.
OpenStreetMap: I mounted the OpenStreetMap on top of the elevation data.

Now I can see the elevation on the map along with street names.

My aim is to be able to export this data to another software which takes shapefiles. How can I package the elevation data and OpenStreetMap as shape file and export it? One thing I did not understand is that I was not able to see the attributes table for elevation data file. Is my method correct? Is there a better way or source to achieve the above? I am not from the GIS field at all, but I need the data.

I am running QGIS 3.10 on an Ubunut18 machine.

Edit version:
After the direction of GBG, I achieved the following:(How to save this as a shapefile?)


texts: Different models, different contours: are they comparable?

texts:
I will get contour lines of electromagnetic field strength around power lines. The contours will be calculated by different agencies using different models. We need to check whether the resulting zones (x meter from power line) around the power lines are similar between models/agencies. There is no 'gold standard' because it is not possible to actually measure the field strength on the ground. Furthermore, the zone is not simply a straight line because in some situations two power lines cross or approach each other and we are most interested in these deviating situations.
Is there a way in GIS to compare the modelling results (polylines/contours) among each other?


texts: What are the characteristics of a raster to Netcdf layer in ArcMap?

texts: Specifically, when trying to edit the attribute table? This is would be for work purely in ArcMap, and converting raster to Netcdf. I have a land cover layer (that I've converted to a raster to convert to netcdf) that has an column from its attribute table that I need in the Netcdf attribute table. I have tried joining the land cover layer with the raster and then converting it to netcdf, but the attribute table for the netcdf only has OID, VALUE, and COUNT.


texts:
Does anyone know of any tools that exist for creating a polygons.

I need to provide some employees with the ability to create shapes that detail the boundary of certain locations. So I just want a simple tool that they will simply be able to use to create polygons. (I'll then let them import the points to my main application).

(It wouldn't take a massive amount of work to create this myself in Google Maps, but would love if that could be avoided and something already exists).

Update with more information: (Sorry for missing out obvious key points).

The main application is very basic. It's really just a frontend to a vacation rental portal that allows travellers to search for locations and then return properties within or around that location. I use a mix of tools to do the searching PostGIS and Solr with spatial search depending on the actual module as the backend changes...

What happens is, I take a location search string, and then attempt to find the best matching WOEID (using Yahoo GeoPlanet) for that location string. At the moment, we've got about 200k different polygons coords which have been assigned to their respective WOEID. The shapes are then drawn on as a polygon to a single Google Map and it's used to help travellers visualize there search boundaries.

Because we have now exhausted  of the available shape files that I can find to easily match up against WOEIDs (such as flickr dataset), we are going to start manually creating the polygons and assign them to each WOEID.

So I'm just looking for a tool that I can set up our employees with to create the polygons. I realise that they won't be perfect, but this is really only for visualization. So it's not important that they are perfect..


texts:
I am trying to add a custom crs in QGIS for GDA2020 using grid file GDA94_GDA2020_conformal.gsb. My parameter is
proj=utm +zone=51 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +nadgrids=GDA94_GDA2020_conformal.gsb

My test point is
North: -21.7403935833333
East: 121.197650527778
When I press "Calculate", it just says "Error" in the result box. I checked and the file is there in C:\OSGeo4W64\share\proj
Is there anyway to know what the error is?


texts: How to count points within adjoining hex grid cells?

texts:
I am currently using QGIS, and have developed a vector grid (hexagonal). I have a separate layer with point data. For each grid cell I need to count the number of points in that cell (but also all adjoining cells).

I searched for this problem on Google and Stack Exchange, but no solutions I've been able to find nothing.

Some diagrams to explain, below.

Image 1: I want to calculate a value for this grid cell

Image 2: Value for selected cell  = the number of points within this cell + all adjoining cells


texts: Handling MOD13Q1 NDVI Product Quality Assessment (QA) flags?

texts:
I am working on MOD13Q1 NDVI Product and i want to apply Savitzky Golay filter for smoothing  Ndvi time series  using Quality Assessment(QA)  flags.

I studied MOD13Q1 NDVI help file  and found that possible range of values for quality layer is  (0-65534) .

For applying SG filter in TIMESAT software i need 3 ranges of QA values for assigning weight of relative importance (1, 0.5, and 0.1) to corresponding ndvi values.

Which ranges of QA values will be useful, moderate useful and useless?


texts: Sorting join layer based on attribute column?

texts:
Im joining dots with a CSV file and map my data with proportional circles. Unfortunately I have some overlapping.

How do I sort a join layer based on attribute column called aPOPCLA72 ?

Filtering virtual layer of my layer based on 1=1 order by "aPOPCLA72" desc do not work. The column is not sorted.
Sorting using MMQGIS gives me invalid sort field name error (field name: aPOPCLA72)

Why I want to sort my data: avoiding proportional circles overlapping. I want to render the big circles first.

I am using QGIS 3.6


texts: How to check if pyqgis layer is in layers panel without knowing its name

texts:
Is there an easy way to check if a layer created in pyqgis and addet to layers panel have been removed? I tried to use the layers = QgsProject.instance().mapLayersByName('my layer name'), but I don't nessesarely know the layer name, and if I use layer.name() after it hase been removed from layers panel, I get an error message.

Edit:

Error message: RuntimeError: wrapped C/C++ object of type QgsVectorLayer has been deleted

This message comes when i try to use any of the layer objects functions, but the layer is stil registerd as a QgsVectorLayer

This is for QGIS 2.18


texts: How to run cmd prompt in the background and show progress bar

texts: I am warping files using gdal. Since I have numerous files to warp how will I run the cmd in background and just have a progress bar to indicate the percentage in warping the files? Because when I tried this the cmd prompt appear consecutively just as how many files I have.


texts: tmap return an error when adding a plot layer over "World" base Map

texts:
I am working on geospatial data :

 place_name    id region_country     x     y
    <chr>      <dbl> <chr>          <dbl> <dbl>
     1 Barcelona  10102 Spain           41.4  2.18
     2 Bilbao     10103 Spain           43.2 -2.97
     3 Cadiz      10104 Spain           36.5 -6.30
     4 Gibraltar  10106 Spain           36.1 -5.35
     5 La Coruña  10107 Spain           43.4 -8.38

I want to plot this cities over the World basemap of tmap.

library(sf)
library(tmap)
data("World")

Here are my data :

x <- c(41.38333, 43.25000, 36.53361, 36.13333, 43.36667)
y <- c(2.18333, -2.96667, -6.29944, -5.35000, -8.38333)

Creating a data.frame object :

 geocode <- data.frame(x,y)

Converting to sf object :

geocode2 <- st_as_sf(geocode, coords = c("x", "y"))

Building the map :

map1 <- tm_shape(World) +  # basemap
        tm_fill() +
        tm_borders() +
        tm_shape(geocode2)+      # dots shape
        tm_dots()

map1

The error :

Error in if (!tmaptools::is_projected(shp)) { : valeur manquante là où TRUE / FALSE est requis


texts:
I'm attempting to create a tool that gives me the next available manhole number in a particular system. The manhole numbering structure looks like 1124-060E and therefore is a string field. Since inheriting, I've discontinued the letter extensions, because the sorting issues they cause. I've written a tool that has been working but is now spitting an error. The code is as follows:

def newMH(subbasin):
lst = []

with arcpy.da.SearchCursor(r"Sewer\Manholes",['OBJECTID','FACILITYID']) as cursor:
    for row in cursor:
        if subbasin.upper() in row[1].split('-')[0]:
            if row[1].split('-')[1].isdigit():
                lst.append(row[1].split('-')[1].zfill(6))
                sorted(lst, key = lambda item: (item.isdigit(), item))
                x = int(max(lst)) + 1
            else:
                pass
return str(subbasin).upper() +'-'+ str(x).zfill(4)

I'm saying if the 'subbasin' is equal to the first have in front of the dash then if the last half after the dash is a number then append it to my list. Then sort that list numerically then assign the max plus one to give me the next number in line then reconstruct the manhole format in my return line. The issue is the error I'm getting now.

"ERROR 000539: Error running expression: newMH('1124')
Traceback (most recent call last):
  File "", line 1, in
  File "", line 13, in newMH
UnboundLocalError: local variable 'x' referenced before assignment"


texts: Are two rasters the same, including null values, in ArcGIS desktop?

texts:
I am working on a task to calculate Height Above Nearest Drainage (HAND) - that is the height above the nearest point along the flow path that is on a drainage - not the geographically nearest drainage point.  I have created a raster where the value is the elevation of the next downstream point, based on a flow direction raster. I plan to populate my result raster by (unless I can figure out anything better):

Start with a raster where all stream cells are assigned a height of 0 above drainage and everything else is null. Use focal statistics to find the height above drainage for each cell's downslope cell (initially the values only along the stream) then calculate the HAND for each cell from:

HAND = cell elevation - downslope cell elevation + downslope cell HAND

This will only populate cells where downslope cell HAND is not null (i.e. adjacent to the stream initially).  Then I will add those cells to the HAND raster and iterate until I have calculated all cells in the area.

How do I know when I'm done iterating?

My current thought is to either make a copy of the HAND raster before each iteration and compare it to the raster after to see if all cells have the same value or if any cell has a value where the old raster has a null, or to count the number of non-zero cells in each raster and stop when the number is the same.  I don't know how to do either or if there is a better way.


texts: Digitizing geological information with MapInfo?

texts:
Due to vast information included in a geological map, I really need to find a way to digitize some features (mostly lines) in a rapid way.

"Click-by-click" method will cost me some millions of years in time.

Any suggestion?


texts: I've never encountered the OGR geometry type MULTILINESTRING25D before (or any of the *25D types). What are these? How do they differ from other types? I couldn't find any useful documentation from the rest of the internet.


texts:
Within the Symbol Selector > Referenced Styles

How does one edit (rename, delete, etc.) the referenced styles in ArcMap 10?

Clarification: I don't mean changing the look of the symbol using the "Edit Symbol" button. I mean editing the name or deleting one of the referenced styles.


texts: quick way to create new dissolved polygon from selected polygons (QGIS Actions)

texts:
Is there a quick way to create a new polygon from selected polygons (dissolve)?

I have 2 or more single polygons from which I want to create a new dissolved polygon within the same layer in one or two click (python action?).

Now I have to follow this steps:

Select polygons
Copy selected polygons
Paste selected polygons
merge pasted polygons and set new attributes


texts: Calculate intersecting areas for different groups of Multipolygon in QGIS

texts:
I have two polygon vector layers in QGIS: black buffers and land-use.
I want to calculate the areas of the black polygons that fall within each of the different land-use classes (these are colour coded according to a column in the attribute table).
I know how to do this if I were to convert one land-use class at a time into a polygon, but given there are 10+ land uses, I was wondering if there was a quicker way?


texts: I have actived and created tiles for layers in my Geoserver 2.5 on windows server 2008 R2. In geoserver it states that I have used 583MB in total after seeding the layer in GWC. Strangely I cannot see where the tiles are being created. I have not set an environment variable to move the tile directory anywhere else, so I expected the tiles to be in the data_dir. However the only directory I have is gwc-layers which contains the xml configuration files. So where is the tile cache?


texts: Converting decimal degrees units to km in R?

texts:
I am trying to figure out how I can transform my unit of length from decimal degree to km when I import the following simple ArcGIS polygon shape.  Ultimately, I would like to have this information available when I assign this spatial polygon to my window in spatstat [i.e., Window area = xx square unit & Unit of length: xx km].

I did came across this post: R-sig-geo

This is what I've done so far, but I don't understand why/if I need to do so. Doesn't my div0A SpatialPolygonsDataFrame already has my CRS argument defined? Why do I need to spTransform then?

library(spatstat)
library(maptools)
library(lattice)
library(sp)
library(foreign)
require("rgdal")
div0A <- readOGR(dsn="Archive", layer="Projections")
> summary(div0A)
Object of class SpatialPolygonsDataFrame
Coordinates:
     min       max
x -80.00 -57.63430
y  66.25  78.16666
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
Data attributes:
       Id      NPAzimutha         UTM20            UTM19           AlberEA
 Min.   :0   Min.   :348877   Min.   :349232   Min.   :349162   Min.   :348656
 1st Qu.:0   1st Qu.:348877   1st Qu.:349232   1st Qu.:349162   1st Qu.:348656
 Median :0   Median :348877   Median :349232   Median :349162   Median :348656
 Mean   :0   Mean   :348877   Mean   :349232   Mean   :349162   Mean   :348656
 3rd Qu.:0   3rd Qu.:348877   3rd Qu.:349232   3rd Qu.:349162   3rd Qu.:348656
 Max.   :0   Max.   :348877   Max.   :349232   Max.   :349162   Max.   :348656

proj4string(div0A) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
xxlcc <- spTransform(div0A, CRS("+init=epsg:4326"))
res1 <- lapply(slot(xxlcc, "polygons"), function(x)
  sapply(slot(x, "Polygons"), slot, "area"))

> summary(unlist(res1))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  89.72   89.72   89.72   89.72   89.72   89.72

newW<-as(xxlcc, "owin")
> summary(newW)
Window: polygonal boundary
single connected closed polygon with 29 vertices
enclosing rectangle: [-80, -57.6343]x[66.25, 78.16666]units
Window area =  89.7204 square units

Does this look correct? How can I specify/know my unit of length in km?
Moreover, this window is ~Baffin Bay, west of Greenland and I am wondering if I should be using a different datum? E.g. EPSG Projection 4747 - GR96? Although, I will have to create ppp object, for which the coordinates have been taken with GPS satellite navigation system.

I would appreciate if someone can walk me through this or suggest an online example or further readings. I read chapter 4: Spatial Data Import and Export in Applied Spatial Data Analysis with R (2008), but I'm obviously very new to all of this.


texts: Process to load elevation data into Oracle 11G

texts:
We have the National Elevation data set from USGS and considering loading it into Oracle (11G) database.
Can someone help with the process/steps?


texts: Connecting/importing PostGIS database in ArcGIS Desktop via localhost?

texts:
I want to import PostgreSQL database in PostGIS to ArcGIS Desktop?

Is it required to install ArcSDE?

I am using ArcGIS 10.3.


texts: How to merge a dataset with geometry an another with number data?

texts:
I want to merge a dataset with geometry and another with the number data to use.

In each one datasets are a field (codparr in one, and cod_parro in another dataset) with 51 values(3601701,3601702,3601702,...) in each one that are repeated in the other, to use this field as (2nd file / key field) merge...

but I can´t merge datasets.

I attach an capture

 ( snap.... ) where you can see that I have not the option to MERGE

Further more, after click in each dataset you can see that there no exist the option of MERGE neither in the dataset with numerical data

nor in the dataset with the geometry

INSTEAD......

a friend of mine, show me that he has a different screen (with the data view / map view in the upper screen)

 and have the option, as you can see in capture YES.... whose is similar/equal to the tutorials/guides captures/videos.


texts: I want to add an "World Image" Datastore to my GeoServer, based on a georeferenced JPEG that I have on my local drive. Can I upload it to the GeoServer data folder from within the GeoServer application?


texts: Exporting all pictures in Feature class as Raster Field in ArcPy?

texts:
I have received file geodatabase (FGDB) and it has images added in raster field type (in attribute table). I wanted to extract these images using ArcPy into the folder.
Does someone have a solution?
I have observed that if we put images in the raster field managed by GDB, the size of GDB increases so much, as per my understanding, it converts images into tables as per pixel value.


texts: arcpy.Near_analysis throws error: feature class does not exist or is not supported or cannot add field NEAR_FID

texts:
I am trying to determine the closest line feature to features in a point feature class by using the arcpy.Near_analysis function, as such:

search_radius = '' #optional
location = 'LOCATION'
angle = 'NO_ANGLE'
method = '' #optional
arcpy.Near_analysis(points, lines, search_radius, location, angle, method)

This in a Python window within ArcGIS Pro returns the error:

Traceback (most recent call last):   File "<string>", line 1, in <module>   File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 1456, in Near
    raise e   File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 1453, in Near
    retval = convertArcObjectToPythonObject(gp.Near_analysis(*gp_fixargs((in_features, near_features, search_radius, location, angle, method), True)))   File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 496, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True)) arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Features: Dataset db_gdb.DBO.GPS_Point does not exist or is not supported ERROR 000732: Near Features: Dataset db_gdb.DBO.GPS_Lines does not exist or is not supported Failed to execute (Near).

both 'points' and 'lines' feature classes exist, and other arcpy functions can perform on these exact variables, and I see nothing about them that would prevent them from being supported by this function.

A stand-alone python script(IDLE) throws this error:

Traceback (most recent call last):   File "C:\Users\me\Desktop\file.py", line 40, in <module>
    arcpy.Near_analysis(points, lines, search_radius, location, angle, method)   File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py", line 1456, in Near
    raise e   File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\analysis.py", line 1453, in Near
    retval = convertArcObjectToPythonObject(gp.Near_analysis(*gp_fixargs((in_features, near_features, search_radius, location, angle, method), True)))   File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 496, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True)) arcgisscripting.ExecuteError: ERROR 000852: Cannot add field NEAR_FID to C:\Users\me\Documents\ArcGIS\Projects\GPS_Point Failed to execute (Near).

I'm not sure why the two different environments throw different error messages, but I am following the documented guidelines set forth by Esri, so am not sure what the issue is.  I have an Advanced license and I see nothing in the documentation that would lead me to think there is something missing in either of these feature classes to prevent this function from working.


texts: Import qgis.core to Python IDLE - No module named 'qgis'

texts:
I'm trying to import the qgis.core module as explained here: https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/intro.html

I'm working on a macOS system and installed python.org Python 3.6. Then I installed qgis.org QGIS 3 and did the following steps:

Check QGIS install path:

  QgsApplication.prefixPath()
  '/Applications/QGIS3.app/Contents/MacOS'

Change PYTHONPATH, PATH, LD_LIBRARY_PATH to:

  export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/Applications/QGIS.app/Contents/MacOS/lib/:/Applications/QGIS.app/Contents/Frameworks/

export PYTHONPATH=$PYTHONPATH:/Applications/QGIS.app/Contents/Resources/python/

export PATH=$PATH:/Library/Frameworks/GDAL.framework/Programs/

Write this code into Python IDLE:

  from qgis.core import *

Then I get this error message:

  Traceback (most recent call last):
    File "<pyshell#1>", line 1, in <module>
      from qgis.core import *
  ModuleNotFoundError: No module named 'qgis'

Could anyone help me?


texts: Multiple ring buffer on point-shape error

texts:
I need to calculate multiple ring-buffers, but every time I try, arcgis produces a failure:

Executing (Multiple Ring Buffer): MultipleRingBuffer Y:\GIS\RDY\clip\sn_clip.shp Y:\GIS\RDY\sn_clip_ring.shp 100;250;500;750;1000;1250;1500;1750;2000;2250;2500;2750;3000;3250;3500;3750;4000;4250;4500;4750;5000;5250;5500;5750;6000;6250;6500;6750;7000;7250;7500;7750;8000;8250;8500;8750;9000;9250;9500;9750;10000;10250;10500;10750;1100 Meters distance ALL FULL
Start Time: Fri Feb 15 12:03:02 2013
Running script MultipleRingBuffer...
ERROR 999999: Error executing function.
The table was not found.
The table was not found. [union]
The table was not found.
The table was not found. [union]
Invalid Topology [Topoengine error.]
Failed to execute (Union).

Completed script MultipleRingBuffer...
Failed to execute (Multiple Ring Buffer).
Failed at Fri Feb 15 13:58:58 2013 (Elapsed Time: 1 hours 55 minutes 56 seconds)

I don't get it. Can anyone help me with this? Note that I'm not used to python.


texts:
How to store jpeg file in spatialite blob column so that spatialite gui and Qgis can recognise the format and show the image.
I have created in memory jpeg from raster bands using jpeglib and stored that in spatialite but spatialite gui can not recognise the format.
I also tried converting the whole jpeg to hexwkb format no results.

Thanks


texts: Showing values in attribute table in scientific (exponential) notation in QGIS3

texts:
One of my virtual field contains too big numbers like "5994865699857.456". I need to show those values in scientific/exponential mode like 5.99E+12. I didn't find proper function in Field Calculator. I searched for solution among others questions but all questions address inverse problem - converting scientific view to number one, and I didn't find a solution among answers.
Is there some trick in Field Calculator to solve the problem?
Edit: I'd like to see value in scientific notation in Feature Attributes form when I use "Identify Features" function.
Edit2: I also tried to create custom function:
from qgis.core import *
from qgis.gui import *

@qgsfunction(args='auto', group='Custom')
def fd_sci(value1, feature, parent):
    vtostring = str(value1)
    a = "{:.2e}".format(vtostring)
    return a

Then created a new virtual field using expression:
fd_sci("virtual_field_with_big_number")
But the result is NULL


texts: ArcGIS mobile to use external GPS device

texts:
Im currently tasked to create an app that would sketch a GPS Device's movement and plot it on the ArcGIS map. By default, the ArcGIS locationDisplay class calls the IOS class CLLocation to retrieve the GPS data in which it only uses the internal GPS of the iPhone.

My question is: Is there any way to force the iPhone to use the GPS of an external device? Or can we configure the ArcGIS to use the data from the external device?


texts: Changing labels from external CSV in QGIS

texts:
I got my hands on a TopoJSON map that unfortunately uses a mix of names for the polygons, sometimes the official name, sometimes not (English). I do, however, have the ISO3 code, which is great.
I also do have a dataset with the official names I would like to use in this TopoJSON file. Is it possible to replace a specific attribute by matching it with a CSV?


texts:
I'm trying to use osmosis to update my data ( http://wiki.openstreetmap.org/wiki/Minutely_Mapnik ) but I'm facing the follow error:

org.openstreetmap.osmosis.core.OsmosisRuntimeException: Unable to read the state from the server.
        at org.openstreetmap.osmosis.replication.common.ServerStateReader.getServerState(ServerStateReader.java:116)
        at org.openstreetmap.osmosis.replication.common.ServerStateReader.getServerState(ServerStateReader.java:50)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.runImpl(BaseReplicationDownloader.java:290)
        at org.openstreetmap.osmosis.replication.v0_6.BaseReplicationDownloader.run(BaseReplicationDownloader.java:383)
        at java.lang.Thread.run(Thread.java:745)
Caused by: java.io.IOException: Server returned HTTP response code: 407 for URL: http://planet.openstreetmap.org/replication/hour/state.txt
        at sun.net.www.protocol.http.HttpURLConnection.getInputStream0(HttpURLConnection.java:1840)
        at sun.net.www.protocol.http.HttpURLConnection.getInputStream(HttpURLConnection.java:1441)
        at org.openstreetmap.osmosis.replication.common.ServerStateReader.getServerState(ServerStateReader.java:97)
        ... 4 more

jul 28, 2016 8:26:55 AM org.openstreetmap.osmosis.core.Osmosis main
GRAVE: Execution aborted.
org.openstreetmap.osmosis.core.OsmosisRuntimeException: One or more tasks failed.
        at org.openstreetmap.osmosis.core.pipeline.common.Pipeline.waitForCompletion(Pipeline.java:146)
        at org.openstreetmap.osmosis.core.Osmosis.run(Osmosis.java:92)
        at org.openstreetmap.osmosis.core.Osmosis.main(Osmosis.java:37)
        at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.lang.reflect.Method.invoke(Method.java:497)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launchStandard(Launcher.java:330)
        at org.codehaus.plexus.classworlds.launcher.Launcher.launch(Launcher.java:238)
        at org.codehaus.plexus.classworlds.launcher.Launcher.mainWithExitCode(Launcher.java:415)
        at org.codehaus.plexus.classworlds.launcher.Launcher.main(Launcher.java:356)
        at org.codehaus.classworlds.Launcher.main(Launcher.java:47)

I'm already done export JAVACMD_OPTIONS="-Dhttp.proxyHost=MY_PROXY -Dhttp.proxyPort=6060 -Dhttp.proxyUser=MY_USER-Dhttp.proxyPassword=MY_PASSWORD" but with no success. I'm able to access this URL using the browser.

I know... 407 is for Authentication. The setting I'm using is correct for sure or my browser will give same error.

One more information: as seen in (How to configure Geoserver to use a proxy when accessing external services?) I'm facing same problem in Geoserver. SOME external WMS services gives me same 407 error, others don't.

I think this is because of proxy auth implementation inside the code (osmosis and geoserver may be using same method): https://stackoverflow.com/questions/17183277/http-error-407-proxy-authentication-required


texts: Adding band value as a property in Google Earth Engine

texts:
I have an image collection of different anomaly values for a specific time period. I'm simply looking to count how many of these values are over 1 in my time period.

This is my code

// create collection
 var collection = ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')
  .select('sst')
  .map(function(image){return image.clip(EEZ)});

// get mean value over ten years
var reference = collection.filterDate('2002-01-01', '2012-12-31')
.sort('system:time_start', false);

//calculate anomaly for 2018
 var series = collection.filterDate('2018-01-01', '2018-12-31').map(function(image) {
return image.subtract(mean)
.set('system:time_start', image.get('system:time_start'));

I thought there may have been a way of adding the band value to a property for each image and then filtering the metadata of this collection for values 1 or over, but I can't seem to find a way of doing it.

Are there other ways of filtering/removing data values below specific targets?


texts:
What is the best solution so far for printing Google Maps as a base layer of Vector layers in Geoserver/OpenLayers application ?

I used to print WMS layers using GeoExt print page example, but when I tried it with Google Maps, I got a blank file !


texts: Populating segment to and from columns with roadnames without manually typing each name in ArcMap?

texts:
I have been given three shapefiles. One is a road centerline file with roadnames containing all the roads in a town (55k+ edges). One is an original file with 9316 edges that contain to and from fields indicating the intersecting road name. The last is a new shapefile that was created containing 1300 edited edges that have been either merged, split or some combination. The to and from fields still contain the original to and from field data, which is now incorrect.

Is there any way to populate to and from field for the new edges without having to manually enter each roadname?


texts:
Do you know if there is a way to block a user after 3 attempts of login?

I was looking for information about but I haven't find anything.


texts: Iterators in model builder - merging nested files with same filename ending

texts:
I am trying to merge files that are nested in different folders by their file names. For example folder 1 contains up to 100 subfolders that each contain a set of shapefiles: subfolder1_A.shp, subfolder1_B.shp, subfolder1_C.shp, etc.

I want to merge shapefiles from each subfolder with the same file name ending (e.g. merge subfolder1_A, subfolder2_A, subfolder3_A, ..., subfolder100_A), so the final output is a series of merged shapefiles with the same file name ending (A, B, C, etc.). So for example, A_merged.shp, B_merged.shp, C_merged.shp, etc.

I am not familiar with python, and would prefer to do this using model builder. I have tried using iterate feature class and merge tools, but not sure how to make this work like I want. Possibly using inline substitution or the iterator wildcard somehow?

Any idea what I should do?


texts: How to write an XML file to config a layer (GeoToolsLayer) fetched from GeoServer?

texts:
I need to connect to GeoServer, fetch the layer and display it in Geomajas. I have followed the instruction in this link http://www.geomajas.org/node/280 but it doesn't work for me.

The map is just blank. I don't find any compilation error or runtime exceptions in Jetty console, and there is also nothing in geoserver.log on the server. Therefore, I have no clue to figure out what the exact problem is.

For other kinds of layer (from shapefile, Google Map...), it works OK.

Where can I find a working XML file to config such layer?
How can I debug or find out the reason why the map is just blank?


texts: Using Near function on multiple variables at same time?

texts:
I am trying to calculate between variables and a Fishnet using the Near Analysis Tool.
I have about 50 variables (Binary (Is there "a" information at this location)) for each XY points (300000) and I am trying to calculate the distance between them and a fishnet. (e.g. What distance is variable "a" to the closest square or points of the fishnet).

What I am currently doing is: Separate each of these variables and then calculate them separately on a specific Fishnet (Variable "a" with Fishnet "a") then merging them together.

Is there a way to do that more efficiently?


texts: Using Configurable Map Viewer Attributes Table?

texts:
I've created a CMV query widget that I am trying to link up with the Attribute Table  widget. The table tab is being created in the bottom pane but that's as far as I'm able to get. I'm unable to get the query to actually execute and I'm unsure of the reason. I have used the Search widget as my starting point.

In my executing function I have the following:

var query = {
    "layerID": "SelectOverlayData",
    "sublayerID": 1,
    "type": "spatial",
    "outFields": [
        "TMK",
        "TPL_ACRES",
        "owner",
        "B_al02"
    ],
    "outputSpatialReference": {
        "wkid": 102100,
        "latestWkid": 3857
    },
    "spatialRelationship": "esriSpatialRelIntersects",
    "geometry": null,
    "where": "TPL_ACRES <= 40 AND B_al02 = 'Yes'"
};

and

var grid = {
    "minWidth": 50,
    "columns": [
        {
            "field": "TMK",
            "label": "Parcel ID",
            "width": 100
        },
        {
            "field": "TPL_ACRES",
            "label": "Acres",
            "width": 100
        },
        {
            "field": "owner",
            "label": "Owner",
            "width": 100
        },
        {
            "field": "B_al02",
            "label": ">>Protect farmlands",
            "width": 100
        }
    ],
    "sort": [
        {
            "attribute": "TMK",
            "descending": "ASC"
        }
    ],
    "defaultColumnWidth": 100,
    "columnHide": true,
    "columnReorder": true,
    "pagination": true,
    "paginationOptions": {
        "rowPerPage": 50,
        "previousNextArrows": true,
        "pagingLinks": 2,
        "pagingTextBox": true,
        "showLoadingMessage": true
    }
};

Into

        topic.publish(this.attributesContainerID + '/addTable', {
            title: layer.name,
            topicID: this.topicID,
            queryParameters: query,
            gridOptions: grid,
            idProperty: layer.idProperty
        });

The tab is created in the bottom pane with the rest of the Attribute Table widget but the query is never run (nothing in the Network tab of Chrome) and no records or columns ever appear in the table. No js errors. I must be missing something small.

Any help?


texts: Joining point data to nearest polyline with additional attribute criterion in ArcMap?

texts:
I am trying to join point data to polylines in ArcMap 10.3.1, for which I have used the Near tool, and that gives me the distance and FIDs of the nearest polylines, and I can then use the FID to join the two datasets.
However, I would like to use an additional selection criterion that is recorded in the point data to choose the polyline to link to. That is, the nearest polyline is not necessarily the adequate line to join to, but it might be a polyline that is further away, if it has the matching attribute I am looking for.
I added a screenshot of the problem: the yellow point feature needs to be joined to one of the nearby line features, which have different attributes according to their colors. I could of course select the proper line manually, but I have about 5000 points to join, so that's not really efficient...
The selection attribute is the line material, which is recorded in both the point and the line features.
Can anyone help?


texts: QGIS 3.12 - added CSV file - only showed 5 points on map

texts:
Already saw the same question 2016 but no answer.
I have loaded my file in QGIS - the layer showed the file and attribute table showed a complete dataset but only 5 points are shown in the map. i made sure that there are no overlapping points.
I'm using Windows 10.
My process:
-saved file as .csv ---> opened QGIS 3.12 --->add layer ---> add delimited text layer ---> file name and layer name automatically filled in ---> encoding is UTF-8 ---> file format: custom delimiters (used tab and comma) ---> record and field options (did not do anything) ---> geometry definition (used point coordinates, filled in x and y fields and used WGS84) --->clicked add


texts: Extracting Elevation (DEM) from air photos using open source software?

texts:
I have a few pairs of stereo air photos and I want to extract a DEM from them.

I think that is possible using PCI (OrthoEngine), but I would like to do it using Open Source (or at least free) software.

I've looking in GRASS GIS and Orfeo Toolbox but I didn't find how to do it?
(It is possible that I have not looked correctly...).

Is it possible?


texts: Downloading the Geoserver SQL Server Extension

texts:
I have successfully installed GeoServer on a 64 bit Windows Server and would like to install the Geoserver SQL Server Extension so that I can hook GeoServer up to SQL Server 2017.

The GeoServer instructions at http://docs.geoserver.org/stable/en/user/data/database/sqlserver.html say that I can download the extension from http://geoserver.org/download/ but I can't see the extension listed there.

Does anyone know from where I can download the extension?


texts: ArcPy Is error 000210 occuring due to permissions

texts:
I have a script that worked fine a month ago but I suspect due to some reconfiguring of my citrix profile by IT as part of a maintennance drill they changed something that affects my ability to create folders and file gdbs from a script. In this script, a for loop is used to create multiple clips of a three large line feature classes into smaller areas. I want the script to create a folder for each location and a generic Features.gdb within that folder and I thought this happens as part of the clip tool being called in the script.

Here is my code:

outputFolder = r"N:\Common\MapClips\CO_181"
env.workspace = r"N:\Workspace\Updates_2017\gisdb\HERE_2017_Q1_3530_Recology2017Q1.gdb"
env.overwriteOutput = True

Cutter = r"N:\Common\MapClips\CO_181\Features.gdb\CO_181_Dissolve"
Dough1 = r"N:\Workspace\Updates_2017\gisdb\HERE_2017_Q1_3530_Rec2017Q1.gdb\here_2017_q1_streets"
Dough2 = r"N:\Workspace\Updates_2017\gisdb\HERE_2017_Q1_3530_Rec2017Q1.gdb\here_2017_q1_streets_alias"
Dough3 = r"N:\Workspace\Updates_2017\gisdb\HERE_2017_Q1_3530_Rec2017Q1.gdb\here_2017_q1_streets_parity"

fcs = [Dough1, Dough2, Dough3]

with da.SearchCursor(Cutter, ["SHAPE@", "Fran"]) as rows:
    for row in rows:
        #1. create fgdb in outputFolder variable
        outGDB = os.path.join(outputFolder,"CO_" + row[1], "Features.gdb")
        geom = row[0]# arcpy.CopyFeatures_management(geom, r"in_memeory\temp")
        for fc in fcs:
            out_feature_class = os.path.join(outputFolder, outGDB, row[1]+ "_" + os.path.basename(fc))
            print out_feature_class
            arcpy.Clip_analysis(fc,geom,out_feature_class)
            del fc
        del row
        print "Created " + out_feature_class

This throws the following error:

 Traceback (most recent call last):
      File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
        debugger.run(codeObject, __main__.__dict__, start_stepping=0)
      File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
        _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
      File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
        exec cmd in globals, locals
      File "N:\Python\Completed scripts\Ratto_Cal_Street_Update_Prep.py", line 3, in <module>
        import arcpy
      File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\analysis.py", line 56, in Clip
        raise e
    ExecuteError: ERROR 000210: Cannot create output N:\Common\MapClips\CO_181\CO_Nov\Features.gdb\Nov_here_2017_q1_streets
    Failed to execute (Clip).

I went and created the first output location manually:

N:\Common\MapClips\CO_181\CO_Nov\Features.gdb

And it worked in the first iteration and then bombed when it got to the next clip feature as there was not yet a folder or gdb created for it yet.

I know this worked fine for me before and I didn't have to manually create the folders and fgdbs. I can do that if necessary but I would prefer to have this process done automatically through the script.


texts: Seeking shapefiles of global climate/climatic zones

texts:
Have been looking around a bit for a shapefile/kml/etc file of global climatic zones (as listed here: https://en.wikipedia.org/wiki/Geographical_zone) -- temperature, subtropical, tropical, etcetera.

Also would be very pleased to find gis overlays detailing photoperiod areas, that sort of thing?

An online GIS service would be fine for this as well for this, as long as its a real GIS with aerial overlays etc, not a plane paper PDF deal, which is about all I can find.

One promising (but dead link) source I found before posting question: https://www.esrl.noaa.gov/psd/data/usclimdivs/boundaries.html


texts: Is it possible to read in specific bands from a multi-band raster with gdal or rasterio?

texts:
I wasn't able to find an answer to this after searching the help files and the forums.

I have a (virtual) raster with 60 bands, where each band represents a realization in time (think of frames of a video). I want to be able to load in cropped, specific band sets from this virtual raster without loading all of them--something like

myVRT = gdal.Open(vrtfile)
bands = myVRT.GetRasterBand([2,4,6,14,35])
cropped_bandstack = bands.ReadAsArray(
   xoff = xoffset,
   yoff = yoffset,
   xsize = xsz,
   ysize = ysz)

which would create an array with only bands 2,4,6,14, and 35.

I haven't been able to find a way to do this without loading each band separately and filling in a 3D matrix.

Another alternative that I also want to avoid is to re-create each of the individual cropped bands (which are also virtual), then create a new virtual raster stack using only the bands I want, and finally read it into an array.

I saw this was listed as a feature to add to rasterio but I couldn't find any evidence of it actually being added.

Edit: After receiving an answer, I think I should clarify that I'm asking if this can be performed in a single call, rather than implementing a loop to read each band individually.


texts: Creating Contours from Raster using PyQGIS

texts:
I am trying to create contours from a raster layer I have loaded in my QGIS project as below using gdal
from qgis.core import QgsProject
import processing

## QgsVectorLayer('Coordinate','Layer name','memory layer')
layer = QgsVectorLayer('LineString?crs=epsg:28356', 'Contour' , 'memory')

for raster in QgsProject.instance().mapLayers().values():
    if 'h_Max' in raster.name():
        processing.run('gdal:contour', raster, 1, 'ELEV' , layer)

I've never had any chance of getting the code to work, the error message I got is
AttributeError: 'QgsVectorLayer' object has no attribute 'feedback'

I am not sure when I should define this feedback attribute, any suggestion?


texts: QGIS v3.10 GRASS plugin on MacOS for hydrology

texts:
Installed QGIS 3.4 and then 3.10, both gave errors "Grass was not found in /usr/local/opt/osgeo-grass/grass-base" In both cases the directory /usr/local/opt/osgeo-grass/ did not actually exist. I tried manually installing grass both with brew install and Grass 7.7 from grass-7.7dev.dmg.zip from the website. I tried setting the path for grass in QGIS to /Applications/GRASS-7.7.app and a few other locations that I forgot now.

All I need to do is create a flow accumulation raster from a DEM I have. I tried using the SAGA plugin in QGIS and that failed, wasn't able to generate a raster. I feel like my plugins have not been compiled correctly. However the base QGIS raster analysis works: I could create slope, aspect rasters sucessfully. The contents of my /usr/local/opt are currently

ant            lcms2           osgeo-libgeotiff
ant@1.10        libcryptopp     osgeo-libkml
armadillo       libdap          osgeo-libspatialite
arpack          libev           osgeo-netcdf
boost           libevent        osgeo-postgresql
boost@1.71      libffi          osgeo-proj
brotli          libiconv        pcre
c-ares          libidn          pcre1
cfitsio         libidn2         perl
cgal            libjpeg         perl@5.30
cryptopp        libjpeg-turbo       pkg-config
curl-openssl        libjpg          pkgconfig
eigen           libkml-dev      popt
epsilon         libmetalink     python
expat           libmpc          python2
fftw            libpng          python3
freexl          libpq           python@2
gcc         libssh2         python@3
gcc@9           libtiff         pzstd
gdal2           libtool         qhull
gdbm            libunistring        readline
geos            libxml2         rtmpdump
gettext         libzip          sfcgal
giflib          little-cms2     sqlite
glib            mdbtools        sqlite3
gmp         mpfr            swig
hdf4            nghttp2         swig@4
hdf5            numpy           swig@4.0
hdf5@1.10       numpy@1.17      szip
hwloc           open-mpi        tcl-tk
icu4c           openblas        unixodbc
isl         openjpeg        webp
jansson         openldap        wget
jasper          openmpi         xerces-c
jemalloc        openssl         xz
jpeg            openssl@1.0     zlib
jpeg-turbo      openssl@1.1     zstd
json-c          osgeo-gdal
krb5            osgeo-hdf4


texts: ArcMap Python Question: Trying to Field Calculate Slope Percent in to Text

texts:
I have a layer of polygons that I spatially joined to a raster to point layer of slope percents. My polygons now have the percent of slope in a field called Avg_point. I want to use the field calculator to take that percent number (something like 22.345 or 56.45 etc.) and turn it into a text in the "slope" field (that was already in the layer of polygons) and I want that text to say something likes "Steep (>25%)"

So I wrote a script

SLOPE=slopecorrect(!SLOPE!, !Avg_grid_code!)

PreLogic:

def slopecorrect(slope, avg):
    if avg > 25:
       return Steep

I am getting no errors but my "slope" field just fills in  I have tried other things but I can't seem to get it to fill in that field

I have looked around a bit and thought this thread was helpful but its still not producing results
ArcMap Calculate Field..if "text" in !Field!


texts: Interpolate value at a specific point using kriging

texts:
Is it possible to interpolate the specific value at a point using the Kriging tool in the spatial analyst toolbox? If not, then, how can this be done in Arcgis?

To my knowledge, Arcgis only produces a raster. This isn't applicable for me because I need to interpolate something quickly and it has to be at a specific point (geographic coordinates). I'm trying to interpolate weather data (temperature, humidity etc).


texts: Selecting multiple variables by Attribute from different target layers

texts:
I have a model with species point data as a source layer and a number of different shapefiles as target layers. To start I run a search by location tool to determine where the source and target layers intersect, and then copy features to create a new target layer. Next I want to use the select by attribute tool to determine where in my original target layer are common attributes to my new target layer.
My question is: Is there an SQL expression that would allow the search by attribute tool to select from an individual field in the new target layer and compare it to the original layer?


texts: Performing a field calculation affects selectable layers in ArcGIS 10.1

texts: In ArcMap 10.1 I have a table open for a layer that I have set as not selectable. Yet when I use the Field Calculator for a record(s) in that table, it automatically makes the associated layer a selectable layer, which I do not want, since I am in the process of editing another layer and only want that layer selectable. Is there a way to keep this from happening?


texts: Extract a polygonal area from TIFF Image

texts:
I have a TIFF image (actually a satellite image) which is consisting of 5 channels (R G B IR NIR) for a specific region and I have a KMZ file which defines a polygon within that specific region.

I want to extract/map that area from TIFF file and want to read channels pixel-by-pixel.

For example imagine that satellite image encapsulates the Germany and KMZ file has the coordinate information of Munich. What I have been trying is to mask/or cut the region defined by KMZ and then be able to read channels separately only for specified region (Munich).

I have read on a library called GDAL however I couldn't find any appropriate solution.


texts: Angular rotation using Proj4, for local coordinate system

texts:
I have a polygon feature on map to which I need to transform the coordinates to local coordinate system with respect to new origin while keeping the same datum.

So far using the Proj4 library and by defining the custom map origin I could able to re-project/transform my global coordinates to local as follows

+proj=etmerc +lat_0=52.481392 +lon_0=13.355795 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs --Using this Proj4 definition

Here in the image you can see the marker with name "New Origin" is the one I want to keep at center and transform the coordinates of highlighted polygons to local coordinate system

These are the original coordinates of selected polygons

These are the transformed coordinates of same polygon w.r.t. new map origin

What I am looking forward now is the case of rotated map origin, as of now the current map origin is facing northwards i.e. to north pole, which I need to be rotated by around 90 degrees. Ideally this rotation should effect coordinates of the polygon (as those -ve now may get positive and vice versa).

+proj=omerc +lat_0=52.481392 +lonc=13.355795 +alpha=90.0 +k=1.0 +x_0=0 +y_0=0 +no_uoff +gamma=90.0 +ellps=WGS84 +units=m +no_defs --New definition I am trying

Regarding the same, I found this question helpful and thought to try the Oblique Mercator projection, however the results are not changing when I did the 90 degree rotation for alpha and gamma


texts: Using gdal to convert from UTM to equidistant lat/lon

texts:
I am getting rather confused with converting a UTM gridded dataset into a equidistant/rectangular latlon gridded dataset. Can anybody help me?

I have a UTM set with the following projection:

PROJCS["WGS 84 / UTM zone 18N",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-75],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","32618"]]
Origin = (556950.000000000000000,5573990.000000000000000)
Pixel Size = (90.000000000000000,-90.000000000000000)
Metadata:
 AREA_OR_POINT=Area
Image Structure Metadata:
 INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  556950.000, 5573990.000) ( 74d12' 0.38"W, 50d18'54.95"N)
Lower Left  (  556950.000, 4905470.000) ( 74d17' 9.82"W, 44d18' 0.91"N)
Upper Right (  655950.000, 5573990.000) ( 72d48'37.20"W, 50d17'50.61"N)
Lower Right (  655950.000, 4905470.000) ( 73d 2'43.67"W, 44d17' 8.74"N)
Center      (  606450.000, 5239730.000) ( 73d35'30.62"W, 47d18' 8.60"N)
Band 1 Block=1100x1 Type=Float32, ColorInterp=Gray

And want to convert this to a rectangular latlon set (with corners just as they are provided in the corner coordinates. What to do? I am sure it should be easy...The biggest problem is that if I directly convert this into latlon, it won't be equidistant...


texts: connection to postgresql; From ArcGIS for Desktop 10.1 OR ArcGIS for Server 10.1

texts:
As it is mentioned in ArcGIS Resources "PostgreSQL client libraries are required on the ArcGIS client application from which you will connect to the database".

My questions are:

by ArcGIS client applications, they mean e.g. ArcGIS for Desktop version 10.1?
Can ArcGIS for server 10.1 have a direct connection to e.g. Postgresql? and if NO why there is pg_client_linux64?

Are we supposed to put the client libraries in both bin directories of ArcGIS for Desktop and ArcGIS for server (on Centos)?


texts: How to get legend image from a WMS using PyQgis

texts:
I want to load a WMS with PyQGIS. I know this is possible on version 2.18
If it is done with the graphical Interface by checking on the Use contextual WMS legend

Now I manage to load a WMS with PyQGIS but I can not figure how to load the WMS and its legend.

Below there is partial working code snippet:

urlWithParams = "url=http://www.geoportaligm.gob.ec/50k/wms&format=image/png&layers=igm:acequia_l&styles=&crs=EPSG:32717"
rlayer = QgsRasterLayer(urlWithParams, 'testing', 'wms')
rlayer.isValid()
QgsMapLayerRegistry.instance().addMapLayer(rlayer)


texts:
I use ArcMap 10.3.

How do I create an expression in RasterCalculator that allows me to sum up 4 rasters (every pixel has an associated value) and the result to be 0 if any of the 4 values is 0?

For example: 1+2+0+1 must be equal to 0 but 1+2+1+1 must be equal to 5


texts: Copy .prj file into QGIS projection folder

texts:
I want a custom projection used in ArcGIS be available in QGIS. I have the .prj file copied out of a ArcGIS folder (containing different projections). Now, I want this file to be added into QGIS. Is there a folder where I can copy this file into?

The custom projection is:

PROJCS["ETRS_1989_UTM_Zone_33N7",
GEOGCS["GCS_ETRS_1989",
DATUM["D_ETRS_1989",
SPHEROID["GRS_1980",6378137.0,298.257222101]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",3500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",15.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0]]


texts: osmpgsql fails to load all data from OSM file

texts:
I have a project where I need to use osm2pgsql to load an OpenStreetMap file (osm.pbf).
(Using a command like this:
osm2pgsql -c -d database_name -U user -H localhost -W -S "path to style file" "path to osm-file")

The command runs with no errors, but whenever I use it for larger areas, data on relations are not imported (the table for relations is created in Postgres, but it is empty). Same thing for the road table. It works fine for smaller areas.
Any ideas as to what the issue might be and, how I can fix it?


texts:
I just can't find the size changing method for the graduated style.

I have an integer containing column, and I would like to display it with icons which are bigger if the value is higher.

So, is this a bug in QGis v2.14.0 or did I miss some config steps?


texts: Cannot edit location of points in Feature Layer in ArcGIS Online after using Collector

texts:
I have students in my GIS class who have created feature layers in ArcGIS Pro to collect points on our campus using Collector on a cell phone.  After going out to collect points they return to the classroom to edit those points.

While they are able to update attributes on ArcGIS Online they are unable to relocate any points without them reverting back to the original location.

Is there something we can do in either building the feature layer or adjusting its settings once it is published that could alleviate this issue?


texts:
I am teaching myself GIS at the moment, and there are still some matters I haven't quite got the hang of.

I'm using ArcGIS 10.1, and have three datasets:

LAYER 1: A population dataset at a high resolution, as a tiff raster. (Edit: the convert to polygon tool does not work on this raster)
LAYER 2: A polygon layer, created as a fishnet, of equal squares, defined by two vectors.
LAYER 3: An outline of the coastline.

What I need to do is to calculate the total population within a buffer of distance X km from the coastline for each polygon defined by the fishnet, that falls on the coastline.

I am currently struggling to combine the population data with that of the polygon layer however - any ideas?


texts:
We have been happily using the 2010 Census Demographic Profile data in Shapefile format for some time. We now want to start using the ACS 5-year Estimates Data Profiles which, sadly, are only available in Geodatabase formats.

What's the best solution?

Do I bite the bullet and buy some version of ArcGIS, are there any services that can convert the files for me or is there some conversion program I can purchase?

Note that I am not a GIS specialist. I've learned enough to do what I want with the 2010 Census demographic shapefiles, but this is not the main focus of my website.

I seek recommendations and advice.


texts: How to add dynamic tree and layer in GeoExt

texts:
I am using GeoExt, OpenLayers and GeoServer for map. Now I want to create dynamic tree and layer using GeoExt. So that any user can create a tree and add its layer (child) to that tree dynamically. I have seen many examples on GeoExt website tree panel like http://api.geoext.org/1.1/examples/tree.html, but this  a static example.

I want it so that a user can add any number of tree and it layer dynamically. It is possible by GeoExt and OpenLayers, If there is any other solution for it please suggest.

How can I do it?

Please help me.


texts: Getting polyline length wise in dataframe in layout view in ArcMap via rotation?

texts:
I'd like to rotate my dataframe until a polyline fits dataframe-width wise in the layout view of a .mxd.  If the lines runs north/south I would like the dataframe rotated until the line is shown east/west (the north arrow adjusts accordingly).

I can rotate the dataframe with arcpy.  But how would I check that my polyline is running dataframe-width wise.

Or would I get the envelope of the feature and then do some math to figure out how much the rotation angle should be to have the data frame rotated to meet my goal?

The end result would be to capture the entire polyline in the dataframe in layout view of ArcMap (dataframe-width wise) at the smallest extent possible.  Like an alignment sheet.


texts: More accurate way to calculate area of rasters

texts: In my daily work, I am constantly asked to calculate areas of global raster datasets in geographic projection at 30 arc second resolution. These datasets are normally the result of a Combine operation (a typical example is a vegetation classes combined with country layer). To do this, our unit created a raster dataset with the area of each pixel in geographic projection at 30 arc second. With this area grid, a zonalstat is performed to sum the areas for each class. Since I am not sure how this area grid was created, I always wondered whether this approach is more accurate of just reprojecting the raster in a equal-area projection (from simple tests the results of the two methods are similar). Does anybody experienced a similar situation?


texts: How to create correct vector grid in meter

texts:
I want to create a vector grid, with squares of 100x100m but have some serious trouble in solving this apparently easy task. I already read all the posts, but didn't find a good solution. Here is what I tried:

Using the research > vector grid tool: I always fail to have the right X and Y Values, it is never 100m, even if I have "m" as my default length.I need 100m but it says always something around 59,... if i measure it with the measurement tool. and also the vertical length seems to be slightly bigger than the horizontal length.
Using processing tools: create grid - I have the same problem with length.
Using the MMQGIS create Grid tool: Same problem as 1 and 2.

how can I change the way i calculates distances? I think I have to look somewhere else, but dont know where.

FYI: Running system is Win 8, 64 bit, QGIS Version 2.4,


texts: How to find point on polyline using a starting point and a distance along the curve using ArcObjects

texts:
I want to query a point on a polyline, using a starting point on the polyline, and a distance from the starting point, using Arcobjects.

My inputs
Point A, Polyline A, Distance from pointA.
Output
Point B.

What I tried

var distanceAlong=0.5;
var distanceFrom=0.0;
var pointB=new PointClass();
polyline.QueryPointAndDistance(esriSegmentExtension.NoExtention,PointA,true,pointB,distanceAlong,distancefrom,false)

This didnt work for me, PointB is turning out to be the same as PointA.

QueryPoint is not what I want either since I have a starting point which I need to query from.

How can I find PointB from PointA using a distance along the curve/line ?


texts:
Through the OpenRouteService, I've generated 2 driving time maps from 2 points approximately 60km from each other. Screenshots below - each colour represents 10 minutes.

I'm trying to find the 'half way point' between Alnwick and Duns. Would anyone have an idea how to make this a cumulative travel time map?


texts: Combining two shapefiles which contain overlapping and duplicate geometry

texts:
Well this could be an interesting one. I've got two shapefiles of polylines (contours). I need them to be one nice and clean shapefile covering the same area.

Only issue is, the shapefiles overlap and cover different areas. As below:

Both together become the shape we're after, but of course there's the bit in the middle which gets very pesky.

The tools at my disposal are QGIS, Manifold 8, and SQL Server (and PostgreSQL). I've combined the layers using QGIS and then used MMQGIS to remove duplicate geometry, but this doesn't quite work either as one line which is present in both shapes, for instance, will continue northwards of shape 1, and southwards of shape 2, meaning it is not the same in either. I would need a way to recognise which lines touch (as no line should touch each other in the finished article) and combine their geometry.

I did also start dissolving the merged shapefile QGIS, but that was taking a very long time. Don't know if that's the way forward - maybe I should try it on a smaller sample. I don't know whether that exists in Manifold, which seems a lot faster at things but restrictive.


texts: Get_feature and Get_feature_by_id problem

texts:
Why does the function works with get_feature_by_id and doesn't with get_feature?
closest_point (
    geometry (get_feature (@layer, "num_plgnal" ,  '1')),
    geometry (get_feature (@layer,  "num_plgnal" , '3'))
)

Result with get_feature, the point is not drawn:

closest_point (
    geometry (get_feature_by_id (@layer, '1')),
    geometry (get_feature_by_id (@layer, '3'))
)

Result with get_feature_by_id:


texts: python-netcdf4: Writing file leaves CRS as 'undefined' and extent not correct

texts:
I am trying to write a netcdf file in my file system and then open it in QGIS for visualization purposes. I have 3 numpy arrays with size 234 rows and 217 columns. One contains the latitudes in each grid cell, the second contains the longitudes, and the third one the values of my represented variable. I define a user variable (type 'i4' in netCDF) called "crs", which contains the EPSG code of interest (EPSG:4326), and the proj4 parameter. The file is successfully written in my disk, but when I open it in QGIS it says that the "CRS is undefined" and then I see the map as blank.

This is how I created the metadata for the netcdf file. Note that the variable "proj_wgs" contains the following string: "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

ds = Dataset(basedir_ou.format(ufile_out), 'w', format='NETCDF4_CLASSIC')

dim_lat = ds.createDimension("x", size=234)
dim_lon = ds.createDimension("y", size=217)
dim_tim = ds.createDimension("time", size=1)
dim_ban = ds.createDimension("bnds", size=2)

ds.createVariable('crs', 'i4')
ds.variables["crs"].grid_mapping_name = 'crs'
ds.variables["crs"].epsg_code = 'EPSG:4326'
ds.variables["crs"].proj4_params = proj_wgs
ds.variables["crs"].spatial_ref = proj_wgs
ds.variables['crs'].valid_min = -200
ds.variables['crs'].valid_max = 200

ds.createVariable("x", np.float32, ("x",))
ds.variables["x"].grid_mapping = 'crs'
ds.variables["x"].grid_mapping_name = 'latitude_longitude'
ds.variables['x'].valid_min = -200
ds.variables['x'].valid_max = 200

ds.createVariable("y", np.float32, ("y",))
ds.variables["y"].grid_mapping = 'crs'
ds.variables["y"].grid_mapping_name = 'latitude_longitude'
ds.variables['y'].valid_min = -200
ds.variables['y'].valid_max = 200

ds.createVariable("lon", np.float32, ("y", "x"))
ds.variables["lon"].grid_mapping = 'crs'
ds.variables["lon"].grid_mapping_name = 'latitude_longitude'
ds.variables["lon"].units = "degree_east"
ds.variables['lon'].valid_min = -200
ds.variables['lon'].valid_max = 200
ds.variables["lon"][:,:] = proj_lons

ds.createVariable("lat", np.float32, ("y", "x"))
ds.variables["lat"].grid_mapping = 'crs'
ds.variables["lat"].grid_mapping_name = 'latitude_longitude'
ds.variables["lat"].units = 'degree_north'
ds.variables['lat'].valid_min = -200
ds.variables['lat'].valid_max = 200
ds.variables["lat"][:,:] = proj_lats

ds.createVariable("vble", np.float32, ("y", "x"))
ds.variables["vble"].grid_mapping = 'crs'
ds.variables["vble"].grid_mapping_name = 'latitude_longitude'
ds.variables['vble'].valid_min = -100
ds.variables['vble'].valid_max = 100
ds.variables["vble"][:,:] = array_variable

# These variables are commented for simplicity
# ds.createVariable("time", np.float64, ("time",))
# ds.createVariable("date", np.int32, ("time",))
# ds.createVariable("hms", np.int32, ("time"))
# ds.createVariable("time_bnds", np.float64, ("time", "bnds"))
# ds.createVariable("height", np.float32)

By doing "gdalinfo destination_file.nc", I see the following data:

Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: none associated
Size is 512, 512
Coordinate System is `'
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:".../output_folder/destination_file.nc":lon
  SUBDATASET_1_DESC=[217x234] lon (32-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:".../output_folder/destination_file.nc":lat
  SUBDATASET_2_DESC=[217x234] lat (32-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"/output_folder/destination_file.nc":vble
  SUBDATASET_3_DESC=[217x234] vble (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Thus, it seems that the CRS is not recognized. Then, the size is incorrect, despite I specified the number of rows and columns (234x217), and when I open the metadata in QGIS, it says that the origin of coordinates is "9.96921e+36,9.96921e+36", which is incorrect as well.

Is there a way of correcting for these issues in python-netcdf4?

How can I make QGIS or GDAL to recognize the CRS?


texts: How to export large attribute table from Arcmap and import it into R?

texts:
I intersected 5 shapefiles (you can find the question and data here) to get polygon shapefile. In which, each polygon has ALL the attributes from the 5 layers.

You can find the shapefile resulting from this intersection in this link.

I want to process the attribute table in R. I exported the attribute table as .dbf file. You can find the .dbf file in this link

I tried to import it to R as below

library(foreign)
df <- read.dbf("5layers1.dbf")

I got errors (below) and RStudio crashed so I couldn't get sessionInfo().

As @Vince pointed out below, it is unlikely I will be able to import it to R because of the large size.

Do you have any suggestions how to export large attribute table and import it into R?


texts:
This is a newbie question that is about converting between wgs84 and nad83. I have aerial photos that I am trying to compute lat long coordinates for as closely as possible. There are several factors that make the usual types of corrections tricky, e.g. questionable parameters, ... so there is some trial and error involved.

Thus I need a way to check my progress as I make corrections. At first just checking against google earth was fine, but at this point the 2.5 foot or so difference between my data (nad83) and google earth(wgs84) has become a limiting factor. So I need to transform nad83 to wgs84.

I have been using rgdal in R and dotspatial in C# to do projections, but I think changing between nad83 and wgs84 isn't considered a "projection" (am I right?). Attempting to call project with "+proj=NAD83 ellps=WGS84" just throws an error.

I am fine with writing a 7 parameter coordinate change myself, except I am not sure how to convert either one to (X,Y,Z) coordinates (which I assume means converting from lat long to the surface of the ellipsoids they are based on). The existing questions on nad83 and wgs84 seem to be about which method to use, rather than answering a very basic question like mine, i.e. how does one do it at all?


texts:
Suppose i want to group 2 or more style of different layers into one SLD file how can i do it ?

A pointer towards the working example or document will do.

Or
any alternative styling example such as css will do it it fits in above requirements.


texts: How can I digitize features on an Y,Z plane in ArcGIS?

texts:
How can I geo-reference and digitize features on a Y,Z plane of a vertical cross-section in ArcGIS?
The surface coordinate system is NAD 83 UTM 15.

This is a simplified version of the question found here:
How to digitize directional drilholes from scanned vertical sections?


texts: Instantiating spatial polygon without using a shapefile in R

texts:
So, the usual way we read a shapefile in R is via the maptools package, like this:

sfdata <- readShapeSpatial("/path/to/my/shapefile.shp", proj4string=CRS("+proj=longlat"))

However, I have a use case whereby I don't have a shapefile.shp but instead I have a series of polygon coordinates

16.484375 59.736328125,17.4951171875 55.1220703125,24.74609375 55.0341796875,22.5927734375 61.142578125,16.484375 59.736328125

and its corresponding projection

coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

How do I "instantiate" sfdata (which will be a "polygon object") directly from this data? (without going in a roundabout way of creating a shapefile with these data and then reading from the newly created shapefile)


texts: Changing TextString field in ArcMap annotations not resulting in changed annotations on map

texts:
I am having an issue with ArcMap while trying to save myself a lot of work.

I recently had to individually place hundreds of annotations for a large map. Without getting into too much detail the labels I placed were at the section level which means something like 1st Street would be annotated as 1ST-01,1ST-02,1ST-03, etc for different sections of the street. However, the client also wants a map at the branch level instead of the section level which means the annotations would just say 1ST without any section. The "-" is pretty useful for this and I opened up the annotation attribute table in excel,  used a formula to remove everything after "-" and bam I got branch labels without having to place them again from scratch.

I then joined this excel sheet to a copy of the section annotation layer and used the field calculator to make the TextString field equal to the trimmed down Section IDs. This worked and all the TextString fields are updated, but on the map they have not changed.

Is it because I did this in excel and tried to join it back?

Is there an easier way to remove everything after a "-" in the field calculator?


texts: Calculating Autocorrelation in Different Zones

texts:
I'm an intermediate-ish user of ArcGIS 10.2 (win7 pro, hardware robust enough to handle very large datasets and computationally intensive analysis) with perhaps a simple problem that's evading me right now.

I have a raster dataset with different classes of land values. I am interested in whether a dichotomous variable corresponding to presence/non-presence of a land classification is spatially clustered within different zones of a polygon feature. Basically I'd like to do something similar to tabulate area or zonal statistics, whereby I could obtain Moran's I for each individual zone in the feature class. To start I'd convert the raster to a polygon, but based on what I'm seeing I would have to compute Moran's I over the entire overlapping polygon with the 'zones,' rather than obtaining values for each individual one. (I could do it one by one, but this seems cumbersome.)


texts:
I am using GRASS GIS 7.4 and I have a DEM as points in a vector format and want to turn it into a raster.

I tried using v.to.rast and v.to.rast3, but the elevation is not one of the attributes so am struggling to set the source of raster values. I get the following options

and selected z, assuming this is the elevation. I run the function and it works, it creates a raster, but the raster is empty.

I did a query in one of the vector points and could see that the elevation is in Point_height, but I do not know how to use this as the source of the values for the raster. Does anyone have any solution for this?