texts: Using QgsFeatureRendererV2 methods writeSld() and loadSld() in QGIS

texts:
I have a vector layer in QGIS and using Python I would like to be able to export the layer styling to a string variable and import styling in a string variable to the layer again (pseudo code):

style_string = layer_1.rendererV2().export_style()
layer_2.rendererV2().import_style(style_string))

I have found the two methods writeSld(...) and loadSld(...) in the QgsFeatureRendererV2 but I cannot find any examples of how to use them together.

The only clue I have found was this blog post which describes how writeSld(...) can be used. But when I try to import the XML document back into the renderer using loadSld(...) QGIS gives me a minidump.

Question:
Can anyone give me an example of how this can be done?


texts:
I have a shapefile that contains all of the property parcel lines in King County, WA. We are currently displaying this data through mapbox on our site with no issues. We would like to be able to show (estimated) parcel dimensions (feet).

Is there anyway to extrapolate that information in qGis?

My thinking was to possibly add it to the attributes table so it could then be displayed as a label?

We'd need it to display something like this:

----------60 ft----------


texts: Should I use a discrete or continuous scale for coloring a chloropleth?

texts:
I am making a choropleth for an upcoming local race using D3. There are two candidates running for an election. For any given polygon on my map, if the red candidate has more votes, the color should be red. If the blue candidate has more votes, the color should be blue. If the blue candidate is winning by a lot, the color should be a deeper blue. If the red candidate is winning by a lot, the color should be a deeper red.

Should I use a discrete or continuous scale for this? Meaning, should I create a color ramp from red to blue for a given saturation or intensity. Or should I create bins such that if red/blue falls in a certain range, the polygon is assigned one of a few colors? From tutorials, it seems like most people make bins.


texts: Finding the optimal boundary polyline in QGIS

texts:
I am looking to find an optimal boundary line between two geographically separated sets of points (a vector layer with S and N points representing two genetic lineages within a species). I would like to do this by minimizing the difference between the sum of inverse square distances to S points and the sum of inverse square distances to N points at a given latitude. I would then connect these points across latitudes. The polyline would end where the sum of inverse sequence distances fell below some threshold value. I may also like to vary symbology along this boundary based on the summed distances. The calculations and output could be a Raster but I would prefer to convert this to a smoothed polyline.

I can do these calculations manually but it is laborious and not-transferable to other projects. Can anyone help with ideas on how to do this within QGIS?


texts: QGIS 'join attributes by location' algorithm produces strange results for min, max

texts:
In QGIS 2.14.11 LTR I have a point shapefile with (besides others) an attribute kd which is either 1 or 0 (never NULL) and a polygon dataset. I try to calculate minimum and maximum values for kd for each polygon containing points and I expect to be minkd and maxkd fields after to be either 1 or 0 also, but never NULL or any other values greater than 1.

My join attributes by location algorithm is set up as follows:

As can be seen below (selected resulting polygon in the middle), values for minkd and maxkd are 6.000 in this case, but NULL and a variety of other values differing from 1 and 0 occur also.

I fairly can't understand how this behaviour can be explained, has anyone observed similar behaviour? Any hints to overcome this?


texts:
I have a piece of code I'm using to replace layer sources and then zooms to the new features. The part where I replace the workplace source works, but I can't seem to get it to zoom to the extent of the new features. I'm running ArcMap 10.2.

import arcpy
mxd = arcpy.mapping.MapDocument('CURRENT')
mxd.findAndReplaceWorkspacePaths(r"\\Client\L$\17TSLB002_Morice\A95440_1_FULT_21\Submission_GDB\A95440_1_TBA_Digital_Sub_Dec11_2017.gdb",
                                r"\\Client\L$\17TSLB002_Morice\A95441_1_NO_6\Submission_GDB\A95441_1_TBA_Digital_Sub_Dec1_2017.gdb")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")
lyr = arcpy.mapping.ListLayers(mxd, "BLOCK", df)
df.extent = lyr.getSelectedExtent()
df.scale = 5000
arcpy.RefreshActiveView()
mxd.saveACopy(r"\\Client\L$\17TSLB002_Morice\A95441_1_NO_6\Submission_PDFs\A95441_1_CHR_11by17_Portrait.mxd")
del mxd

I get this error though when I run the script.

Runtime error
Traceback (most recent call last):
  File "<string>", line 6, in <module>
  File "e:\sw_nt\arcgis\desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
  File "e:\sw_nt\arcgis\desktop10.2\arcpy\arcpy\mapping.py", line 1500, in ListLayers
    result = mixins.MapDocumentMixin(map_document_or_layer).listLayers(wildcard, data_frame)
  File "e:\sw_nt\arcgis\desktop10.2\arcpy\arcpy\arcobjects\mixins.py", line 820, in listLayers
    raise TypeError(str(type(dataframe)))
TypeError: <type 'list'>


texts:
Are there any free, vector, mercator projection maps available online? I'm not wanting satellite imagery, but more a stylized map along the lines of this:

However, I need it to be vector-based such that I can create sharp tiles for greater levels of zoom.

Ideally, I would like an Illustrator file for something resembling the above map.


texts: How to Convert from GMT Line files to Shapefile for QGIS?

texts:
I have several line text-files of (used in GMT) saved in the format (this is a small sample):

-111.29 42.28
-111.29 42.31
-111.29 42.33
>
-111.37 41.99
-111.38 42.02

I am trying to convert this to a Shapefile to use in QGIS. I know people convert from shp to GMT with ogr2ogr can you do this in reverse GMT to Shp? I tried using ogr2ogr and had no luck. Maybe I did something wrong or it is not possible with that tool.

I have tried to import this format in QGIS via the delimited text layer tool, but can only figure out how to get points when I need lines.


texts: Converting spherical coordinates to geographic coordinates in complicated netcdf file

texts:
I am very new with netcdf file. And currently have to extract the data from such netcdf which looks like below:

8 variables (excluding dimension variables):
        double time_bnds[bnds,time]
        int rotated_latitude_longitude[]
                grid_mapping_name: rotated_latitude_longitude
                grid_north_pole_latitude: 39.25
                grid_north_pole_longitude: -162
                north_pole_grid_longitude: 0
            float lat[rlon,rlat]
                standard_name: latitude
                long_name: latitude coordinate
                units: degrees_north
                bounds: lat_vertices
            float lon[rlon,rlat]
                standard_name: longitude
                long_name: longitude coordinate
                units: degrees_east
                bounds: lon_vertices
            float lat_vertices[vertices,rlon,rlat]
                units: degrees_north
            float lon_vertices[vertices,rlon,rlat]
                units: degrees_east
            double height[]
                units: m
                axis: Z
                positive: up
                long_name: height
                standard_name: height
            float tasmax[rlon,rlat,time]
                standard_name: air_temperature
                long_name: Daily Maximum Near-Surface Air Temperature
                comment: daily-maximum near-surface (usually, 2 meter) air temperature.
                units: K
                cell_methods: time: maximum
                history: 2015-05-06T15:11:26Z altered by CMOR: Treated scalar dimension: 'height'. 2015-05-06T15:11:26Z altered by CMOR: replaced missing value flag (-1e+20) with standard missing value (1e+20).
                coordinates: height lat lon
                missing_value: 1.00000002004088e+20
                _FillValue: 1.00000002004088e+20
                associated_files: gridspecFile: gridspec_atmos_fx_MPI-CSC-REMO2009_rcp45_r0i0p0.nc
                grid_mapping: rotated_latitude_longitude

5 dimensions:
        time  Size:1825   *** is unlimited ***
            bounds: time_bnds
            units: days since 1949-12-01T00:00:00Z
            calendar: proleptic_gregorian
            axis: T
            long_name: time
            standard_name: time
        rlat  Size:103
            units: degrees
            axis: Y
            long_name: latitude in rotated pole grid
            standard_name: grid_latitude
        rlon  Size:106
            units: degrees
            axis: X
            long_name: longitude in rotated pole grid
            standard_name: grid_longitude
        bnds  Size:2
        vertices  Size:4

Although, there are many explanation of different netcdf files but I found my situation is a bit complicated in terms of number of variables and dimensions.
Here I am confused with float lat[rlon,rlat]; why rlon and rlat is here with normal latitude? Similarly with the case of float lon[rlon,rlat]. And why rlon and rlat separately present in dimensions also?

Eventually I need to extract the temperature with geographical coordinates (latitude North and longitude East).


texts: Performance of Shapefile and PostGIS with GeoServer?

texts:
I'm going to build a web application, which will render around ten maps using GeoServer. I will be also changing the attribute data in the application for some layers.

What is the recommended choice for performance and deployment: A map rendered from GeoServer using Shapefile or PostGIS database ? ( or they don't really matter)

I don't deal with any spatial query.


texts: How to calculate area in square kilometers of a OSM boundary

texts:
I have tried using measurement plugin in JSON.

It does not work if the boundary is made of multiple selected lines. That plugin works when it is a single line for the whole boundary.

Is there any way to do this?


texts: Create 3D-SHP / MIF from 2D-SHP and height attributes

texts:
I have two .SHP-files:

Contains building footprints as polygons with IDs in the attributes
Contains building height attributes (+ other information not relevant for me) as points with the same IDs in the attributes

I want to combine the building footprints with the height information to create a 3D-Shape or MIF/MID-file.

Is there a simple GUI-way for noobs like me or do I have to do some fancy Python stuff?


texts: Automated quality control and assurance - geodata?

texts:
I've recently started looking at the issue of automated quality control and assurance of geographic data.

What software exists that can be used for accomplishing this? Esri has got their product 'Data Reviewer' which seems interesting, are there any comparable alternatives to this product?

I would like to tag this question with quality-control and quality-assurance, but can't do so.


texts:
I'm working with ArcGIS v10.1 and MS Access 2010.  Can you export a datalayer directly into Access 2010 (.accdb)? This was possible under using older versions of ArcGIS with Access 2003.

I created a DataBase Connection using ArcCatalog to reference the Access 2010 database and I can successfully import a table from Access into ArcGIS. However, when I try to export and point back to the same Access 2010 database using the same DataBase Connection, the .ODC I created is not available.

I can work around this by exporting the data into a spreadsheet and then importing the spreadsheet into Access.


texts:
Access to XMLHttpRequest at 'http://localhost:8080/geoserver/decuong_doan/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=decuong_doan%3Aho_chi_minh_quan&STYLES=&LAYERS=decuong_doan%3Aho_chi_minh_quan&INFO_FORMAT=application%2Fjson&FEATURE_COUNT=50&X=50&Y=50&SRS=EPSG%3A3405&WIDTH=101&HEIGHT=101&BBOX=106.5312897896764%2C10.584902670383462%2C106.72903557538959%2C10.782648456096654' from origin 'null' has been blocked by CORS policy: Response to preflight request doesn't pass access control check: No 'Access-Control-Allow-Origin' header is present on the requested resource.

How can I resolve this problem? I write code in order to make a Web GIS-based on GeoServer and OpenLayers. But when I get features of layer I can't do it.
This below is this code get features:

map.on('singleclick', function (evt) {
document.getElementById('info').innerHTML = "Nhìn thấy dòng này thì làm lại đi :) !!!";
var view = map.getView();
var viewResolution = view.getResolution();
var source = hcm.getSource();
var url = source.getFeatureInfoUrl(
evt.coordinate, viewResolution, view.getProjection(),
{ 'INFO_FORMAT': 'application/json', 'FEATURE_COUNT': 50 });
if (url) {
$.ajax({
type: "POST",
url: url,
contentType: "application/json; charset=utf-8",
dataType: 'json',
success: function (n) {
var content = "
";
for (var i = 0; i < n.features.length; i++) {
var feature = n.features[i];
var featureAttr = feature.properties;
content += "
Tên quận:" + featureAttr["name_2"]
+ "	Số hiệu quận:" + featureAttr["id_2"]
+"
"
}
content += "
";
$("#info").html(content);
}
});
}
});


texts: Extracting features with Export Training Data for Deep Learning tool in ArcGIS Pro

texts:
The stats file generated from the tool:

images = 3400 *3*30*30

features = 6412

features per image = [min = 1, mean = 1,89, max = 10]

classes = 1

cls name cls value images features min size mean size max size

Building 1 3400 6412 0,00 0,00 0,00

It is clear that I have 6412 features, but what exactly are these features and what does it mean for there to be a min 1 feature in the image or max 10 features in an image?


texts:
I am having problems using the raster calculator. Some pixels of my GeoTIFF file have the value 255 on bands 1, 2 and 3. I want to replace this value with 0 on the three bands.

I tried with this formula but it does not look good:
("test@1" != 255)+("test@2" != 255)+("test@3" != 255) * "test"


texts:
I need to do the following in QGIS and hope so much anyone could help.

I have a vector layer with points. Each point represents a firm with the number of employees in the attribute table. I want to buffer around every point, so that the buffer's attribute table contains the sum of employees within the buffer.

It would be so great if the sum to each point could be written to the attribute table.

Thanks for your help.


texts: In ArcMap 10.3, I am making a map, which contains a road layer represented by dual lines. Unfortunately, at some intersection points, the lines are interrupted and tend to create wedge-like gaps (see red circles in figures below). I am struggling with it for quite a long time, but I still have not found any viable solution. One idea was to slightly change the position of the vertices at the intersection points. Indeed, it worked, but this approach is not very well suitable for the relatively large dataset (with numerous intersections) I have. Further, and more importantly, I don't want to disrupt the topology. Do you have any suggestions how to resolve this matter?


texts: Achieving Super or Subscript graticule labels in QGIS Composer windows?

texts:
It is quite common for paper topographic maps to label the map graticule with ordinate values in a mixture of superscript and standard digits, eg. ⁴123⁰⁰⁰ to emphasise the significant digits.

The QGIS Map composer allows graticules (termed grids) to be drawn on the map, and the ordinate labels to be customised, but does not recognise HTML tags such as <sup>.

Is there any other way to achieve this effect?


texts: Using Visual Studio 2008 or 2010 for ArcGIS Desktop?

texts:
Is there any differences, pros or cons between Visual Studio 2008 and 2010?

I am going to buy Visual Studio and just want to be sure which one will work the best and have the most support going forward.

I am guessing Visual Studio 2010 but just want to be sure.

Will I need the professional version or will standard work fine?

I am a beginner using ArcGIS Desktop 10.


texts:
I’ve got a shapefile of Switzerland. Shown Polygons are 6-digit ZipCodes. Whenever I dissolve the file for example by 4-digit ZipCodes, it results in incorrect polygons.

Unfortunately some polygons borders aren't touching. If I zoom in on the orignal shapefile, it seems like there are empty spaces between the polygons.

I'm using QGIS and R and I've tried a lot.
This is my first time using GIS data.


texts: Add multiple fields to attribute table in ModelBuilder

texts:
I want to add a number of columns to an attribute table in a few steps using Model Builder.

I think this link gives the answer https://community.esri.com/thread/182530-how-do-i-add-multiple-fields-to-a-table-in-model-builder-arcgis-104
but I am having trouble replicating it, because I don't understand the basics.

I created a table that has the attributes I want, but I now need to figure out how this fits in the model builder.

I have the layer I am interested in Property.shp,

The input to Iterate is MyTable, which gives the 2 outputs I_MyTable and Value.

I don't understand how this links in to Add Field - or if that is even the correct next step.


texts:
How open layer's attribute table as tab?

Before reinstall OS on QGIS 3.8 I opened layer's attribute table as tab, after on QGIS 3.10 I can't find this option.


texts: Semi-Automatic Classification Plugin - Crush dumped

texts:
I installed Semi-Automatic Classification Plugin but it doesn't work.
After starting QGIS, program displays a message "Please restart QGIS to executing the Semi-Automatic Classification Plugin. Possible missing dependencies."
Closing the program gets an error "Crush dumped", every time.

I can't uninstall or reinstall this plugin from plugins repositories interface, this also causes "Crush dumped".

I use QGIS 2.14.7 and 2.16.3. Windows 7 Pro 64bit.


texts: Inconsistent raster image display in ArcMap

texts:
I loaded a WorldView 2 satellite image into ArcMap 10.5, however the tiles display in an inconsistent color and shade.  All 9 tiles are from the same WorldView image.

How can I get them to display in a consistent pattern in ArcMap?


texts: How do I access my newly added fields in a spatialite layer from QGIS?

texts:
Forgive my language, I am a newbie in QGIS. I created a spatialite layer and added 4 fields. As I was working on the project I realised that I didn't need one of the fields, and I needed to add 3 different fields, to make them 6 in total.I encountered the limitation of SQLite not being able to delete created fields. So I looked for a GUI option of modifying the database. I successfully deleted that one column and added 3 new ones. Now the problem is when I open the layer in QGIS to edit the feature attributes, only the old 4 columns are still available. My spatial layer initially had these attributes:
.
After updating the SQLite database using DB Browser for SQLite, below is the structure of my table from DB Browser for SQLite:

Coming back to QGIS to add more information to newly created fields, I still find the old 4 fields:
.

Anyone to help so that my layer loads all fields for feature attributes please?


texts: GEOTIFF Binary File with Color Display Issues

texts: I have a GEOTIFF in QGIS that shows values of 0 and 1 and I want the 0 values to have no fill so that background vector layers can show through and so I can easily overlay things in a print composer for final production of the map. This is an easy task in ArcMap but doesn't seem to be possible with QGIS raster layers because they require you to have a color ramp/a color for each class when you change styles. How do I get this to work? Having 0 values be white isn't good enough for what I am mapping.


texts: Is it possible create a custom UTM zone?

texts:
I'm trying to detect the CRS (coordinate reference system) of a layer, but it is between

WGS 84 / UTM zone 32N

and

WGS 84 / UTM zone 33N

In fact the definition shows srsName="GEOPOIUTM32.5" srsDimension="2"

But I can't find any SRS like this in QGIS.


texts: How to generate map tiles including world files in QGIS?

texts: Is there any way to export bitmap tiles including world files from a QGIS project? I've tried QTiles but it seems not to support world files.


texts: EE Profiler shows by far biggest compute time on "Reprojecting geometry to SR-ORG:6627"?

texts:
EE Profiler shows by far the biggest compute time on "Reprojecting geometry to SR-ORG:6627"

Is there a way to import a shapefile into EE already in SR-ORG:6627 projection?

I've tried using QGIS to save a shapefile already in SR-ORG:6627 projection (aka Google Mercator), but without success.  The Profiler still shows the same compute time whether the shapefile (FeatureCollection) is in Geographic Decimal Degrees WGS 84 or in SR-ORG:6627 projection (aka Google Mercator)

Has anyone found a way to import a shapefile as a FeatureCollection already in SR-ORG:6627 so this high Compute time shown in the Profiler can be avoided?


texts: Get tile images directly from sqlite database

texts:
I'm trying to get map tile images to display in my application. I created a sqlite database and used WMS Lite server to request map tiles using GetMap syntax.

wmslite?SERVICE=WMS&VERSION&REQUEST=GetMap&BBOX=xmin,xmax,ymin,ymax&CRS=EPSG:3857&STYLES=&FORMAT=image/png&WIDTH=512&HEIGHT=512&LAYERS=test

But WMS Lite is not reliable because it crashes so often. So now I'm trying to get the tiles directly from the sqlite database. Is there anyway that I can download the tiles similar to the GetMap syntax.


texts: Adding text labels to longitude and latitude points in shapefile map by plot()?

texts:
I have read my shapefile, and have also read my longitude-latitude points (in a csv file), that I want to overlay on the map. I have further created a character vector with the location names that I want to appear as labels on point on the map. I have tried to use the text function but the labels do not display.

Below is my code:

# Read shapefile
map1 <- readOGR(dsn=".", layer="mymap")

# Read CSV file
cvsfile <- read.csv('./Hcenters.csv', sep=',', header=TRUE, as.is = TRUE)

# Convert CSV data frame to SpatialPoints class
coordinates(cvsfile) <- ~ lat + lon

# Make sure the two files share the same CRS
cvsfile@proj4string <- map1@proj4string

# Create label name text vector
names <- c("Bugono HC IV", "Namungalwe HC III", "Lubira HC III", "Busesa HV IV", "Busowobi HC III", "Baitambogwe HC III", "Mayuge HC III", "Kityerera HC IV", "Kigandalo HC IV")

# Plot shape file
plot(map1, border="red", xlim=c(33.0, 33.7), lwd=2, axes=TRUE)

# Overlay longitude-latitude point locations to the map
points(cvsfile$lat ~ cvsfile$lon, bg="red", pch=15)

# Add label names to latitude-longitude points
text(cvsfile, labels = names, cex=0.7, pos=3)

However the labels do not display.

Where I am I going wrong?


texts: I have several adjacent (county level) TIGER road files.    I want to merge them and create one contiguous file of snapped roads for the purpose of doing network analysis.   Merging them is easy but how do I get the ends of the roads to snap together properly?


texts: Adding address (newly-built house) to existing street on Google Maps

texts:
How can I add my address (newly-built house) to an existing street on Google maps?

It shows on Bing maps OK but not on Google maps.

Existing answers to this question seem to be outdated.

Surely this is a common problem every time a new house is built.


texts: Make a Loop for Monthly NDVI in Google Earth Engine?

texts:
I've computed maximum NDVI for each month. but I don't know how can I write this code through a loop in Google Earth Engine?

code link: https://code.earthengine.google.com/99fd025cdb0688dce74ea5fa0966a2c2

Map.centerObject(table);
Map.addLayer(table);

// ndvi jan

var jan = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-01-01','2018-02-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('jan');

// ndvi feb

var feb = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-02-01','2018-03-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('feb');

// ndvi mar

var mar = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-03-01','2018-04-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('mar');

// ndvi apr

var apr = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-04-01','2018-05-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('apr');

// ndvi may

var may = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-05-01','2018-06-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('may');

// ndvi jun

var jun = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-06-01','2018-07-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('jun');

// ndvi jul

var jul = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-07-01','2018-08-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('jul');

// ndvi aug

var aug = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-08-01','2018-09-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('aug');

// ndvi sep

var sep = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-09-01','2018-10-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('sep');

// ndvi oct

var oct = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-10-01','2018-11-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('oct');

// ndvi nov

var nov = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-11-01','2018-12-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('nov');

// ndvi dec

var dec = ee.ImageCollection("COPERNICUS/S3/OLCI")
.filterBounds(table)
.filterDate('2018-12-01','2019-01-01')
.map(function(img){
  return img.normalizedDifference(['Oa17_radiance','Oa08_radiance']);
})
.max()
.clip(table)
.rename('dec');

var ndvistack = jan.addBands(feb).addBands(mar).addBands(apr).addBands(may)
.addBands(jun).addBands(jul).addBands(aug).addBands(sep).addBands(oct)
.addBands(nov).addBands(dec);


texts: Changing default precision and scale of number fields in ArcSDE geodatabase table?

texts:
After I change the precision and scale for my short, long, float, and double field types when creating a table in ArcMap SDE geodatabase, ArcMap will reset my precision and scale to default values. For instance, I wanted my short integer to have a precision of 3 in one table and 2 in the other so I could test some python code. ArcMap changed them both to precision of 5 automatically and the python tool isn't given a proper test.

Please note that this is a different but related issue that was addressed here: Specifying precision and scale of new field using ArcPy? I know how to specify the precision, but ArcMap isn't keeping my settings.

If this makes sense stop reading here. If you want more detail, please continue reading:

I am creating an enterprise geodatabase in ArcMap 10.2.1 for the purpose of testing two tables: a source table and a destination table. I am building a tool that will describe the tables, put the field names as keys in a dictionary and the field objects as values, populate a list of all matching fields, then a series of if statements remove the fields that may have the same name, but incompatible data types, precisions, scales, and field lengths. Incompatible is different than non-matching. I created if statements that allow a short integer in a source table to be compatible with a long integer(or a float, or a double) in a destination table. But this cannot be reversed. Only a short can go into a short. But if they have different precisions, a smaller or equal precision can go into a larger precision, but a larger precision cannot go into a smaller one.

This is why I need ArcMap to keep my varying precision values across different fields. I need to be able to test every possibility.

#           Pass SmallInteger types if going into a SmallInteger.

            if source_type == 'SmallInteger' and destination_type == 'SmallInteger':
                pass

#           Filter out incompatible precision.
                if source_precision <= destination_precision:
                    pass
                if source_precision > destination_precision:
                    matched_column_list.remove(name_1.lower())
                    continue

That is my issue. If anyone can get my enterprise geodatabase to behave, please let me know!

I'm using ArcGIS 10.2.1.


texts: How to update ArcGIS Server MapServer service without re-cache?

texts: I published a cached map service with set of attributes. I need to add / remove published attributes, but I want to use old cache, because there was no change in the geometry. Is there a way how to use old cache? I tried to create new service and copy cache into its cache directory, but ArcGIS seems to recache anyway.


texts: Parameterized insert query to PostGIS fails

texts:
I have a PostGIS database table that I need to insert data into. I'm using Python 3 and the pg8000 interface to Postgres. Somehow, parameterized inserts do not work. If I do

iq="INSERT INTO points(location) VALUES(ST_GeographyFromText('SRID=4326;POINT(:lat :lon)'))"
cursor.execute(iq, {"lat":lat,"lon":lon})

or

iq="INSERT INTO points(location) VALUES(ST_GeographyFromText('SRID=4326;POINT(? ?)'))"
cursor.execute(iq, (lat, lon))

Then the query fails with a ProgrammingError:

pg8000.errors.ProgrammingError: (b'ERROR', b'XX000', b'parse error - invalid geometry')

This, however, works:

iq="INSERT INTO points(location) VALUES(ST_GeographyFromText('SRID=4326;POINT({0} {1})'))".format(lat,lon)
cursor.execute(iq)

What am I doing wrong in the parameterized query?

Edit:
Forgot to mention that, the database query works, if done with PgAdmin against the PostGIS base with some random point coordinates.


texts: Polygonize tool not giving output in QGIS?

texts:
I am very new to QGIS (long time ESRI user that is making the transition).

I am using QGIS 3.0.3 and I am trying to make polygons from a raster layer by value.  Unfortunately every time that I attempt to do this I get the following error in the log.  I am certain (given my time learning new programs) that this is a silly human error but I also know that those can eat up an entire day.


texts: Visualizing OpenStreetMap tiling scheme in QGIS

texts:
OpenStreetmap uses a tiling scheme which is adapted to usual access patterns, as shown in the figure below. I would like to visualize it, as shown in the figure, in QGIS.
Any suggestions on how can I do it?


texts:
I have yield data of a field. The final data are represented by quarter million of polygons.

The width of the harvester was 30 ft and it measured and stored the yield data after each 1.2 ft. So, the raw data contained locations of some points (where the yield was measured) and the yield. For practical use, the points were converted to polygons with the size of 30ft by 1.2ft. I have been asked to concentrate (dissolve) every 5 polygons and produce one, thus the resulted polygon will be 30 ft by 6ft (1.2*5). The new polygon will also contain the average yield of the old 5 polygons. The polygons are in a sequential order from left to right and I don't have to do change the polygons length which is 30 ft (harvester width).

The polygons are shown in the image. However, in this image, the polygons have some blank space between them, please assume there will be no blank space.


texts: PostGIS - default spatial reference and Z value on insert or geometry update?

texts:
I have a spatial table with a geometry column geometry(MultiPolygonZ,4326).

We are trying to handle inserts and updates of geometries as GeoJSON from other systems of varying vintage. Everything is in reality 4326, but sometimes a geometry arrives that lacks the "crs" member, or the geometry might be XY instead of XYZ.

So those cases fail with either

ERROR: Geometry SRID (0) does not match column SRID (4326)
ERROR: Column has Z dimension but geometry does not.

I'd like to know if I can somehow set defaults on the geometry column for the SR (4326) and Z (0.0) to handle those cases.


texts:
I want to join shapefiles (or a shapefile and a data file) based on their mutual FIDs and can't see the FID in the "Target Field" dropdown on the join form.
I believe this is a basic need to have joins based on FIDs, so either I am doing it the wrong way or maybe something is missing here. {QGIS 2.8.1 on Linux}


texts: I'm trying to import united_states.osm into Spatialite and spatialite_osm_net just quits with no error

texts:
When I try to run spatialite_osm_net using united_states.osm as input, the program simply stops after Windows presents a pop-up saying the program ended with an error.

No other information is presented.  The failure occurs in the first phase of spatialite_os_net processing after about 3 hours of processing.

I merged the downloads from cloudmade and unzipped the combined .bz2 file.

I created the output DB with spatialite_gui.

The program level is 3.0.0.

There is more than 250GB of space on the disk and wen te failure occurs the output db is only 7GB.

Has anyone ever successfully done this?

Regards,
Jim


texts:
I a little bit stumped.
I have the entire world as a raster with values ranging from 0-7.
I am wanting to extract the highest values (ideally only those from 6-7) as points.
Is there any simple way i can do this?
Current i am manually putting points in the high areas that i can visually see and then using "Extract values from raster" to give the points a value.
This is tedious and also isn't 100% accurate and i also have to produce it on many rasters...
Here is an image so you can see what i mean (Just put in NZ for an example) Basically just wanting points in all the red areas (values 7-8) and then i can Extract values from Raster.


texts: CalculateField_management: Different results than in the field calculator

texts:
I`m using ArcMap 10.2.2 and Python 2.7.5

I´m tryig to calculate values with the CalculateField_management tool.

my script:

def mthree(mthre, mtwo):
    if mthre >= 1:
        return -3
    if mtwo >= 2:
        return -5
    else:
        return None

try:
    arcpy.CalculateField_management(FC, "EPK_Gesamt", mthree("countmd","countmz"),"PYTHON_9.3")
    print "done"

The Problem:
It writes in every row the same value (-3)
However if I do the same calculation within the Field Calculator it works great!
Has anyone an Idea what I´m doing wrong?


texts:
I am using QGIS.
I have a .csv file containing 3 pieces of information on each row; name, longitude and latitude.
I would like to turn each row of the .csv file into a vector point.

The file will load into QGIS when I open it as a vector layer, but whatever I do, the screen remains blank.

Here is the metadata for the file:

General:
Storage type of this layer: CSV
Source for this layer: G:/VECTOR FILES for qgis/Maritime Traffic Visualisation/ships2.csv
Geometry type of the features in this layer: No geometry
The number of features in this layer: 20189
Editing capabilities of this layer: Add Features
Extents:
In layer spatial reference system units : unknown extent
Layer Spatial Reference System:
+proj=longlat +datum=WGS84 +no_defs

I'm fairly new to GIS and this is probably something very simple, but your help would be appreciated.

Here are the first 10 rows:

name    lon lat
NAPOLEON I  -79.5231    8.66208
ZZ-9 VERTROUWEN 3.91968 51.64689
MIN LONG YU 66628   118.4183    24.12905
ZZ-5 CONSTANT PETRUS    3.92033 51.64691
MIN JIN YU 8401 118.574 24.33855
MIN LONG YU 60209   118.4849    24.15916
JINTANGYU03817  117.7681    39.07499
GOLDEN SCEPTRE PD50 -1.7716 57.5045
MIN  LONG  YU  62503    118.4402    24.27062
THILANT 4.43238 51.88438
VIKING  4.431762    51.23607
TYCHA   4.12002 51.9541
ILHAN AGA   27.981  43.199
ANKIE   -0.7144283  53.60629


texts: Updating the legend of a joined SDE feature layer in Python

texts:
Is it possible in Python to apply the legend from a template layer file to a SDE based feature layer with a GDB table joined to it where the value field being used for the unique values symbology is in the joined table? The script I have worked in ArcGIS 10.2 but since we moved to 10.4.1 the 'arcpy.mapping.UpdateLayers' tool is no longer functioning as it did before.

The problem I see is with the field name after the join has taken place. For example, the template layer and respective GDB table both have the field "GDE_Code", but once the table is joined to the feature class this field becomes something like "Building.GDE_Code". So when I run the 'arcpy.mapping.UpdateLayers' tool in the script it will apply the legend from the template but no features will be displayed.

If I were to export the Building join to its own template layer and then use that one in the script on the joined Building feature layer everything works perfectly; however, is it possible to make the template generic so that it will work on other feature classes?

I've also tried using the 'arcpy.ApplySymbologyFromLayer_management' tool but it doesn't do anything.

This has been attempted as both a stand-alone geoprocessing tool and within the Python window of a Current ArcMap session.

One resolution I have seen is to export the joined feature layer to its own feature class so that the field name is just "GDE_Code", but ideally the goal is to perform this operation on a feature layer that is linking a SDE feature class to a file geodatabase table.

Here is the offending portion of script:

layers = arcpy.mapping.ListLayers(mxd, "*", df)
for layer in layers:
    layer2update = arcpy.mapping.ListLayers(mxd,layer,df)[0]
    layer2update.visible = False
    layerName = layer2update.name
    if layerName in FTRTables:
        GDBtable = os.path.join(FTRepGDB,layerName)
        arcpy.AddJoin_management(layer2update,"facilityIDFK",GDBtable,"facilityIDFK","KEEP_ALL")
        desc = arcpy.Describe(layer2update)
        geometryType = desc.shapeType
        if geometryType == 'Polygon':
            arcpy.mapping.UpdateLayer(df,layer2update,Polygon,True)
        elif geometryType == 'Point':
            arcpy.mapping.UpdateLayer(df,layer2update,Point,True)
        elif geometryType == 'Polyline':
            arcpy.mapping.UpdateLayer(df,layer2update,Line,True)
        else:
            pass
    else:
        pass


texts: AGOL "Elevation Layer (Hosted)" is not displaying outside of a narrow scale range

texts:
28/08/17 Edit: Solved. None of the errors or tools will specify this, but you need to manually generate pyramid layers on the input raster to have it display correctly as an elevation surface. Certainly works on Server, testing in AGOL now.

Edit 2: AGOL Approved.

I am attempting to build a local scene in Pro, host it on AGOL, and embed it it in a Story map. Unfortunately I'm having major problems displaying custom elevation layers.

The end goal is something similar to the scene found here:

https://learn.arcgis.com/en/projects/get-started-with-the-scene-viewer/lessons/explore-a-scene.htm

The issue is that while their elevation layer ("Palm Beach County LADS Surface") displays at all scales, mine only loads at extreme levels of zoom:

http://imgur.com/a/xExA0

i.imgur.com/kvCx02M.gifv

Elevation services can't be produced with the usual tile publishing workflows (i.e. right click, share as web layer, specify your configuration) because they need to be in a very specific format (LERC). If a specific process isn't followed the tile layer publishes, but is treated by AGOL as a flat layer to be draped rather than an elevation surface.

My workflow is that defined in the Pro documentation:

Select a cache tiling scheme (or create one using "Generate Tile Cache Tiling Scheme") - this specifies the tile format, the scales at which tiles should be created, and so forth.
Use "Manage Tile Cache" to build a cache of tiles according to the specifications of the tiling scheme defined above from an input raster.
Export the tile cache to a tile package using "Export Tile Cache".
Upload to AGOL and publish as a hosted layer.
Add the resulting "Elevation Layer (Hosted)" to a scene.

I've experimented with local vs. global scenes, custom tile cache schemes (building caches at specific scales, duplicating the scale breaks & other settings found in both the LADS surface linked above and AGOL basemaps, producing a single scale, etc.), aggregating the raster to 5m cells from 0.5, and publishing smaller features. Nothing seems to work.

One thing to note is that while my features are hosted on AGOL ("Elevation Layer (Hosted)"), the Esri example is being served from a server and lacks the "Hosted" post-fix. While hosting it on a server would certainly fix the issue, it's not something I'm able to do for this project.

What's really fun is that opening a hosted elevation layer in AGOl and hitting the settings tab displays an "Elevation Layer (hosted) Settings" section, which fades out within half a second. If you open a normal tile layer (hosted) in the same way, that tab allows you to re-publish tiles at specific zoom levels.

Does anyone have experience with this? Specifically: do you know if this is a built-in limitation, or something I can overcome? I'm beginning to suspect it's a symptom of the limited functionality of AGOL.


texts: Convert 16 bit raster into 8 bit using rasterio reproject

texts:
I am trying to reproject a raster image from EPSG:4326 to EPSG:32737, and I believe I am having trouble because the original image is 16 bit- UINT16.

When I try and do the reprojection, I am ending up with an image that is all black, which suggests some possible truncation of the values at 255.

The code I am using is the same one from the Rasterio website. But I am not sure how to add the step of changing the datatype from Uint-16 to Uint-8 bit.

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = 'EPSG:32737'

with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    with rasterio.open('/tmp/RGB.byte.wgs84.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

I checked and the reproject function does not take the standard dtype argument. So I imagine this conversion has to happen in the src.meta metadata. But I was not sure if I can just change the setting, and that will take care of the conversion, or if I have to loop over the layers and change them, etc.

UPDATED:

So I wanted to show some of the output of this code. Mind that I tried this conversion with both a UINT8 image and a UINT16 image. The UINT8 image reprojection went just fine, but the UINT16 image had a problem. Here is an example:

This is a piece of the original raster drone image.

And here is the image after using this same code:


texts: Which one of these azimuth values is the right one? And why?

texts:
First of all, let's make clear that I want the Azimuth on the surface of the Earth, i.e. the angle between two locations, for example New York and Moscow.

I am testing some azimuth calculations with my JS functions (shown below). For the points A(-170, -89) to B(10, 89), I get ~90º.

JS function for Azimuth on sphere (from Wikipedia)

var dLon = lon2 - lon1;
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(dLon);
var angle = Math.atan2(y, x) * 180 / Math.PI;

JS function for Azimuth on oblate spheroid (from Wikipedia)

var dLon = lon2 - lon1;
var f = 1 / 298.257223563;    /* Flattening for WGS 84 */
var b = (1 - f) * (1 - f);
var tanLat2 = Math.tan(lat2);
var y = Math.sin(dLon);
var x;
if (lat1 === 0) {
    var x = b * tanLat2;
} else {
    var a = f * (2 - f);
    var tanLat1 = Math.tan(lat1);
    var c = 1 + b * tanLat2 * tanLat2;
    var d = 1 + b * tanLat1 * tanLat1;
    var t = b * Math.tan(lat2) / Math.tan(lat1) + a * Math.sqrt(c / d);
    var x = (t - Math.cos(dLon)) * Math.sin(lat1);
}
var angle = Math.atan2(y, x) * 180 / Math.PI;

In Calculator 2, I get 90º.

In PostGIS, I get 270º

In Calculator 1, I get 180º.

I know the Azimuth gets more and more distorted near the Poles, but that's exactly why I am testing at these spots. This variety of different solutions are confusing me. Could you please help me getting the right answer for this?


texts:
My setup:
Windows 10 64 bit (including all patches)
QGIS 3.0.0-Girona

My issue:
Installing numba not working

why:
I want to use GPU to process my python scripts

Details:
I followed the following:

OSGeo4W: typed "curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py"
OSGeo4w: typed "python get-pip.py"
OSGeo4w: typed "python -m pip install numba"
OSGeo4W: typed again "python -m pip install numba" and get...

Requirement already satisfied: numba in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages

Requirement already satisfied: funcsigs in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages (from numba)

Requirement already satisfied: singledispatch in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages (from numba)

Requirement already satisfied: numpy in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages (from numba)

Requirement already satisfied: llvmlite>=0.22.0.dev0 in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages (from numba)

Requirement already satisfied: enum34 in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages (from numba)

Requirement already satisfied: six in c:\progra~1\qgis3~1.0\apps\python27\lib\site-packages (from singledispatch->numba)

OSGeo4w: typed "python" and from prompt typed "import numba"
*no issues
QGIS python plugin: typed "import numba" and get...

      import numba
      Traceback (most recent call last):
        File "C:\PROGRA~1\QGIS3~1.0\apps\Python36\lib\code.py", line 91, in runcode
          exec(code, self.locals)
        File "", line 1, in
        File "C:/PROGRA~1/QGIS3~1.0/apps/qgis/./python\qgis\utils.py", line 664, in _import
          mod = _builtin_import(name, globals, locals, fromlist, level)
      ModuleNotFoundError: No module named 'numba'

I'm I doing something wrong here?  If I install numba in OSGeo4W, is that supposed to be represented in QGIS?

From OSGeo4W, from...

import numba
print(numba.__version__)

I get 0.37.0


texts:
In GeoServer (2.13.1) I'm trying to connect to Oracle DB 11g and 12c  and  connecting successfully both.  But in 12c I couldn't see my tables and can't add layer.

There isn't any error message or somethng like it.

Is there any way to reach them or how can i solve it.

-- Solved--

Oracle DB 12c has ojdbc8 driver in library file path  but GeoServer is tested and working correctly  with ojdbc7. With ojdbc8  successfully  connected and listed 11g tables but for 12c just connected to database, didnt list tables who in databese. Changing ojdbc8 to ojdbc7 connected and listed tables from 11g and 12c.


texts: Leaflet + GeoJSON and ArcMap : Map fails to load if missing variable

texts:
I am working on my first web map using the Leaflet API.

I have a shapefile that was converted to the GeoJSON format (within ESRI's ArcMap software). The Shapefile contains a field attribute named "Link_Ref" which contains a string (web address hyperlink) as a value. However, if there is a Link_Ref with no value (just blank space as the value) as shown in the attribute table in ArcMap, and then I export it to GeoJSON for Leaflet, Leaflet will fail to load and I get a blank white browser screen. I see in the browser console the following error:

No value provided for variable {Link_Ref}

So, if one value is missing, it seems to keep my web map from loading.

I am trying to clean the data, but unsure how to proceed. Will adding the text "null" to all the empty value fields allow my map to load? A quick test of this failed to work.

Why do I get this error and how to fix it?


texts:
So, I have an interesting problem... I have a PDF file which has over 100 pages and each of these pages is displaying a map in vector format.

By using some third party program I was able to find out that one of these layers is a parcel layer.

Now, what I am trying to do here is extract just this layer for every page in this PDF document and then stitch each of the parcel layer pages back together.

I am looking for some program or script which will automate some or all parts of this without screwing up the parcel line geometry. Additionally, I am assuming that I probable will need to manually align each of the extracted vector parcel pages which I guess is fine but what I think I need for this to work is some type of basemap background, I was thinking if one script can extract all the parcel vectors and then some other program can load all of them into some mapping program, me manually alligning each page and then finally being able to "join" the seperate pages into one large parcel vector map that the same (or yet another gis program) can append the correct lat/long data and turn it into a shapefile.

Uggghhh, so... Does anyone know of a solution for this?


texts: Creating link between geopackage and pdf

texts:
I'm trying to create a link between geopackage and a folder that contains pdf files. There are points with the same name and they should open the same pdf. I want name 100 to open 100.pdf. Name 101 should open 101.pdf

My folder looks like picture 2.

In QGIS 3 I've been trying to create this link under "actions" but don't understand how to connect this. Is it possible?


texts: Arcpy RuntimeError: An invalid SQL Statement on UpdateCursor

texts:
I'm trying to populate a field based on the values of another field. CalculateField was buggy, so I switched to da.UpdateCursor. This initially worked when I ran it as part of a script tool on ArcGIS Pro 2.4.3, but when I returned to this project after putting it aside for a few weeks and then making a few minor modifications to my code, ArcGIS Pro gave me an error:

 Traceback (most recent call last):
 File "Y:\AGUSPDB\States Procedures\Normalizer.py", line 126, in <module>
      for value in APN_Update:
 RuntimeError: An invalid SQL statement was used. [SELECT APN, Y:\AGUSPDB\03 States Compiled\Compiled_Ashley_AR_20_Parcels.gdb\Ashley_AR_Parcels, OBJECTID, Shape_Area, Shape_Length FROM Ashley_AR_Parcels]
 Failed to execute (Normalizer).

Incidentally, this code works in the Python IDLE shell, so I suspect the issue may have something to do with ArcGIS Pro. The offending code is prefaced with the comment "Problem begins here". Because the feature class's file path and other fields I never coded somehow ended up in the UpdateCursor SQL statement, I've added the a few other steps for context.

import arcpy

#Joins a standardized feature class template to the feature class the user loads.
FC_toNormalize = arcpy.SpatialJoin_analysis(inFeatures, normalizedTemplate, filepath,\
"JOIN_ONE_TO_ONE", "KEEP_ALL", "#", "HAVE_THEIR_CENTER_IN")

#Create list of fields to reverse-index from the template's attribute table
fields = [f.name for f in arcpy.ListFields(FC_toNormalize)]

#Problem begins here. fields[-13] is called "APN".
with arcpy.da.UpdateCursor(FC_toNormalize, [str(fields[-13]), "userSelected_field"]) as APN_Update:
    for value in APN_Update:
        value[0] = value[1]
        APN_Update.updateRow(value)

As for the minor modifications I made to my code before re-running the tool, they shouldn't have caused problems, but here they are:

Renamed a few variables
Reorganized the placement of some GetParameterAsText functions within the code (but not the index positions)
Added a new function a few lines below the for loop in which this error occurs

Even though the code works in the Python shell, I'll obviously have to modify it so ArcGIS Pro to digest it. Any thoughts?


texts: Postprocessing with NMEA-data if no rinex-data is available

texts:
I have an UAV which does not give me the raw sensor data (e.g. rinex or ubx). The only format I get is NMEA. Now I want to know how could I use the NMEA-data to get the highest precision of my postitions via postprocessing.
Are there any known methods to do this?
I thought of something like a virtual sensor at the positions from the NMEA-data to maybe get the rinex-data afterwards from a service like SAPOS.


texts: GDAL: can I rasterize a layer's FIDs without having an explicit field id column?

texts:
I would like to rasterize a layer by FID without adding an FID field to the layer.  Is this possible with GDAL?  Something like like this:

raster = gdal.Open(raster_path, gdal.GA_Update)
vector = ogr.Open(vector_path)
layer = vector.GetLayer()
gdal.RasterizeLayer(
    raster, [1], layer, options=["ATTRIBUTE=FID"])

But of course this doesn't work because FID is not defined.  I've tried using layer.GetFIDColumn() but I get an empty string back.  That's not surprising since this is the result if  GetFIDColumn is not supported.

Otherwise, my approach is to create a new DataSource with an FID field only.  Copy the original's layer geometry to this one while updating the FID value.  Rasterize that layer, then delete it.  I'd like to avoid this solution since in general I don't have control over the size of geometries coming in.  It seems wasteful to create a whole copy to add an FID field.  And dangerous to add a field to the original input only to delete it later.


texts: Identifying clusters of points within a given distance

texts:
So I have a database of lats/longs representing incidents occurring on or near roadways.

My goal is to identify the street block areas with the highest number of incidents in the database. I don't have any other map layers I'm using for this.

What I'd like to do is, for each point, see how many other points are a given defined distance away (a quarter of a mile, for instance). The result would then show where the clusters are.

Any recommendations on the best way to do this? I have copies of both QGIS (with the spatial query plugin) and ArcMap, although I'm much more familiar with QGIS.


texts:
I have a polyline that I need to make some gaps in at some locations, by cutting this part inside the red box, and leaving a gap in this particular location for instance.
How can I do that in QGIS?


texts: Where can I find better data than yahoo geoplanet or factual's world geographies?

texts:
Currently I am using Factual's world geographies api to do a location search so if you search london, london england should appear or if you search new york city, new york city will appear.

I am having trouble because there is no filters like population that I can use to sort by the best match.

Here is factuals data. I found a download for geoplanet here but I can't seem to find any data set with these sorts of filters.

Do I need to find a better data set or am I just not structuring the query correctly?


texts: How to run equivalent of gdal2tiles using GDAL C/C++ API?

texts:
GDAL includes a comprehensive low-level API in C and C++.  It also includes some high-level APIs that effectively run the equivalent of some of the GDAL command-line utilities.  (See:
https://gdal.org/development/rfc/rfc59.1_utilities_as_a_library.html ).
The supported command-line utilities that can be effectively run directly in C/C++ using this API include gdalwarp, gdalinfo and gdal_translate.  However, gdal2tiles is not available this way (presumably because it's a Python script, rather than a C/C++ utility program).  Additionally, because I'm runnning on iOS, I do not have access to simply run the actual gdal2tiles utility itself (eg, with a system call).
Using the GDAL C APIs, what is the best way to run the equivalent of gdal2tiles?
Or is my only option to look up the actual source code of the gdal2tiles utility itself, and translate the relevant portions into C code?  (Would be waaaaay overkill if there is any other way!)
My goal is to produce a TMS directory, including TMS metadata, from a user-supplied image file and user-supplied (interactively generated) georeferencing information.  Or something similar that can be managed by GDAL and read by OpenLayers.  (I did briefly consider MBTiles instead. However, the GDAL driver for MBTiles specifies that it only works for pseudo-mercator and I need to work with all SRSs.)
(I did something similar without GDAL a few years ago, including pyramid tiles. But was hoping to come up with a more standardised, flexible and maintainable method. I may have to look up my old Objective-C code and see if it is applicable to my current situation, and maybe translate it into Swift.)


texts: How can I add the outline of Europe countries in this figure in R?

texts:
Here is my code to create this distribution map of Europe in R. So how can I add the outline in this figure to divide different countries and also add the countries' names in relevant regions?
>library(sf)
>library(rgbif)
>all <- st_read('~/Desktop/amre2935_supplementary_atlas_s4_shapefiles_reptiles/allcountries_rep.shp')
>all <- subset(all, select='emy_orb')
> summary(all)
    emy_orb                geometry
 Min.   :0.0000   POLYGON      :5453
 1st Qu.:0.0000   epsg:4326    :   0
 Median :0.0000   +proj=long...:   0
 Mean   :0.1423
 3rd Qu.:0.0000
 Max.   :1.0000
> all
Simple feature collection with 5453 features and 1 field
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -14.3358 ymin: 34.30561 xmax: 66.7329 ymax: 71.65074
geographic CRS: WGS 84
First 10 features:
   emy_orb                       geometry
1        0 POLYGON ((-14.3358 47.4027,...
2        0 POLYGON ((-13.6733 47.39694...
3        0 POLYGON ((-13.01106 47.3873...
4        0 POLYGON ((-12.34923 47.3739...
5        0 POLYGON ((-14.29924 50.1012...
6        0 POLYGON ((-14.30573 49.6515...
7        0 POLYGON ((-14.28579 51.0004...
8        0 POLYGON ((-14.29261 50.5508...
9        0 POLYGON ((-14.2788 51.45003...
10       0 POLYGON ((-14.32421 48.3023...
>pal <- c('grey', 'red')
>plot(all['emy_orb'], key.pos=NULL, reset=FALSE, main='All Countries', pal=pal, border=NA)


texts: Splitting feature into fixed measurements

texts:
For my thesis I'm trying to split featuring in QGis into multiple new features of a fixed size (5.5 by 2 meters). I think it would be the equivalent of ArcGIS' Parcel Division tool, but unfortunately I don't have the ArcEditor/ArcInfo license. I already downloaded some plugins but I can't manage to find the solution. Hope you can help!

I've uploaded some pictures I edited with paint to visualize what I mean. NB: they're not only straight squares, but also parallellograms.

[3


texts: Clipping/cutting more than one polygon layer in one step using ArcGIS for Desktop?

texts: currently I am working with ArcGIS 10.4. My project contains several layers such as streets and buildings, as well as retail and residential areas for a large-area state in Germany. As this is a lot of data and I am currently focusing on a smaller region within this state, I am looking for an easy way to “cut” all data outside my region of interest, but retaining the symbols and the labels. To sum it up: Is there a way to do something like zooming into my project and cutting everything else outside in one step?


texts: Defining a cell style of an attribute table in a layout using PyQGIS 3

texts:
I have an attribute table in my layout, defined as follow :
layout = QgsPrintLayout(project)
layout.initializeDefaults()
QgsProject.instance().layoutManager().addLayout(layout)

pdf = QgsProject.instance().mapLayersByName('plan_de_ferme')[0]
pdf_table = QgsLayoutItemAttributeTable.create(layout)
pdf_table.setVectorLayer(pdf)
pdf_fields = ['champ','sup','supepnd']
pdf_table.setDisplayedFields(pdf_fields, True)

I would like to set a background color for the header row and all even rows, as follow :

I can use pdf_table.setBackgroundColor(QColor(166, 206, 227, 255))  to apply a blue color for the whole attribute table but it's not what I am looking for.

I think I should use QgsLayoutTable Class and QgsLayoutTableStyle Class. According to its definition, the class QgsLayoutTable has a method called setCellStyle  that requires two inputs :
void setCellStyle( CellStyleGroup group, const QgsLayoutTableStyle &style )
CellStyleGroup is a public type of QgsLayoutTable, defined ad follow in the first link:

I don't know by what to replace "const QgsLayoutTableStyle &style"
I tried
pdf_table.setCellStyle(qgis.core.QgsLayoutTable.CellStyleGroup.EvenRows,
                        cellBackgroundColor(QColor( 166, 206, 227, 255 )))

but I get the following error: "NameError: name 'cellBackgroundColor' is not defined",
and
pdf_table.setCellStyle(qgis.core.QgsLayoutTable.CellStyleGroup.EvenRows,
                       QgsLayoutTableStyle().cellBackgroundColor(QColor( 166, 206, 227, 255 )))

but I get the error "TypeError: 'QColor' object is not callable",
and
pdf_table.setCellStyle(qgis.core.QgsLayoutTable.CellStyleGroup.EvenRows,
                       QgsLayoutTableStyle.cellBackgroundColor(QColor( 166, 206, 227, 255 )))

returns "AttributeError: 'QgsLayoutTableStyle' object attribute 'cellBackgroundColor' is an instance attribute".
Like I'm not familiar with the several Python definitions, I don't understand what these means.


texts: QGIS Pre-set CRS/projections missing/unavailable

texts:
I am using QGIS 2.0.1 running off of a Ubuntu computer.

I have been going through some basic tutorials (http://fossgeo.org/free-qgis-courses/) and have come up against a small and silly little obstacle. Some of the problems and examples involved in these tutorials require specific CRS presets. When I first started working, my copy of QGIS had all of the different CRS types that came up. However, after a few days something appears to have happened and I no longer have access to a bunch of the CRS types referenced in the tutorials. I have tried looking them up by name and EPSG# to no avail. I have tried reinstalling the program.

I have not been able to find a download or plugin to restore all of the CRS presets, does such a thing exist?


texts:
I have two stars-objects and I simply would like to calculate the mean for each pixel.
I am not really sure how to do that and if my approach is correct. At the moment I'm doing the follwing:
I have
> jan
stars object with 2 dimensions and 1 attribute
attribute(s):
  mean_no_na
 Min.   :0.27
 1st Qu.:0.82
 Median :1.01
 Mean   :1.15
 3rd Qu.:1.44
 Max.   :3.69
 NA's   :194505
dimension(s):
  from  to  offset delta  refsys point values
x    1 641  605720   250 unnamed    NA   NULL [x]
y    1 643 5220237  -250 unnamed    NA   NULL [y]

and
> feb
stars object with 2 dimensions and 1 attribute
attribute(s):
  mean_no_na
 Min.   :0.28
 1st Qu.:0.64
 Median :0.80
 Mean   :0.85
 3rd Qu.:1.06
 Max.   :1.93
 NA's   :194505
dimension(s):
  from  to  offset delta  refsys point values
x    1 641  605720   250 unnamed    NA   NULL [x]
y    1 643 5220237  -250 unnamed    NA   NULL [y]

and now I'd like to simply compute the mean. The two two rasters cover exactly the same geographic extent and are just for the month of january and february. For that I'm doing:
c = c(jan, feb, along=3)
mean_two_month = st_apply(c, c("x", "y"), mean)

But I'm not really sure if this is the correct approach.


texts: I have a set of shapefiles (Kabul, Afghanistan) that are currently in the WGS84 geographic coordinate system (GEOGCS), which I believe has a SRID of 4326.  In ArcMap, this coordinate system is decimal digit based.  How would I determine the appropriate meters-based coordinate system for this set of shapefiles, and thereby a meters-based SRID for this set of shapefiles?


texts: Importing shapefiles to a postgis on a remote server

texts:
I have two servers X and Y. X contains a postgis instance whereas Y doesn't have postgis installed, but connects to the postgis instance on X through a connector in nodejs.

I want the node script to import some shapefiles which is normally done with the pgsql2shp command. However since pgsql2shp isn't installed on server Y I can't just call the command.

Is there a way to call the pgsql2shp command (or do it's equivalent) through a query? Or should I just have the script on Y ssh to server X and run the command?


texts:
Is there a way to report a point layer's labels from the map to the legend in the Print-Composer ?

I have a quite dense point layer for which I need to display names. They wont fit on the map, so I though I might just assign them numbers, and display their names alongside the numbers in the legend.

Is there a way to achieve this in QGIS (currently using 2.18.10) ?


texts:
I am busy building a web map app using boundless SDK.
I am trying to implement a file up-loader, however, the majority of examples use
new Ext.ux.form.FileUploadField
When I try to use it in my project it says that the xtype is not recognised.

I am currently using this:

    xtype : 'textfield',
    inputType:'file',

But I find it cumbersome, ugly and have to add allot of extra code to make sure validation fires when file is selected etc...

How can I make the FileUpload UX plugin part of my project?

I am quite new to the environment.


texts: Can't find MapFrame of GraphicsContainer in PageLayoutControl (C#) ?

texts:
I am getting a problem: I added element (text, line) in a ActiveView.FocusMap of a PageLayoutControl, then all elements are showed in the MapFrame of PageLayoutControl.GraphicsContainer.

If I don't select elements (in MapFrame) I can find the MapFrame by the loop:

IGraphicsContainer graphicsContainer = myPageLayout.ActiveView.GraphicsContainer;
graphicsContainer.Reset();
IElement pElement = graphicsContainer.Next();
 while (pElement != null)
            {
                if (pElement is IMapFrame)
                {
                      // do something
                 }
            }

But, when one element is selected by using "Select Elements" command, the program won't enter the "if" function (can't find MapFrame).

Are there anyone know how to get The MapFrame in this case ?

It's very urgent

Thank you very much !


texts:
I have a lot of polygons which look like this:

In QGIS geometry checker there is no error even with a tolerance of 10^-3. But in an Online-ArcGIS System (10^-3 tolerance) it seems that it detects an self-intersection at this point. Is there any possibility to find this points in QGIS?

The files are Shape-Files from an external Programm. I have to take them like they are. So the correction have to be in QGIS. Most Polygons have more that 60k Points and in the shapefiles are sometimes more that 50k features.


texts:
Trying to perform some Linear Spectral Unmixing on a LandSat 7 scene. I have 6 bands, and after MNF-transformation have kept 4 of them.

I want to perform a Vegetation, Impervious surface, Soil (V-I-S) based classification. The end goal is to generate a map, with each pixel showing the percentage of each pixel that is predicted to be impervious. These will be sampled from 3x3 (90x90m) pixel squares that I have randomly distributed throughout the study areas.

Thus far, I've had limited success. My best result has had an R^2 of 0.64, when compared to the ground-truth digitized 90mx90m squares. I've achieved this by using the Scatter Plot function of MNF Band 4 vs. MNF Band 3, and selecting groups of pixels based upon their appearance on a ground truth image. Here's what the scatter plot that resulted in the best R^2 looked like: https://drive.google.com/open?id=0B-qDNCi6oOLIcUxseG9OYzIyRkk

I've tried using the Spectral Hourglass Wizard, but it stops during the n-D Visualization stage. The error warning I get is: "SVDC: Argument A must be a 2-D matrix: MEANS". When I hit "OK", ENVI aborts.


texts: Can anyone tell how to reproject VIIRS products (interested VNP09 and VNP14IMG) into standard projection and reference system? I managed to process some of products like CLDMASK_L2 and even VNP14 (750 m) using SNAP and HEG tool. But some of products cannot be processed in the mentioned software despite of georeferencing metadata.


texts: How to access webmap in ArcGIS Portal from Angular app on localhost

texts:
I'm testing an Angular CLI application that's running on my PC at http://localhost:4200 and which needs to access a webmap hosted on an ArcGIS Enterprise Portal.

The typical method for doing so within the ArcGIS Server JavaScript API is:

EsriConfig.portalUrl = "https://url-to-portal";
const webMap: esri.WebMap = new EsriWebMap({
    portalItem: {
      id: "webmap-item-id"
    }
});

This is failing for me with a CORS error because the Portal security has been restricted to a particular domain:

  Access to fetch at 'https://url-to-portal/web/sharing/rest/portals/self?f=json&culture=en-gb' from origin 'http://127.0.0.1:4200' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource. If an opaque response serves your needs, set the request's mode to 'no-cors' to fetch the resource with CORS disabled.

How can I access this webmap from my PC during development?

I tried to add my localhost to the list of approved domains on the Portal security page using all combinations of localhost, http://localhost, http://localhost:4200, http://127.0.0.1 etc but they all return an error message from the Portal security page:

Alternatively, how can I try the "set the request's mode to 'no-cors'" suggestion within the ArcGIS Server JS API syntax?


texts: MapBox GL JS Draw plugin: how to modify point dimensions?

texts:
I need to modify the point dimensions in an application that use MapBox GL JS Draw plug-in (rif. https://github.com/mapbox/mapbox-gl-draw/).

I've tried to produce a little example with this code:

    <meta charset='utf-8' />
    <title></title>
    <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link href='https://api.tiles.mapbox.com/mapbox-gl-js/v0.24.0/mapbox-gl.css' rel='stylesheet' />

<link rel="stylesheet" type="text/css" href="https://api.mapbox.com/mapbox-gl-js/plugins/mapbox-gl-draw/v0.11.15/mapbox-gl-draw.css" />

This code works fine but the points dimensions are too small and I'd like to try to increase them.

I've tried to use a new definition of the "draw" object like this ...

var draw = mapboxgl.Draw({
  drawing: true,
  displayControlsDefault: false,
  controls: {
    'point': true,
    'trash': true
  },

  styles: [
    {
     'id': 'points-are-blue',
     'type': 'circle',
     'paint': {
      'circle-radius': 7
     }
    }
  ]
});

It works but I lose the original styles ....

Suggestions / examples / useful link are useful!


texts:
I am trying to to write a custom feeder tracing routine, but facing some issues.
I am using the following ArcFM method:

IMMFeederTracer.TraceFeeder(IGeometricNetwork, IMMFeederSource)

MMFeederSource is simply initiated as new MMFeederSourceClass();.
There does not seem to be any QI from IRow and IFeature to set the source.

In the current form, the Feeder Manager starts but after a while throw an HResult 0x80044020 error.

Any ideas how to initialize the feeder source?
I could not find anything on the ArcFM Resource Center.

Thanks in advance
Ujjwal


texts: Resetting measure tool units back to decimal degrees in ArcMap?

texts:
I was using the measure tool and it was displaying in decimal degrees as I wanted.  Then I changed it to meters temporarily and then wanted to go back to decimal degrees, but no matter what I do I can't get it back.  I followed these instructions from Arc desktop help online:

Specifying a coordinate system—Help | ArcGIS Desktop

and my data frame is set to display in decimal degrees, but the measure tool is still displaying meters.  What am I missing?


texts: Feature property uploads as null- Google Earth Engine

texts:
I am attempting to upload a feature collection which has a SURVEYDATE column. However, when I ask for the property, it returns null.
var table = ee.FeatureCollection("users/shopvac18/zm/NAS_ZM_2")

print(table.first())

// ee.Feature with a series of properties, including "SURVEYDATE"

print(ee.Feature(table.first()).get("SURVEYDATE"))

// null

Here is a link to a small example
Link to the original csv.


texts:
Let's say we have the following line layer.

Why is that when we perform routing, it doesn't go from point A to point B like below:

Instead, it goes in the other path. I've sent the one way direction to false, yet it still chooses the other longer path. I'm using QGIS


texts:
After a day of unsuccessfully searching, I ask for help -how can I use symbol levels in PyQGIS plugin . Here is my code that display the layer:
layer = self.iface.activeLayer()

symbol =QgsLineSymbol.createSimple ({'color':'black', 'width':'0.8', 'line_style':'solid'})
symbol_l2 = QgsSimpleLineSymbolLayer.create ({'color':'white', 'width':'0.5', 'line_style':'solid'})
symbol.appendSymbolLayer (symbol_l2)

rendererr = QgsSingleSymbolRenderer (symbol)
layer.setRenderer(rendererr)

In QGIS DESKTOP, you can see the intersection of roads clean without lines,
after executing the following commands :
layer properties -> symbology -> line -> Advanced -> Symbol Levels -> Enable symbol level`


texts:
I want to calculate distance between smaller dots (they represent towns, villages etc.) and these larger red dots (capital cities of admin. region).
My goal is to show wich location is the furthest from the capital city for each region.

What tool/plugin is best to do it? I was trying with 'HUB distances' but it doesn't work properly I think.


texts: Clipping raster python script stops abruptly?

texts:
I have created a tool in ArcGIS based on a python script that will clip a large raster dataset into smaller tiles based on a fishnet polygon feature class. The script iterates through each feature in the fishnet and uses the selected feature to clip the raster. For some reason it keeps stopping after the 58th feature (there are 89 total). I tried using it on a different raster but with all the same parameters and it stopped after the 53rd feature. Any idea why this is happening?

Here is the pertinent code

for feat in fishnet:
    arcpy.SelectLayerByAttribute_management("fishnetlayer", "NEW_SELECTION", '"FID" = ' + str(select))
    arcpy.Clip_management(inRaster, "#", "BE_Seg_" + str(tile) + ".png", "fishnetlayer", "0", "ClippingGeometry")
    arcpy.AddMessage("Tile " + str(tile) + " successfully created")
    rastnodata = arcpy.GetRasterProperties_management("BE_Seg_" + str(tile) + ".png", "ALLNODATA")
    rastempty = rastnodata.getOutput(0)
    arcpy.AddMessage(rastempty)
    if rastempty == "1":
        arcpy.Delete_management("BE_Seg_" + str(tile) + ".png")
        select = select + 1
        tile = tile + 1
    else:
        select = select + 1
        tile = tile + 1


texts: Python and Data driven pages - Advanced use

texts:
I need to make a number of maps where only one input changes (a layer with the range of a radio tower) and I want to add the layer, apply symbology, zoom to it, change the title to include the name of the tower and then export it as a pdf with the tower name.

Rather than adding the layers and applying symbology I can just put all the towers in a Group Layer and then the script can just turn all but one of them on, zoom, print and then turn it off and continue to the next tower.

One long way is using: Data Driven Pages --> Page Definition option in the Definition Query tab of a Layer's Properties and specifying this but I was wondering if there was an easier work around?

The issue with using Page Definitions is that some towers overlap and it means that I have to add in a column which will have have the tower name and signal strength repeatedly for each row. Given that we have more than 30 towers, even in a model this is laborious process that has to be then repeated when any base data changes.


texts: Split a feature when intersecting with a feature of another layer using PyQGIS/Python?

texts:
I have a buffer layer (green polygon) which I want to split to two polygons whenever it crosses a barrier (blue line). I have been trying to use "splitGeometry" method, but I just can't get it to work. My code so far is this :

while ldbuffprovider.nextFeature(feat):
  while barprovider.nextFeature(feat2):
    if feat.geometry().intersects(feat2.geometry()):
        intersection = feat.geometry().intersection(feat2.geometry())
        result, newGeometries, topoTestPoints=feat.geometry().splitGeometry(intersection.asPolyline(),True)

Which returns 1 for result( error) and an empty list for newGeometries.
Any help is greatly appreciated.


texts: What causes RuntimeError: ERROR 010240 saving after CellStatistics?

texts:
I am having trouble saving outputs from CellStatistics but only when I run it on a list of rasters. This was working on python 26 but not 27.

I know this example is nonsense, but:

env.workspace =    "G:\\Faculty\\Mann\\Historic_BCM\\Aggregated1080\\"
relevant_tifnames_year_i= ['cwd1975jan.tif', 'cwd1975feb.tif', 'cwd1975mar.tif', 'cwd1975apr.tif', 'cwd1975may.tif', 'cwd1975jun.tif', 'cwd1975jul.tif', 'cwd1975aug.tif', 'cwd1975sep.tif', 'cwd1975oct.tif', 'cwd1975nov.tif', 'cwd1975dec.tif']

mean = CellStatistics(relevant_tifnames_year_i, "MEAN", "DATA")
mean.save(str(env.workspace)+"outa.tif")

sd = CellStatistics(relevant_tifnames_year_i, "STD", "DATA")
sd.save(str(env.workspace)+"outb.tif")

Atest= CellStatistics([mean,sd], "MEAN", "DATA")
Atest.save(str(env.workspace)+"outc.tif")

Traceback (most recent call last):
  File "<pyshell#81>", line 1, in <module>
    Atest.save(str(env.workspace)+"out3.tif")
RuntimeError: ERROR 010240: Could not save raster dataset to G:\Faculty\Mann\Historic_BCM\Aggregated1080\Scratch\meanc_ras3 with output format GRID.

Notice that I can save CellStats run on a list of files from my working directory, but not on a list made up of my intermediate inputs "sd" and "mean". I know that I can run it if I pass the file paths, but in my actual code I need to run it off of intermediate inputs. Am I missing something here?

Note: Overwrite = True has been set. There are no files with the same name in the output folder etc.

Here is the actual code (storing rasters to a list and then running CellStats)

for year in base_year:
    # paths for files
    relevant_tifnames_year = [v for v in relevant_tifnames if int(re.sub('[^0-9]','',v) or 0) >= year and int(re.sub('[^0-9]','',v) or 0) <= year+years_to_end ]
    #list to store annual means
    decadal_holder_Amn = []

    # for each period of interest
    for year_i in range(year, (year+years_to_end+1)):
       # get relevant paths
       relevant_tifnames_year_i = filter(lambda x: str(year_i) in x,relevant_tifnames_year)
       #calc annual mean
       mean = CellStatistics(relevant_tifnames_year_i, "MEAN", "DATA")
       #append mean mean to list of annual means
       decadal_holder_Amn.append(mean)
    #calc mean of annual means for the decade
    Amean=CellStatistics(decadal_holder_Amn, "MEAN", "DATA")
    Amean.save(str(env.workspace)+"outd.tif")


texts: Changing sequence of point features into order along line using ArcPy?

texts:
The picture below shows that the order of attribute record for point features is different to their actually distributed location.

How can I change or sort the attribute records into a proper order by using Python for ArcGIS?

There also a line under the point features.


texts:
I'm new to ArcGIS and am trying to determine sites that meet a series of geographic criteria including: within various distances of certain features, as well as within a distance of one layer but not within that layer, and within a distance or inside of the last layer.
I'm having trouble figuring out the last two criteria right now without getting too complicated with the geoprocessing tools. So far I've just buffered the layers where I just need to be within a distance of them, and then I intersected these to cut down on potential sites. I tried to use the Select by Location tool after with selection method, "select features from", where target layer = the intersected layer, and using the layer with the bold criteria as the source. I'm not seeing any options for the spatial selection methods that can help me; the closest one seems like, "are within a distance of the source layer feature"
Also, if anyone could tell me why my intersection from buffered layers doesn't have any records in the attribute table even though I can see the different polygons on the map - I was confused on that, too.


texts: Store GPS trace with timestamps, then interpolate location by time


texts: Recommended NoData value for DEM/DTMs created using gdal

texts:
Background info: I've been generating some raster DTMs from vector points using gdal_rasterize, and was wondering if there is any consensus as to what is the best NoData value to use. The Python GDAL script I'm using appears to default to 0, but obviously that might be problematic if my country-wide dataset has areas close to sea level.

Sample of script I'm using:

try:  from osgeo import ogr, osr, gdal
except: sys.exit('ERROR: cannot find GDAL/OGR modules')

nztmSRS = osr.SpatialReference()
nztmSRS.ImportFromEPSG(2193)
xRes, yRes =(5,5)  #Spatial resolution defaults to 5 (metres)

src = "Points.shp"
dst = "Output.tif"

gdal.Rasterize(dst, src, outputSRS=nztmSRS.ExportToWkt(), xRes=xRes, yRes=yRes, noData="0",attribute='ELEVATION')

I've seen NoData values out there choosing -9999, -3.4028230607371e+38, etc. There are probably valid arguments for using any of those, and I'd like to have a feel of A) what the NoData options are and B) what seems to work better in different use cases. E.g. when:

The output will need to be reprojected using gdalwarp
The outputs will need to be mosaicked together with other geotifs
The output will need to be loaded into Spatialite/Postgis down the line


texts: Polyline (shp) layer not displaying correctly in QGIS

texts:
I'm a new QGIS user. I'm trying to convert a polyline vector layer to polygon. I tried to it in ArcMap 10.1 but it misses a lot of polygons. So, I'm attempting to this in QGIS running the process "Lines to Polygons". However, when I load my shapefile into QGIS, it doesnt' display well. Only half of my map and some random spots are shown (I also noticed that in QGIS the attribute tables has way less lines).
I attach two screenshots from the same layer in ArcMap and QGIS.
Thanks in advance,

JF


texts: Flow Direction created from depressionless DEM inaccurate in flat areas

texts:
I am trying to build a flow direction raster for an area from a 1m DEM created from LiDAR. I need this in order to build a flow accumulation raster, so I can use that to build a stream network. I have followed ESRI's guide (found here http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?id=3914&pid=3910&topicname=Creating_a_depressionless_DEM) on how to create a depressionless DEM so that flow direction is accurate. This works very well except in some areas that are quite flat. In these areas the Fill tool used in the depressionless DEM workflow fills in the streams on the DEM. This results in an inaccurate flow direction because the streams are filled to the same elevation of the surrounding area. As a result, flow accumulation is inaccurate in these areas and as a result of that the stream network built from it is as also inaccurate in these areas.

Inaccurate Flow Direction made from depressionless DEM with inaccurate streamlines (black lines) created from resulting inaccurate flow accumulation (shown to display inaccuracy of flow direction):

1m DEM with clearly visible streams. The blacklines are the stream lines created from the resulting inaccurate flow accumulation. Notice how they are very inaccurate compared to where the stream actually is on the DEM. This is due to the Fill tool filling in the stream, causing the flow direction to not flow into the right area.:

Does anyone know of a way to create a depressionless DEM without inadvertently filling in streams?


texts: Extracting polygons touched by lines using ArcGIS Desktop

texts:
I am using ArcGIS Desktop 10.5.

I made an hexagonal grid layer and also I have a road layer (lines). The end goal is to calculate every hexagon's distance to particular variables (sport facilities, commercial places, important squares, trees, etc) and ponderate them to get a "problematic" index. The closer to those variables, the more problematic the area in the hexagon will be considered. Then, via an Spatial Join (unless I find a better tool), that index will be transferred to the line layer where every segment will get as a value being the mean of the index in the hexagons it touches/crosses.

The thing is, the grid has 200 square meter hexagons that add up to 481338 polygons, and the road layer has (for now, most likely to be reduced) about 7600 segments. My PC isn't very powerful so I would like to reduce the number of elements to be processed, especially since many of the hexagons will be worthless as they are not touched by any road, and work with only those polygons that I know will be able to transfer their index to the road layer.

I already achieved that using Select by Location with intersect, but I would like to find a tool that, with the right parameters, does this without the need of manual intervention. I have intention of creating a model so that this process can be executed in further analysis without having to do extra work, just entering the initial layers and clicking Run. Also I don't want to increment the hexagon's area to reduce their number (I might even reduce the area to have more hexagons).

I also thought that there might be a parameter in the Generate tessellation tool, or another different tool, that creates a grid consisting only of polygons along a line, those being the roads. Still haven't found this yet.

I should add that I'd like to be the grid to be hexagonal, as i consider its shape to be better for this particular task. I know Grid Index Features tool creates the grid I want, but it only creates squares.

The thing is, I know how to select the overlapping features and then copy them to another feature/layer with a process consisting of various steps, I even found a new tool to add to the toolbox that lets me create the exact thing I need. I already got exactly what I wanted, but I still have the doubt if there is a way within, lets call it "vanilla", ArcGIS to do this exact same thing with just one step, one single tool that you input your roads and creates these hexagons along them. Without having to add new tools or needing more than 1 step and manual intervention.

Does this tool exist?

Something you can show to a person who barely knows ArcGIS and say "you execute this, this happens, period".


texts: I want to import 100 tif images into QGIS. For each tif I have to click 'OK'. Is it possible to batch import all 100 without clicking ok every time?


texts: Domains can only be delete or modified by the owner?

texts:
Arcmap States that "When a new domain is created, the owner of that domain—that is, the user who created it—is recorded. Only the owner of an attribute domain can delete or modify it."

So say for example a usercreates a domain, then this user leaves the organization. Would it be impossible to then edit or delete that domain?


texts:
I am new in ArcMap while using R Crisis software.

I need to convert a shapefile into a .grd file.

I read back the previous questions and answer and tried to execute and pop up that it was successfully saved as .grd.

But why it does not appear as .grd file, all files are in type .adf?


texts: Using gdalwarp to remove one row of pixels or add one row of pixels

texts:
I am trying to use GDAL_calc but my two rasters are off by one pixel in the y dimension.  One is 256005 the other is 256004 so GDAL_calc fails with a multiplication.

Is it possible in gdalwarp to keep the same pixel resolution but to remove the extra row in the one with 256005 or add a row of 0's (or nodata) to the one with 256004?

  C:\temp>gdal_calc --co="COMPRESS=LZW" -A a2000mfw.tif -B wdpa_eq_1_2000.tif -
  tfile=mfw_pa_00.tif --calc="A*B"

  > Error! Dimensions of file wdpa_eq_1_2000.tif (1295778, 256005) are
  different from other files (1295778, 256004).  Cannot proceed


texts: Tearing of polygons using ggmap and readOGR

texts:
I'm trying to incorporate a shapefile of US states with ggmap and keep getting "tearing" of my polygons.  I tried changing group=id to group=group which helped, but didn't seem to resolve the boundary.

library(rgdal)
states <- readOGR(dsn="shapefiles", layer="states")
proj4string(states)
\[1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"
states <- spTransform(states, CRS("+proj=longlat +datum=WGS84"))
states <- fortify(states)
sstates <- get_map(location = c(-81, 35, -69, 45), zoom = 4, maptype = "watercolor")
sstates <- ggmap(sstates)
sstates <- sstates + geom_polygon(aes(x = long, y = lat, group=group),
                                  data = states, color ="white", fill ="orangered4",
                                  alpha = .4, size = .2)

Original shapefile

https://www.arcgis.com/home/item.html?id=f7f805eb65eb4ab787a0a3e1116ca7e5


texts: ArcGIS 10.0, Field Calculator Nested "if" statement with strings

texts:
I need to assign a value depending on other two columns "FNP" and "Teilort". I wrote this script in phyton for field calculator. It actually works (no error messages) but the results in every cell is 99.

I checked the two columns "FNP" and "Teilort" and the problem doesn't look like to com e from there (no extra blank or white spaces).

What can be the problem?


texts: Create KML from MXD that contains map service layers like feature layers

texts:
How can I create a KML (tools suggestion), provided that I have a map (.mxd) file which layers are service layers of hosted services and not local data coming from a geodatabase.
I know that the map to kml tool works only for local data and not with data coming from services.
Any ideas ?


texts: Open LiDAR DEM of the UK with *.asc extension

texts:
I've downloaded a tile of a freely accessible LiDAR DEM of the UK from http://environment.data.gov.uk/ds/survey#/download.

But I don't know how to open it. The files inside have the *.asc extension.
I can see that this extension is not among those supported by Gdal. And, neither Sentinel Application Platform SNAP nor ENVI did manage to open it.

Which software is used to open these DEM files?


texts: Service or way to input a table of lat long and return an output of street addresses?

texts:
The following website can be used to enter a single pair of latitude and longitude coordinates and return a street address:

https://www.latlong.net/Show-Latitude-Longitude.html

I want to get street addresses for a table of lat long coords. Any ideas?


texts: Counting occurrences across many rasters using ArcGIS Desktop?

texts:
I am interested in drought occurrences in the US, especially how often droughts of a particular severity occur in particular regions.  I’ve got nice vector maps from http://droughtmonitor.unl.edu showing drought severity for each year (scale of 0 - 4, with 0 = minor drought and 4 = extreme drought) for 2006 - 2017. I’ve converted the vector maps to rasters, all with the same extent and cell size.

Using ArcMap 10.5 (and prob Spatial Analyst), I want to count how many times a drought index of each level has occurred in a given region (i.e. raster cell) over the 10 year period.  I’m imagining inputting the 10 yearly maps and ending up with 5 output rasters, each holding the count of drought level for each cell. So far I haven’t found a way to count occurrences across multiple rasters.  I don’t want sum or weighted average, I want counts.  Any suggestions?


texts: Connecting to SQL Server Geodatabase using ArcObject in C#

texts:
I am going to save a polygon that is created from points to remote geodatabase. I am new with ArcObjects and can't understand how it works... at this point I want to open a connection to database using ArcObjects and save geometry in it. I am using this code:

public IWorkspace open_ArcSDE_Workspace(string server, string instance, string user,
string password, string database, string version)
{
    // Create the workspace factory.
    Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.SqlWorkspaceFactory");
    IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance
        (factoryType);

    // Create the connection properties.
    IPropertySet connectionProps = new PropertySetClass();
    connectionProps.SetProperty("dbclient", server);
    connectionProps.SetProperty("serverinstance", instance);
    connectionProps.SetProperty("authentication_mode", "OSA");//The type of authentication to use. Valid values are DBMS and operating systems authentication (OSA). DBMS is the default mode and is not required.
    connectionProps.SetProperty("user", user);
    connectionProps.SetProperty("password", password);

    // Open the workspace.
    //IWorkspace workspace = workspaceFactory.Open(connectionProps, 0);
    return workspaceFactory.Open(connectionProps, 0);
}

But it throws errors:

enter code hereError HRESULT E_FAIL has been returned from a call to a COM component.

Can you help me and give me some examples that I can understand... The only source I have found for ArcObjects is ESRI and I couldn't find ant solutions for it.


texts: ArcGIS Server Data Extraction - Job does not exist or is inaccessible

texts:
I have followed the documentation and video to set up a data extraction geoprocessing service, which I am accessing with the ArcGIS Viewer for Flex (v. 2.3.1).  When I run the Extract Data widget in the viewer I get an error: Job '(job folder name)' does not exist or is inaccessible.  When I go to the arcgisjobs folder on the server the folder is there, along with the scratch folder and the output .zip file containing the data that was clipped.  According to the messages.xml file in the scratch folder, the job finishes successfully.

I have checked the ArcGIS Server and Windows logs and haven't found anything obvious that would point to an issue.  Folder permissions should be correct, since the task can create the job folder, clip the data, and create the output .zip folder.  The task will run successfully when used through ArcGIS Desktop.

So has anyone experienced similar behavior with the data extraction tool or any other GP tasks?  Any ideas of where I can look next?

Environment:
ArcGIS Server 10, .NET, SP2 - Distributed installation
ArcGIS Desktop 10, SP2
Windows Server 2008 x64 R2
Thanks,
Tim

Update:
I did some testing today (in Firefox) and found a few more things:
 1.  The GP Service does occasionally work, it seems.  Every few tries I do get a prompt to download the extracted data, only to get a security sandbox violation:
http://SERVER IP/Samples/DataExtract/Flex/widgets/DataExtract/DataExtractWidget.swf cannot load data from http://SERVER NAME/arcgisjobs/sumner.dataextract_gpserver/j77643b1b114e458f82e6f416262aaf58/scratch/output.zip.
 2.  When I am prompted to download the file, Firebug shows that it cannot download my crossdomain.xml file.  It shows the address as http://SERVER NAME/crossdomain.xml

In both cases, it is trying to access a http://servername URL, and can't get there.  Problem is, I can't figure out where it's getting this URL.  I tried modifying the virtual directory URL's and republishing the services using the ArcGIS Server IP instead of machine name, with no change in behavior.


texts:
I would like to count the number of points that intersect with a list of polygons.

The polygons are very complex (county/state boundaries).

I'm using this query:

SELECT gid, county_name,
   (SELECT COUNT(*) FROM points WHERE
     st_intersects(ST_Transform(points.location::geometry, 4326),
     ST_Transform(counties.geom::geometry,4326))
   ) AS c FROM counties;

In a query where there are 630 counties and 100 points, the query takes 290,000ms.

The difficulty is that I have 100,000 points and this query will then take 4 days.

Needless to say, this isn't viable.

The outcome would be:

County    | count
------------------
Wiltshire | 356

Do you have any advice on how to dramatically improve the efficiency of this query?


texts:
I'm trying to publish Mapbox vector tiles from my GeoServer using OpenLayers. I'm a total novice when it comes to everything around web mapping, web design etc. I'm using GeoServer deployed locally on my Windows 10 machine through Tomcat. I´ve gone through this official tutorial, but pointing to my own data, code looking this this:
<!DOCTYPE html -->

  <title>Vector tiles</title>

  <link rel="stylesheet" href="ol.css">

Mapbox Protobuf - vector tiles

When I try to open the map, the data does not display, and inspecting the page shows the error message:

Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at file:///geoserver/gwc/service/tms/1.0.0/SKS:sksbiotopskydd@EPSG%3A3006@pbf/2/0/2.pbf. (Reason: CORS request not http).

I´ve tried to read up on CORS, but not sure I totally grasp the concept. I've also tried activating CORS in Tomcat (by uncommenting the filters in the web.xml), but that doesn´t help (gives the same error message). What do I need to do to get this working?


texts: Is there a setting in QGIS that can stop auto-docking of dragged windows?

texts: I apologize if this seems like a trivial question, but here goes: Often I find myself dragging windows to get them out of the way of the map, and QGIS decides to dock them in very inconvenient ways. Can be a real pain sometimes, as I have to repeatedly click to undock the window, and wait for QGIS to pop it back out. Is there a setting which can turn this behavior off?