texts: Get tile URL with projection different than 900913

texts:
I have the following code, which should display a slippy map showing tiles producted by me with mapnik:

Proj4js.defs["EPSG:31258"] = "+proj=tmerc +lat_0=0 +lon_0=13.33333333333333 +k=1 +x_0=450000 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs";
var latLonProj = new OpenLayers.Projection("EPSG:4326");
var salzburgProj = new OpenLayers.Projection("EPSG:31258");
var bounds = new OpenLayers.Bounds(426249,294999,426251,295001);
var maxZoom = 19;
var minZoom = 11;

var map1 = new OpenLayers.Map({
               div: "map1",
               allOverlays: true,
               maxExtent: bounds,
               maxResolution: "auto",
               projection: salzburgProj,
               displayProjection: salzburgProj,
               allowSelection: true,
               units: 'm',
               controls: [new OpenLayers.Control.PanZoomBar(), new OpenLayers.Control.Navigation()]
    });

var layer1 = new OpenLayers.Layer.XYZ(
                 "layer",
                 "tiles/",
                 {
                     sphericalMercator: true,
                     crossOriginKeyword: null,
                     buffer: 1,
                     getURL: get_my_url
                 });

map1.addLayer(layer1);
var center = bounds.getCenterLonLat();
map1.setCenter(center, minZoom);

function get_my_url (bounds) {
    var z = this.map.getZoom();
    var x = long2tile(bounds.transform(salzburgProj,latLonProj).left, z);
    var y = lat2tile(bounds.transform(salzburgProj,latLonProj).top, z);

    var path = z + "/" + x + "/" + y + ".png";
    var url = this.url;
    if (url instanceof Array) {
        url = this.selectUrl(path, url);
    }
    return url + path;
}

function long2tile(lon,zoom){
    return (Math.floor((lon+180)/360*Math.pow(2,zoom)));
}

function lat2tile(lat,zoom){
return (Math.floor((1-Math.log(Math.tan(lat*Math.PI/180) + 1/Math.cos(lat*Math.PI/180))/Math.PI)/2 *Math.pow(2,zoom)));
}

I have taken get_my_url() from here (and customized it) and long2tile and lat2tile from here.

The problem is that I do not get the correct tile urls, probably because long2tile() and lat2tile() refers to projection 900913 and not to my salzburg projection 31258.

How can I modify these functions in order to get correct tile urls? Maybe can I use some tuning parameter given in 31258 definition at the top of my source?

Thanks in advance.


texts:
I'm having trouble to load a GeoTIFF file using Geotools (never used it before).
What I've tried so far:
File koppenClassificationTif = new File("mongodb-repository/src/main/resources/land_shallow_topo_2048.tif");
GeoTiffFormat format = new GeoTiffFormat();
Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
GeoTiffReader reader = format.getReader(koppenClassificationTif, hints);

The GeoTiffFormat could not be initialized because of java.lang.NoSuchMethodError: 'javax.measure.Unit javax.measure.Unit.multiply(java.lang.Number)'. My .m2 repository contain 2 versions javax/measure/unit-api but I don't see why this would be a problem.
I've also tried to directly use GeoTiffReader reader = new GeoTiffReader(koppenClassificationTif); which throw this exception:

org.geotools.data.DataSourceException: class
org.apache.sis.referencing.factory.GeodeticObjectFactory cannot be
cast to class org.opengis.referencing.Factory
(org.apache.sis.referencing.factory.GeodeticObjectFactory and
org.opengis.referencing.Factory are in unnamed module of loader 'app')

I'm using newly released 25.0 but I've been experiencing the same issue using 24.2
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-geojsondatastore</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-referencing</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-coverage</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-opengis</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-geotiff</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-main</artifactId>
        <version>${geotools.version}</version>
    </dependency>
</dependencies>

<repositories>
    <repository>
        <id>osgeo-alt</id>
        <url>https://repo.osgeo.org/repository/release</url>
    </repository>
</repositories>

<properties>
    <geotools.version>25.0</geotools.version>
</properties>


texts:
I have problem with generating charts in GEE. I am trying to calculate NDVI for more than 100 Shapefiles (FeatureCollection) and export the results to a chart and CSV file. I got an error:

  Error generating chart: Collection query aborted after accumulating over 5000 elements

How can I revise my code to avoid this error?

// Define a geometry.
var Akdarya = ee.FeatureCollection("users/xodjayevshovkat/Akdarya_GEE_35");

// Make the geometry a feature collection.

// Import S2 image collection.
var S2 = ee.ImageCollection('COPERNICUS/S2')
  // Filter collection.
  .filterDate('2018-01-01','2018-12-31')
  .filterBounds(Akdarya)
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20));

// Define a cloud masking function.
var maskcloud1 = function(image) {
  var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1));
};

// Define an index function (return only NDVI).
var NDVI = function(image) {
  return image.expression(
    '(NIR - RED) / (NIR + RED)',
    {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
    }).rename('NDVI').copyProperties(image, image.propertyNames());
};

// Apply the cloud mask function.
var NDVIcol = S2.map(maskcloud1)
 // Apply the NDVI transformation.
 .map(NDVI);

// Calculate the per-pixel median NDRE for all images in collection.
var NDVImed = NDVIcol.median();

// Define vis params.
var NDVIvis = {
  palette: ['#d73027', '#f46d43', '#fdae61', '#fee08b', '#d9ef8b', '#a6d96a'],
  min: -0.5,
  max: 0.9};

// Plot a time series of
var plotNDVI = ui.Chart.image.seriesByRegion({
  imageCollection:NDVIcol,
  regions: Akdarya,
  reducer: ee.Reducer.mean(),
  //scale: 1000,
  xProperty: 'system:time_start',
  seriesProperty: 'system:index'})
  .setChartType('LineChart').setOptions({
    title: 'NDVI Akdarya 2018',
    hAxis: {title: 'Date'},
    vAxis: {title: 'NDVI'},
});

// Display the results.
print(plotNDVI);
Map.addLayer(NDVImed.clipToCollection(Akdarya), NDVIvis, 'NDVI');


texts: Creating network simulation / temporal models for videos?

texts:
I know of ArcGIS Tracking analyst and QGIS Time Manager. I have a transportation network with time based schedule for buses. I will create time aware shapefiles from this data. Is it possible to create a simulation that shows buses moving along a track (as per the time ) using some plugin? This is necessary as i want to convert it to a video.

Essentially can you all list all Add-ins / Plugins that deal with time based data in popular GISystems.


texts:
I need to setup a continental geoportal, that would be based on a series of geonetwork servers disseminated in the continent (about 6 servers).

In the geonetwork documentation, http://nationaalgeoregister.nl/geonetwork/docs/eng/users/admin/harvesting/index.html
They give an example of a hierarchical harvesting: C harvest B, B harvests A, D harvest C, B and A. According to the documentation, there is no problem due to the UUID.

In my case, the harvesting would be circular: say we have S1, S2, S3, S4, S5, S6,
then S1 harvests S2, S3, S4, S5, S6; S2 harvests S1, S3, S4, S5, S6, etc.
So the user can interrogate the closest geoserver, and get information for the whole continent.

Will this kind of circular approach work, or shall I fall back to a hierarchy?


texts:
I've a polygon layer for parcels in a landscape protection area and a point layer, contains zoological data (e.g species, entities, etc.).

How can I make a statics for each polygon, contains the list of observed species, the count of observations by species and the sum of entities by species?

I would like to make it in QGIS with shapefiles or spatialite database.


texts: Converting CSV file of WKT data to shapefile using QGIS 2.01

texts:
Is it possible to create a shapefile directly from a CSV WKT file in QGIS 2.01?
WKT
POLYGON(...)
POLYGON(...)

Save as shapefile doesn't seem to work.


texts: Installing rgdal - cannot solve dependencies

texts:
I am trying to build a docker image based on osgeo/gdal:ubuntu-XXXXX. However it keeps failing when trying to install rgdal. I've tried this on gdal-3.1.1/proj-7.1.0 and gdal-3.0.3/proj-6.3.0, as well as from ubuntu-full-v2.5.0RC1.
After some deep digging, it looks like the error stems from this:
configure:2847: g++ -std=gnu++11   -DPROJ_H_API -I/usr/local/include  -Wl,-Bsymbolic-functions -Wl,-z,relro conftest.cpp -L/usr/local/lib -lproj    -L/usr/lib -lgdal >&5
/usr/bin/ld: cannot find -lproj
collect2: error: ld returned 1 exit status
configure:2851: $? = 1
configure:2889: result: no
configure: failed program was:
| /* confdefs.h */

resulting in a status that looks like this:
...
configure: rgdal: 1.5-12
checking for /usr/bin/svnversion... no
configure: svn revision: 1018
checking for gdal-config... /usr/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 3.0.3
checking GDAL version >= 1.11.4... yes
checking GDAL version <= 2.5 or >= 3.0... yes
checking GDAL: linking with --libs only... yes
checking GDAL: gdal-config data directory readable... yes
checking GDAL: /usr/share/gdal/stateplane.csv readable... yes
configure: pkg-config proj exists, will use it
configure: PROJ version: 6.3.0
configure: PROJ CPP flags: -DPROJ_H_API -I/usr/local/include
configure: PROJ LIBS: -L/usr/local/lib -lproj
checking PROJ header API:... yes
checking whether the C++ compiler works... no
configure: error: in `/tmp/RtmpxFc0qV/R.INSTALL1e86d5b26dc/rgdal':
configure: error: C++ compiler cannot create executables
See `config.log' for more details
ERROR: configuration failed for package 'rgdal'

Some googling led me to this stackoverflow.
I tried to munge things by ln -sv /usr/local/lib/libproj.so.19 /usr/local/lib/libproj.so but that just results in:
configure: error: proj_context_create not found in libproj.

I then tried some
export PKG_CONFIG_PATH=/usr/local/lib/pkgconfig
export LD_LIBRARY_PATH=/usr/local/lib
autoconf
./configure
Still complaining of /usr/bin/ld: cannot find -lproj !
Here's a Dockerfile you can use to reproduce:
ARG BASE_IMAGE=osgeo/gdal:ubuntu-small-latest
FROM ${BASE_IMAGE} as base

# Setup tzdata and set timezone so we don't get prompted for it later
RUN     echo 'Etc/UTC' > /etc/timezone \
    &&  export DEBIAN_FRONTEND=noninteractive \
    &&  ln -sf /usr/share/zoneinfo/America/New_York /etc/localtime \
    &&  apt-get update -qq \
    &&  apt-get install -q -y tzdata \
    &&  rm -rf /var/lib/apt/lists/* \
    &&  dpkg-reconfigure --frontend noninteractive tzdata

RUN     apt-get update -qq \
    &&  apt-get install -qq autoconf curl git build-essential dirmngr software-properties-common \
    &&  rm -rf /var/lib/apt/lists/*

RUN     apt-key adv --no-tty --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9 \
    &&  add-apt-repository 'deb https://cloud.r-project.org/bin/linux/ubuntu focal-cran40/' \
    &&  apt-get update -qq \
    &&  apt-get install -y r-base \
    &&  rm -rf /var/lib/apt/lists/*

RUN echo 'install.packages("rgdal")' | R --no-save


texts: Manually moving multipart label results in disappearing labels

texts:
When I attempt to move a multipart label, all of the labels except for the one I am moving disappear. For example, when I move the lower "2-C" label (Figure 1) the upper two "2-C" labels disappear from the display and map layout (Figure 2).

I have encountered a related, but different question that addresses single labels disappearing in QGIS 2, though not multipart labels in QGIS 3.

How can I manually move a label and preserve all of the corresponding multipart labels using QGIS 3.10?

Figure 1

Figure 2


texts: QGIS Transforming from ETRS89 to WGS84 - why are Helmert parameters 0?

texts:
If ETRS89 (EPSG:4258) is tectonically fixed and WGS84 (EPSG:4326) isn't, the Helmert coefficients shouldn't be zero, right? In fact these parameters should diverge more and more with every new release of QGIS

I understand they're based on the same ellipsoid, but I think their central meridians are diverging, no? Why is the transformation (0,0,0,0,0,0,0)??


texts: Finding Overlapping/Close Points of Same Layer in ArcGIS Pro

texts:
I have a dataset of point features with photo attachments of street signs from volunteers through the Survey123 app. The app records their location which ultimately creates the point feature.
Part of the project involves identifying and removing photo submissions of the same signs, which I think can be accomplished by identifying the points that overlap each other or are near each other - this would indicate that a photo was taken at the same location which means the chance of the same sign being photographed is higher and could be manually reviewed.
The select by location tool doesn't seem to work properly as any configuration just selects all the records.
Is there a tool that could accomplish what I'm trying to do?


texts: How to find income data by census tract or block?

texts:
I am trying to find a file with income data for every census tract.

On the American FactFinder/Census websites, I can find data for each individual state or county but I want all census tracts in all 50 states.

Does this mean I have to download all 50 files and combine them?


texts: Exporting an ImageCollection to a CSV in Google Earth Engine

texts:
I am pretty new to using the code editor. I am trying to export ten-day NDVI composites (images from Sentinel-2) to a table (as a CSV file). I think my problem lies with not knowing how to convert my ImageCollection to FeatureCollection (btw. in this version survey data is simplified).
var surveydata = ee.Geometry.Point (16.13298, 77.85595);

// study area
var studyArea = ee.FeatureCollection(surveydata);
var surveyData = ee.FeatureCollection(surveydata);

// Load Sentinel-2 data.
var s2 = ee.ImageCollection('COPERNICUS/S2');

// Load cloud masking tool
var s2mask = require('users/fitoprincipe/geetools:cloud_masks').sentinel2;

// function to apply NDVI calculation to every image in collection
function addNDVI(image) {
  var NDVI = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(NDVI);
}

//set a list of start dates
var startDate = '2020-05-01';
var endDate = '2020-09-30';
var interval = 10;
var increment = 'day';
var tendayDifference = ee.Date(startDate).advance(interval, increment).millis().subtract(ee.Date(startDate).millis());
var listMap = ee.List.sequence(ee.Date(startDate).millis(), ee.Date(endDate).millis(), tendayDifference);

//function to get a tenday composite
function gettendaySentinelComposite(date) {
  // Only include the NDVI
  var s2 = ee.ImageCollection('COPERNICUS/S2')
                      .filterBounds(studyArea)
                      .filterDate(date, date.advance(interval, increment))
                      .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 60)
                      .map(s2mask())
                      .map(addNDVI)
                      .select('NDVI');
  var composite = s2.max()
                      .set('system:time_start', date.millis(), 'dateYMD', date.format('YYYY-MM-dd'), 'numbImages', s2.size());
  return composite;
}

var ndvi_tenday = ee.ImageCollection.fromImages(listMap.map(function(dateMillis){
  var date = ee.Date(dateMillis);
  return gettendaySentinelComposite(date);
}));

print(ndvi_tenday);

var tempTimeSeries = ui.Chart.image.series({
      imageCollection: ndvi_tenday,
      region: studyArea,
      reducer: ee.Reducer.max(),
      scale: 100,
      xProperty: 'system:time_start'
    })
    .setChartType('ScatterChart')
    .setOptions({
    title: 'NDVI-over-time',
    vAxis: {
      title: 'NDVI'
    },
    lineWidth: 1,
    pointSize: 4,
    series: {
      0: {color: 'blue'},
      1: {color: 'red'},
      2: {color: 'green'}
    }
});

print(tempTimeSeries);

// Convert the object to an image with one band for each tenday period
var ndvi_tenday = ee.ImageCollection(ndvi_tenday)
.toBands();

// code to plot one band (one tenday period)
var image = ndvi_tenday.select("9_NDVI");

Map.centerObject(studyArea, 9);
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
Map.addLayer(image, ndviParams, 'NDVI image');

print(ndvi_tenday);

I thought to use this code for exporting:
// Export the image to your google drive
Export.table.toDrive({
        collection: ndvi_tenday,
        description: 'NDVI',
        folder: 'EE',
    fileNamePrefix: 'NDVI',
    fileFormat: 'CSV'
    });

Here is the link for my code for clarification:
https://code.earthengine.google.com/?noload=1&scriptPath=users%2Fanjadamlisby%2FTest%3AS2_NDVI_GPS_Table


texts: Find and replace unicode character in a shapefile

texts:
I've found several example of script and blog posts on working with unicode characters, but I haven't been able to make anything works so far... It's a little frustrating

I have a shapefile encoded in utf-8 (exported from QGIS) and there's many É, È, À, Ô, etc. in the values of some fields.

For a geocoder, I need to normalize my fields values.

This is what I have so far:

# -*- coding: utf-8 -*-
###############################################################################
import sys
from sys import argv
import osgeo.ogr
from amtpy import EndOfScript
###############################################################################

script, src_shp, fld_shp = argv

#- Opening the shapefile
shapefile = osgeo.ogr.Open(src_shp)
layer = shapefile.GetLayer(0)
spatialRef = layer.GetSpatialRef()

#- Going through each feature, one by one
for i in range(layer.GetFeatureCount()):
    print "Normalisation de la ligne %i" %(i+1)
    feat = layer.GetFeature(i)

    texte_norm = feat.GetField(fld_shp)
    texte_norm = texte_norm.encode('utf-8')

    texte_norm = texte_norm.upper()
    texte_norm = texte_norm.replace(u'\00c2', 'A') #À
    texte_norm = texte_norm.replace(u'\u00C9', 'E') #É
    texte_norm = texte_norm.replace(u'\u00C8', 'E') #È
    # I've remove 16-17 characters to replace, for the example...
    texte_norm = texte_norm.decode('utf-8')

    print texte_norm
    #feat.SetField(fld_shp, texte_norm)

#- print EndOfScript message
print EndOfScript()

I've put a print instead of the SetField just so I could see if the replace works.

Then the script encounters the first unicode character, I have this error message:
'ascii' codec can't decode byte 0xc3 in position 5: ordinal not in range(128).

Well, I know that, since the character is a "À". I can't seem to find the right way.

Any tips?

Thanks!


texts: Preserving origin IDs when creating Origin-Destination (OD) Cost Matrix?

texts:
I have created an Origin-Destination cost matrix showing the distance in feet between several thousand home addresses (origins) and about 24 rental vehicles (destinations) at various points throughout the city. Each address point that was input into the OD matrix analysis also originally contained information on the number of times the resident at each address had reserved a particular vehicle.

My goal is to create a table that will show, for all home addresses, the distance from that address to each of the 24 rental cars and the number of times the resident reserved each of the 24 cars. I could then say, for example, "40% of reservations made by residents of Neighborhood A are for cars within 1/2 mile" or "Residents who reserved car X lived an average of 2000 feet away".

Here is the problem. In conducting the OD matrix analysis, it seems that all the unique identifiers of the data that I input have been obliterated. It would appear that there is no way to associate the distances between origins and destinations back to the other attributes associated with the home address points except for some sort of spatial join. Yet, between the inaccuracies introduced by geocoding thousands of addresses and the fact that many home locations are apartments with multiple residents in the same building, it is impossible to distinguish which resident made which trips by using a spatial join.

How have others typically handled the origin-destination problem when origin points overlap and still need to be associated back to other attributes?


texts:
In ArcGIS I have used to share a shapefile (for example beach areas) by sending a package that contains different files (.shp, .shx, .prj, .cpg, .dbf, .sbn, .xml) which contains all necessary properties for the shapefile. So far everyone has been happy.

But now I got NetCDF file of average wind speed. I successfully imported it from ArcToolbox -> Multidimensional Tools -> Make NetCDF Raster Layer. From layer properties I changed the "Color ramp" to be "red to blue", so that blue color is 3 m/s and red is 8 m/s.

Now, what is the proper way to export/save that wind speed raster layer? So that I could send it by e-mail and everyone could see the same colorful picture and the "color ramp." Is GeoTIFF-format mostly used? Or is there something better for sharing raster layer?


texts: What are considered invalid fieldnames with arcpy Cursor classes?

texts:
Using this code, I get the following error.
infeature = r"C:\pathto\file.gdb\file" #location of feature class
fields = ["Name:", "Location", "Email Address", "Cell Phone Number:",
          "Home Phone Number:"] #these are the fieldnames, not the aliases

with arcpy.da.SearchCursor(infeature, fields) as cursor1:
    for row in cursor1:
        print(row)

Traceback (most recent call last):
  File "pathto\DuplicateFields_updated.py", line 18, in <module>
    for row in cursor1:
RuntimeError: An invalid SQL statement was used. [SELECT OBJECTID,Name:,Location,
Email Address,Cell Phone Number:,Home Phone Number: FROM backup_02232021_NEW]

In the help docs and all the cursor-related questions here, I haven't found anything telling me how the fieldnames in a table should look. I have spaces and colons in mine. But, what if the fieldnames cannot be changed-- as in my case? I'm using a SearchCursor here, but my question goes for the other cursors too.


texts: Running QChainage from python script in QGIS?

texts:
I'm fairly experienced in both python and QGIS but not so much together.

I have a series of vector lines I need to run the very useful QChainage plug-in on to get points every 1000m. I have about 200 of these .shp files and have already written a rough script to reproject and dissolve them using the processing module, but I am not sure how to do that in the script. I am guessing the code would look something like this:

Qchainage.run(input='input.shp', distance=1000, output='output.shp')

But obviously that doesn't seem to be working. I cna't find any documentation on running qchainage from python either.


texts: Identifying feature OID from Geometry object using ArcPy?

texts:
I know ArcPy allows you to access geometries directly using the CopyFeatures_management method. Note in their code sample for the Geometry object, it is very simple to copy the geometries to in-memory objects of type Geometry (and extremely fast). The geometry object of course relates to the geometry/shape and completely bypasses attribute information.

My question is whether there is a way to identify the OID of each geometry cheaply?

A little background in what I am trying to perform:

For all feature classes in a geodatabase
Read through each feature's geometry
Go over the points in the geometry
Calculate the distances between the points
If the distance is below a threshold, that specific feature needs to marked for review (meaning I need either the OID/FID or any other field to identify it).

As I am still in ArcGIS 10.0, I would prefer to not use cursors that seem to work a lot slower than simply dealing with the geometries.


texts: Type mismatch error in Calculate Field tool in ArcGIS 10.1

texts:
I'm trying to run a rather lengthy model (built in ArcGIS 9.3.1) in ArcGIS 10.1 and I'm recieving an error on a Calculate Field tool. The error is - Type mismatch: 'iif'.

This error is refering to an expression in the tool: iif(IsNull( [PID] ) ,NA,[PID]).

I tried "if" instead of "iif" but then recieved a generic error.

I never had this issue when running the model in 9.3.1. I'm not sure what to do here.  Any suggestions?


texts:
I'm attempting to iterate through several shape files within a folder and only move some of them to a new folder.  I'm doing this based on contents of the DBF file where I'm using dbfread to determine whether is should be moved or not.  I have this all working as expected except it's only moving the DBF file and not the .shp, .shx, and .prj. files since the logic is only on the DBF.  One way I was thinking I could move all files would be to add the basename of the ones I want to move to a List but I seem to be failing at that.  Here is the code I'm using.

from dbfread import DBF
MoveFolder_path = r"C:\Data\Move"
move_list = []
for file in directory_list:
    if file.split('.')[1] == 'dbf':
        file_basename = file.split('.')[0]
        table = DBF(file, load=True)
        if table.records[0]['name'] == None:
            move_list.append(file_basename)
for movefile in move_list:
    for file in directory_list:
        file_basename = file.split('.')[0]
        if file_basename == movefile:
            shutil.move(file, MoveFolder_path)


texts: QGIS connection to ESRI feature server: Failed to retrieve service capabilities: Error 499: Token Required

texts:
I've created a survey in Survey123 Connect, published it to ArcGIS online and am attempting to connect to the feature server to view the collected survey data in QGIS 3.10.2. To add the feature service I'm using Layer, Add Layer, Add ArcGIS FeatureService Layer.

The feature service link is posted below the following error:

Failed to retrieve service capabilities: Error 499: Token Required

https://services8.arcgis.com/8tX9OH8doNBidChm/ArcGIS/rest/services/service_edbe38394aa44de28045d39f27ffaecc/FeatureServer/0?token=YyrjuSR45GWMoH0A0WIifewkW_LAiYvu1OVr1HXJLNJSRO9Lphb48Jr8fxjofTEeRf3UPMoIIHJVIe-DSpp6vuNVlaB_GzvfykiA8Mm9autHVoWvZjWKI1W0uKhcV2lwiSoPQCKG42KtKKz0lc2ZyjAxVti8WJPhYvhKxAR-wrgOPOoXc17xEwfXo8bPN3nijUVnP83459I3D9GKSe79j3L3gYwPRHM_ONO-BtJCCQlAPwiqS9R1OOKeclbJzcMV

To establish the connection I've only used the following portion of the URL: https://services8.arcgis.com/8tX9OH8doNBidChm/ArcGIS/rest/services/service_edbe38394aa44de28045d39f27ffaecc/FeatureServer/

I've added the same feature layer in ArcMap 10.5 with no issues (using add data from ArcGIS online).

I've tested another survey feature service layer that I built in Survey123 Connect under a different ESRI username and organization and it loads with no issues in the same QGIS project: https://services2.arcgis.com/NlsizNmbMFiinWw4/ArcGIS/rest/services/service_0ca13752c9fb461f9b1d1b2490f79504/FeatureServer/

If you have any ideas I'd love to hear them.


texts: Possible to set default attribute fields for vector creation?

texts: I'm using the latest release of QGIS (2.14). Is it possible to set default attribute field parameters to be used automatically during the creation of new shapefiles/vector files? I'd be specifically interested in a way to do this on the project level, but if there is a way to specify this on the system level that would also be relevant to my needs.


texts: Attribute data is not supported by the GeometryProvider in th Sharpmap

texts:
Attribute data is not supported by the GeometryProvider.

SharpMap.UI ver 1.2.0.0

I want to make a click event on the image("maker"). so I used the following sample code URL.
https://xaml.dev/post/2006/05/18/Using-NetTopologySuite-in-SharpMap

but I face the exception.(SharpMap/Data/Providers/GeometryProvider.cs line 211)
What do I make a source code in function("OnExecuteIntersectionQuery", "ExecuteIntersectionQuery")?

  //Convert click-point in image coordinates to world coordinates
    Envelope ClickPnt = new Envelope(mapBox.Map.ImageToWorld(new System.Drawing.Point(e.X, e.Y)));
    SharpMap.Data.FeatureDataSet ds = new SharpMap.Data.FeatureDataSet();
    //Execute click-query on first layer in layers collection
    (mapBox.Map.Layers[1] as SharpMap.Layers.VectorLayer).DataSource.ExecuteIntersectionQuery(ClickPnt, ds);
    if (ds.Tables.Count > 0) //We have a result
    {
        //Add clicked features to a new selection layer
    }

public override void ExecuteIntersectionQuery(Envelope box, FeatureDataSet ds)
{
    throw new NotSupportedException("Attribute data is not supported by the GeometryProvider.");
}


texts: Determining raster intersection using ArcGIS for Desktop?

texts:
Using Arc 10.2.2 with Spatial Analyst, I need to determine the spatial intersection of two rasters (both rasters are in the same coordinate system, and are aligned to the same snap raster).  The two rasters have different values.

Specifically, I want to determine which cells are occupied by both raster A and raster B.

Here's my method, using the Raster Calculator:

I create an intermediate raster each for raster A and B, calculating the output to the value 1.  This identifies, for each raster, which cells are occupied:

intermediate_A = (A * 0) + 1

intermediate_B = (B * 0) + 1

Add the two intermediate rasters together:

final_output = intermediate_A + intermediate_B

The final_output raster will have the value 2 for the intersecting cells, otherwise NoData.

Everything seems to work correctly.  However, I'm wondering if there's a more efficient way to determine the spatial intersection with fewer steps or without creating the intermediate rasters?


texts:
Reading a DGN file is a big problem for me. I tried more and more but can't find a way to read a dgn(v8) file by VB.Net. After many search, i found out that GDAL can read and write dgn file, but i 'm not able to use it.
Can anyone help me to do that ?


texts:
How can we integrate ArcGIS maps with RaVis BirdEye's VisualGraph?
My intention is to place the RaVis VisualGraph over the map and the nodes of the visual graph  must be placed on their specific location on the map...

Can anybody help me in this regard?


texts: QGIS creating a DEM using OS Terrain 50 data - component tiles not matching

texts:
I'm attempting to create a DEM from OS Terrain 50 data (ASCII Grid data) in QGIS.

The data is provided as multiple zipped folders representing sub divisions of OS grid squares. Each zipped folder contains an .asc file, a .gml file, and two .xml metadata files.

Dragging and dropping an .asc file into a layer gets me a DEM raster. The problem is that when I add adjacent .asc files, the created rasters don't blend. The edges have different brightness levels giving a tiling effect.

What am I doing wrong?


texts: Indexing to Open Data Cube not working for land cover map

texts:
I'm trying to index a clipped version of ESA's 2015 land cover map into an Open Data Cube instance. I've written a product specification yaml and a yaml for te actual product itself and it indexes into the data cube fine. However, when I try to load the data later on using dc.load, it gives the following error:

~/miniconda3/envs/cubeenv/lib/python3.6/site-packages/datacube/api/core.py in <listcomp>(.0)
643
644 def get_bounds(datasets, crs):
--> 645     left = min([d.extent.to_crs(crs).boundingbox.left for d in datasets])
646     right = max([d.extent.to_crs(crs).boundingbox.right for d in datasets])
647     top = max([d.extent.to_crs(crs).boundingbox.top for d in datasets])
AttributeError: 'NoneType' object has no attribute 'to_crs'

This makes me thinks that it is an error with the yaml coordinate specifications as it isn't picking up any of the actual data at all when trying to define the bbox coordinates of the data. I can only think that this is because I'm trying to index data in EPSG:4326 coordinate system, as the yaml is basically identical to previous yamls which have worked, only those yamls involved indexing data in UTM coordinate systems.

Both the product specification and individual product yamls are attached.

Product specification link: https://slack-files.com/T0L4V0TFT-FGXRVUWUV-b4eea86062

Product yaml link: https://slack-files.com/T0L4V0TFT-FGYJ8N223-ebe40ec2b0


texts:
I want to use mdenoise / r.denoise but I can't get it working. I have Xubuntu 12.04 and Grass 6.4.1

GRASS tells me after installing the plugin:

g.extension extension=r.denoise svnurl=https://svn.osgeo.org/grass/grass-addons/grass6
Fetching <r.denoise> from GRASS-Addons SVN (be patient)...
A    r.denoise/description.html
A    r.denoise/r.denoise
A    r.denoise/Makefile
 U   r.denoise
Ausgecheckt, Revision 51599.
Compiling <r.denoise>...
Makefile:5: /usr/lib/grass64/include/Make/Script.make: Datei
oder Verzeichnis nicht gefunden
make: *** Keine Regel, um
»/usr/lib/grass64/include/Make/Script.make« zu erstellen.
Schluss.
ERROR: Compilation failed, sorry. Please check above error messages.

I got no solution for that so I tried to use mdenoise in the terminal but here I got trouble with the file extensions: ASC-is the only extension I can work with but I got a "memory access error". Converting to other extensions failed (using SAGA).

root@MAP-PC:~/mdenoise# /home/map/mdenoise/mdenoise -I esri.asc -o output.asc
Input File: esri.asc
Neighbourhood: Common Vertex
Threshold: 0.400000
n1: 20
n2: 50
Speicherzugriffsfehler (Speicherabzug geschrieben)
root@MAP-PC:~/mdenoise#

with xyz-extension:

root@MAP-PC:~/mdenoise# /home/map/mdenoise/mdenoise -I esri.xyz -o output.asc
Input File: esri.xyz
Neighbourhood: Common Vertex
Threshold: 0.400000
n1: 20
n2: 50
Read Model...
Triangulation...
    13.620 seconds
Denoising Model...     0.000 seconds
Saving Model...No .prj file is found.
     0.120 seconds
root@MAP-PC:~/mdenoise#

I'm working now 10 hours non-stop on it! I'm going crazy!

I really hope, someone can help me, to get mdenoise running!


texts: OpenStreetmap Multipolygon with Members without role

texts:
I found this relation in OpenStreetMap

http://www.openstreetmap.org/relation/1561932

http://www.openstreetmap.org/api/0.6/relation/1561932/full

It has only one outer member but it needs the other members to be a valid multipolygon. Is it common that Members without role are treated as outer?


texts: What affects the kind of date (local vs. UTC vs unspecified) on JTXJobClass.StartedDate?

texts:
I am accessing the JTXJobClass.StartedDate for a JTX Job (with IJTXJob3).  I am getting an unspecified DateTime.Kind value.  The value appears to be UTC (iirc, WMX always stores UTC in the database but I'm questioning that since the DateTime.Kind is unspecified). The value comes from the started_date field of the JTX_JOBS table and its values looks like UTC to me.

Is there a WMX configuration that would make the DateTime.Kind on the JTXJobClass.StartedDate return as Utc? Or Local?

Or is it safe to assume that the unspecified value is Utc?  Instinctively, this is my first approach (because it is easy), but then I started to think about what could go wrong if that assumption doesn't hold true.


texts: Transforming aspect values in order to use them in multiple regression models

texts: I want to use aspect as an independent variable in a multiple regression model but I would like to transform the aspect values so that from degrees 180 to 360, I would get their equivalent to degrees from 0 to 180. How could I do that in ArcGIS or R?


texts: Bulk importing raster data into GeoServer

texts:
I have more than 1000 world images in tiff format. I need to create the raster layers for those images in GeoServer 2.4-RC2 on Windows XP OS and set it as a backdrop in my web application.
Is there any way for bulk importing in GeoServer? or
any way to group those images into one raster store? or
Is there any opensource to do this job?

Please suggest me.thanks in advance.


texts:
Licence: Advanced
  Version: ArcMap 10.3

In my shape are thousands of overlapping polygons. My goal is, to generate a new output in which the biggest overlapping value represents a new polygon/raster cell.

I've tried the Conversion Tool "Polygon to Raster" and set diffrent "cell assignment type" configurations. It works fine, if only 2 polygons overlapping each other. In the output raster represents each pixel value the biggest value of the polygon. However, in cases with >2 overlapping polygons it doesnt work for me.

Maybe there is a solution with QGIS too? The output should be a raster or polygon shape.

  Input: polygon shape
  Output: Raster or polygons


texts: Downloading files with file paths informed by attribute information

texts:
I am looking to download files that are saved in a server location with a file path that is informed by attribute information created in the graphical modeler.

The information that I would like to add to the file path is stored as values in each row of the attribute table under the field 'TILENAME'.

I would like to input these values (e.g. 'SK74NE' in the field 'TILENAME') to inform the file path:

'file:///X:/GIS/Location Map/OS Open Map Local Apr 2019/SK/SK74NE.tif'

I have tried various appropriates using select by attribute, extract by attribute and refactoring fields but haven't had any luck.


texts: I'm doing a project in PostGIS where I need to know what is the conversion of meters to radians. I would like to know how to do or what to look for to carry out the doubt.


texts:
When adding a layer to a custom style and symbolizing the layer, the new layer doesn't seem to show up in the published style. We have waited a day to ensure the new tiles have had time to render.

For example, here is a layer in Studio:

And in a leaflet map:

Notice that the black dashed line, which is the layer that was added to the style is not sowing up. We've cleared the cache in the browser. Any suggestions on how to proceed?


texts: How to fix QGIS Desktop startup crashes with Plugins already disabled?

texts:
Yesterday everything was working fine, but since today I can't start Qgis anymore, because it crashes.

This affects all the versions that I have currently installed (2.8, 2.10. 2.12 and 2.14). The program loads and when it says "Starting Python", it crashes and writes a minidump (see screenshot below). QGis Browser launches normally.

What did I do yesterday: I used Grass (v.clean) and changed/cleared the folder paths for the plugin as explained in the answer to this question. Also I tried out something with a function in field calculator. But can all this cause all installed versions to crash at startup?

What could be causing this? Is there a way to manually fix this? I also tried a restart, but this didn't work. Edit: I use Windows 7 64-bit.


texts: Convert planar distance to ellipsoid (real) distance

texts:
I have the following linestring which approximately marks an airport runway:

LINESTRING(-381578.819887 6574044.360475, -378513.534443 6574799.475787)  # EPSG 3857
LINESTRING(-3.427781 50.732117, -3.400245 50.73641)                       # EPSG 4326

If I use the GLength function, I get a length of

SELECT GLength(GeomFromText('LINESTRING(-381578.819887 6574044.360475, -378513.534443 6574799.475787)'))
3156.925m   # EPSG 3857
SELECT GLength(GeomFromText('LINESTRING(-3.427781 50.732117, -3.400245 50.73641)'))
0.027869°   # EPSG 4326 ~ 0.027869 x 111319.9 = 3102m

The actual length of this runway is 2076m. QGIS will display this correct distance if I check the ellipsoidal box using the measure tool, so I assume this performs a projection/trigonometry conversion. Yet, I can't seem to find the correct calculation to replicate this conversion and go from the (planar?) distance returned by the GLength function to the real distance.

Can anyone direct me to a tutorial for this?

This question is a continuation of my pervious question How can I get the ellipsoidal values for length and area using spatialite? which got side tracked into correcting a different issue.


texts: Create IPoint with Lat/Long in decimal degrees and Z in feet

texts:
I have written an ArcMap v10 AddIn that reads a folder of geotagged JPEG's and creates a new feature class where each point is the geotag information from on of the JPEG's.  The issue is:

The GPS information from the exif information in each JPEG is in Latitude/Longitude with Degrees, Minutes, Seconds and I have successfully converted to decimal degrees.  The Altitude information in the exif properties is in feet.

I recently discovered that ArcMap is apparently interpreting the Z value I provide when setting point.Z = zAltitude as being in the decimal degree unit instead of feet.  This happens despite the fact that I have set the spatial reference information of the feature class and the IPoint to a geographic coordinate system of WGS84 and a Vertical Coordinate System of NAVD88 and a ZCoordinateUnit of feet.

Below is the code that creates the spatial reference and the code that creates the points:

Private _spatialRef As ISpatialReference3 = CreateSpatialRef()

Private Function CreateSpatialRef() As ISpatialReference3

    Dim COORDSYS As ESRI.ArcGIS.Geometry.esriSRGeoCSType = ESRI.ArcGIS.Geometry.esriSRGeoCSType.esriSRGeoCS_WGS1984
    Dim VERTCOORDSYS As ESRI.ArcGIS.Geometry.esriSRVerticalCSType = ESRI.ArcGIS.Geometry.esriSRVerticalCSType.esriSRVertCS_NAVD1988

    Dim spatialRefFact As ISpatialReferenceFactory3 = New SpatialReferenceEnvironment
    Dim geogCS As IGeographicCoordinateSystem2 = spatialRefFact.CreateGeographicCoordinateSystem(CInt(COORDSYS))
    Dim vertCS As IVerticalCoordinateSystem = spatialRefFact.CreateVerticalCoordinateSystem(CInt(VERTCOORDSYS))

    Dim spatialRef As ISpatialReference3 = CType(geogCS, ISpatialReference3)
    spatialRef.VerticalCoordinateSystem = vertCS

    Dim spatialRefResolution As ISpatialReferenceResolution = CType(geogCS, ISpatialReferenceResolution)
    spatialRefResolution.ConstructFromHorizon()
    spatialRefResolution.SetDefaultXYResolution()
    spatialRefResolution.SetDefaultZResolution()

    spatialRef.ZCoordinateUnit = CType(spatialRefFact.CreateUnit(CInt(esriSRUnitType.esriSRUnit_Foot)), LinearUnit)

    Return spatialRef

End Function

Private Function CreatePoint(ByVal xLong As Double, ByVal yLat As Double, ByVal zAlt As Double) As IPoint

    Dim pnt As IPoint = New ESRI.ArcGIS.Geometry.Point
    pnt.SpatialReference = _spatialRef
    pnt.PutCoords(xLong, yLat)

    Dim zAwre As IZAware = CType(pnt, IZAware)
    zAwre.ZAware = True

    pnt.Z = zAlt

    xLong = Nothing
    yLat = Nothing
    zAlt = Nothing
    zAwre = Nothing

    Return pnt

End Function

Any ideas how to do this without converting my Altitude from feet to decimal degrees?

Thanks in advance for your help.


texts: Converting GeoTIFF from EPSG:4326 to Spherical Mercator using QGIS

texts:
I'm using QGIS 3.12. I'd like to reproject NASA's Black Marble GeoTIFF (3600x1800 (0.1 degrees)) from https://earthobservatory.nasa.gov/features/NightLights/page3.php. The original projection is EPSG:4326. I would like to reproject it to EPSG:3857. Ultimately I'm going to take this reprojected GeoTiff into MapTiler to create basemap tiles.

I opened a new project in QGIS and dragged the GeoTiff onto a layer. I went to Raster > Projection > Warp, which pulled up this dialogue window:

There is no dropdown option for EPSG:3857 as a Target CRS. How might I set this target projection?


texts:
So, I've basically got what amounts to a vector-based raster set.  It's a series of regular-sized squares that, for all intents and purposes, are raster pixels. The only reason it isn't a raster is because I need to store multiple attributes for each polygon.
I want to aggregate these polygons into larger, hexagonal tiles.  My plan right now is to just do a merge and use the "has center in" option to aggregate the centers of the rectangles into whatever hexagon their center exists in.
Are there other, more accurate ways to do this?  Such as splitting the rectangles along the hexagonal edges and averaging out the attribute values based on the percentage split?


texts:
I am completely new to QGIS. I am on a Mac using OS 10.13.6. I have installed QGIS 3.6, and installed Python 3.6.8 per install instructions. I opened a terminal window to check python install:

EricsMacBookPro:~ ericfoertsch$ python
Python 3.6.8 (v3.6.8:3c6b436a57, Dec 24 2018, 02:04:31)

So I think I am good and on with the show.

I have 10 files I am trying to merge using Raster/Merge tool.

When I enter the file names in Raster/Merge dialog box I get this error:

  'Incorrect parameter value for INPUT' .

If I click on Run as Batch Process and enter same values I get this error message:

  'Wrong or missing parameter value:Input values(row 1)'.

I put double quotes around the file name strings in both input layers and merged file name and clicked Run again. This time I got logged output but with an error:

  GDAL command:
  gdal_merge.py -ot Float32 -of -o "\"/Users/ericfoertsch/Google
  Drive/Documents/Condor Scenary V2/New England/Merged Bil file/NewEngland.bil\"" --optfile
  /var/folders/sy/mbkmjwg12lb1df08dgjw7psh0000gp/T/processing_7e37da37ef8f4cfd89156358b42637eb/fa859f8db49c4623b356e04a58da9fdb/mergeInputFiles.txt
  GDAL command output:
  /bin/sh: gdal_merge.py: command not found

I opened a terminal window and typed in gdal_merge.py and it executed so I am not sure why it shows gdal_merge.py not found:

  EricsMacBookPro:~ ericfoertsch$ gdal_merge.py
  No input files selected.
  Usage: gdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*
                   [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]
                   [-ul_lr ulx uly lrx lry] [-init "value [value...]"]
                   [-n nodata_value] [-a_nodata output_nodata_value]
                   [-ot datatype] [-createonly] input_files
                   [--help-general]

Continuing the saga, I cut pasted the batch command into the terminal window:

  gdal_merge.py -ot Float32 -of EHdr -o "/Users/ericfoertsch/Google Drive/Documents/Condor Scenary V2/New England/Merged Bil file/NewEngland.bil" --optfile /var/folders/sy/mbkmjwg12lb1df08dgjw7psh0000gp/T/processing_7e37da37ef8f4cfd89156358b42637eb/58cfba9edf9f435cb75ca505c11d85ae/mergeInputFiles.txt

And got an error msg:

  No input files selected.

So I checked the ..../mergeInputFiles.txt file and it was empty. I guess I expected to see the 10 file names I had entered.
So I entered the 10 file names (I am really desperate by now) and tried running the gdal_merge command in the terminal window and got a series of error messages similar to:

  ERROR 4: `/Users/ericfoertsch/NewEngland/USGSBilFiles/n42_w074_1arc_v3.bil' not recognized as a supported file format.

Any ideas?


texts: Disadvantages to ESRI Plug-In Data Sources

texts:
I wanted to know what disadvantages ESRI plug-in data sources have when compared to custom layers.  I am looking to implement a fairly large (in terms of average dataset size) read-only custom type for ArcMap and ArcScene.

I have read the chart here.  I was wondering if there are any other disadvantages or major issues with plug-in data sources.

In particular, has anybody experienced issues related to performance, integration and interoperability with the Desktop tools, or any other significant issues?


texts: "Point in polygon" counts that respond to filters?

texts:
I have a list of addresses in one layer and a polygon layer of census block groups for the local town, all stored in postgis.  I have a nice display right now that labels each polygon with a count of houses within that block group by creating a postgis view along the lines of:

create view block_groups_count as
  select g.id, count(1) from
  block_groups g join addresses a
  on st_covers(g.geom, a.geom)
  group by g.id;

I can join this with the polygon layer itself and get the labeling I
want.  This works great with respect to modifying the polygon layer: I
can make changes, click "save", and everything updates automatically.

But now that I've got this far, I am looking for a mechanism by which
I can get the labels in the polygon layer to update dynamically in
reaction to changes to qgis filters applied to the address layer.

That is, if I apply a "show only blue houses" filter to the list of
addresses, I want to counts to reflect only the number of blue houses
within in each polygon.

Is it possible to set up this sort of behavior with qgis?

One alternative I've considered is to simply to create a new count
view for each attribute, but this results in several additional
layers, and switching between different views isn't particularly
graceful.


texts: Replacing 0 values of a raster with values from another raster in QGIS?

texts:
What I want to do is pretty straight forward: I have two raster with the same crs, resolution and extent. Both rasters show the number of days an area is flooded so they just have two values. Raster a consists of 0 and 183; raster b of 0 and 149. Logically, the area that is flooded only 149 days (b) is a little bigger than the area that is flooded 183 days (a).

Now I want to combine both maps. The formula I think should work using the saga raster calculator is:

ifelse(eq(a,0), b, a) .

In my opinion it translates to "If a pixel of map a is equal to 0 then use the pixel value of map b if not, keep the pixel value of a" So the result should be a map having just three values: 183, 149 and 0, right? But no, the resulting map has all kind of weird (i cannot explain them) values where map a is 0 and map b is 149. It seems like some transition from 183 over 149 to 0.

Why is that, or is it a bug?

Ps: Here is an example for one of the grey pixels of the calculated map: a=0, b=149 but the calculated map has a value of 128.30


texts: Summing pedestrians who pass through every road using Network Analyst with ArcMap

texts:
I have a network/routing challenge.
The main aim of this task is to estimate the flow of pedestrians/people who move on the road network from their start locations toward a final destination. The origins is a points layer where every point includes a field that counts the number of people who live in that point location. Further, we have the destination point (another points layer) where all these people will travel to via a road network layer. Using the closest facility tool with the road network dataset, origins and destination we could obtain the closest route from every origin to the destination as can be seen below:

After obtaining the closest facility above, we did an intersection between the original network (dark blue network) and the closest facility route (Purple network) to find how many times every segment has been traversed by every point. However, this method assumes that every origin point has the same number of people which is not true. Every point has a different number of people who live in that location/origin.  So instead of calculating the number of times every road has been passed by every route, we need to calculate how many people/pedestrians who passed every road towards the destination. As the closest facility route table result does not contain the number of people who are occupied in every point, I am thinking that there is a way (within Arc GIS) which can calculate the number of people who pass via every road. Personally, I think there is some way to calculate it within the closest facility to know how many people have traversed via every road.
How can do it?


texts:
I am trying to create custom mbtile for Cesium consumption.

I was able to create a mbtiles file in QGIS (zoom 0 to 16) using the QTiles plug-in. I tested the generated mbtiles file using mbviewer from mapbox and it worked fine. mbviewer

node cli.js --port 9000 ~/tmp/starpeak_low_res.mbtiles Listening on http://localhost:9000

So to serve it to Cesium, I used mbtiles-server from npm to create a WMTS compatible server. Here is the link to the server mbtiles-server

So in Cesium, all I had to do is add a a new imagery provider:

var high_res = layers.addImageryProvider(new Cesium.WebMapTileServiceImageryProvider({
    url : 'http://10.10.3.13:8030/starpeak_low_res/wmts',
    layer : 'starpeak_low_res',
    style: 'default',
    format : 'image/jpeg',
    tileMatrixSetID : 'GoogleMapsCompatible',
    maximumLevel: 16
}));

however, it takes a very very very long time for the images to load on Cesium but it's almost instant using mbviewer. The mbtiles file is 200 MB but I assume at different zoom levels only a portion of that file is fetched.

Just to make sure the issue is not with mbtiles-server, I downloaded the free low res sattelite mbtiles files from openmaptiles link here. link. That mbtiles file is served fine it's instant loading at different zoom levels.

for reference I uploaded the mbtiles file here. Credits to  NAIP GEOTIFF from earthexplorer.usg.gov generated mbtiles file


texts: Filter FeatureCollection by geometry type in GEE python API

texts:
I have a featurecollection with 141 objects inside (135 polygon, 4 geometrycollection and 2 multipolygon) an I have problems to indicate in a loop the region for images exportation in GEE, so I need to filter in the collection only the objects of polygon type for example for run the code easily.

table= ee.FeatureCollection("users/tpc/capas")
lista = table.toList(table.size())
m= lista.getInfo()

for i in range(0,141):

    predio_id = ee.Feature(lista.get(i)).getInfo()["properties"]["id"]
    predio = ee.Feature(lista.get(i)).geometry()
    predio_region= predio.getInfo()["coordinates"]

    S2_NDVI_Col= ee.ImageCollection('COPERNICUS/S2_SR').filterBounds(predio).filterDate(dateNDVI,today).filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 30)).map(addNDVI).map(addSAVI).map(addNDMI).map(addEVI)
    S2_NDVI = S2_NDVI_Col.sort("system:time_start")

    n_img=4000
    colListNDVI = S2_NDVI.toList(n_img)
    n= colListNDVI.size().getInfo()

    for j in range(0,n):
        NDVI = ee.Image(colListNDVI.get(j)) #obtiene la imagen n
        id_NDVI = NDVI.id().getInfo()
        x1 = len(id_NDVI)
        dateNDVI = id_NDVI[x1-38:8].replace('_','')
        int_dateNDVI = int(dateNDVI[5:7])
        print(dateNDVI)
        export = ee.batch.Export.image.toDrive(\
                                               image=NDVI.select("NDVI").float(),
                                               description='S2_NDVI_'+predio_id+"_"+dateNDVI,
                                               maxPixels= 220582361000,
                                               region=predio_region,
                                               folder='S2_FNDR_NDVI',
                                               scale = 10,
                                               crs='EPSG:4326')#,
        export.start()
    print('Descarga MODIS finalizada')


texts: Embedding interactive Google Maps into a Google Form

texts:
I have a Google form with some questions. People can answer the questions about (future) trees they want to have in their neighbourhood. But there also must be the possibility that people mark the place of this tree on a (google) map in the same Google form.
So my question is: how can I add (embed) a Google Maps into a Google form so that people can mark a location and so that the latitude and longitude are included together with the other answers of the form?
I really want to have the possibility that people can mark this location on a map, because the address is not detailed enough and it also costs a lot of time otherwise afterwards.
I searched on the internet but can't find the answer!


texts: Create a custom buffer in Google Earth Engine

texts:
I am working to make mosaics of Landsat 5 scenes that cover watersheds in Montana, and am currently struggling to deal with the "rough edges" on the east/west sides of scenes that get incorporated into the mosaics.
My current solution is to buffer the primary scene in the mosaic (the one that covers most of the watershed) by a negative amount sufficient to remove the edges of the scene prior to making the mosaic. This works in many cases but in some, it also removes valuable data from the top and bottom of the scene that could otherwise be used.
Is there a way to specify a custom buffer that would trim the scene by a variable amount depending on the side of the scene (i.e. trim the east/west sides by enough to remove the bad pixels but leave the north/south edges alone)?
===================== VISUALIAZTION PARAMETERS ================================
=============================================================================*/

// viz parameters for displaying scenes at Landsat 8 TOA scale
var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B3','B2','B1']
};

/*=============================================================================
====================== USER CONTROLS ==========================================
=============================================================================*/

var huc08 = "10050004" // ID of the watershed of interest at the HUC08 level

/*=============================================================================
====================== VECTOR DATA ============================================
=============================================================================*/

// Define bounds of the state of MT
var mt = ee.FeatureCollection("TIGER/2018/States").filterMetadata('NAME', 'equals', 'Montana').geometry()
// print(mt)
// Map.addLayer(mt)

// import all of the HUC08 boundaries and retain only those that intersect MT
var fc = ee.FeatureCollection("USGS/WBD/2017/HUC08").filterBounds(mt)
  .map(function(feature){
    return feature.intersection(mt)
  });
// print(fc)

var feature = fc.filter(ee.Filter.eq('huc8', huc08)) // identify individual HUCs by Name

var ROI = feature.geometry() // turn the desired watershed into a geometry object

/*=============================================================================
=================================== LANDSAT SCENES =============================
=============================================================================*/

var primary = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_038026_20090902')

var secondary = ee.Image('LANDSAT/LT05/C01/T1_TOA/LT05_037026_20090725')

Map.addLayer(primary, rgbVis, 'least cloudy - best scene', false)
Map.addLayer(secondary, rgbVis, 'least cloudy - best scene_1', false)

/*=============================================================================
=================================== MOSAIC ====================================
=============================================================================*/

var mosaic = ee.ImageCollection([secondary,primary]).mosaic() // default mosaic opertor = "LAST"
Map.addLayer(mosaic, rgbVis, 'mosaic l5 scene for roi', false)

/*=============================================================================
=============================== SOLUTION ====================================
=============================================================================*/

var fp = primary.geometry() // get the footprint of the primary scene
var fp_buff = fp.buffer(-8000)

var primary_clip = primary.clip(fp_buff)

var newMosaic = ee.ImageCollection([secondary,primary_clip]).mosaic()

Map.addLayer(newMosaic, rgbVis, 'clipped mosaic l5 scene for roi', false)

Map.centerObject(feature,9); // zoom the map to the HUC of interest
Map.addLayer({eeObject: feature, name: 'watershed bounds'}) // watershed of interest


texts: Trying to get right combination of Leaflet libraries

texts:
I've inherited a project which was using an ancient version of Leaflet and some associated libraries, and have now managed to add the Google Mutant to it. It used to use:

<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.1/leaflet.css"  />
<link rel="stylesheet" href="assets/vendor/leaflet.draw.0.2.1/leaflet.draw.css" />

And

<script src="assets/vendor/leaflet.draw.0.2.1/leaflet.draw.js" type="text/javascript"></script>

Now, I've upgraded a few things to:

<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet/1.6.0/leaflet.css" />
<link rel="stylesheet" href="/web/assets/vendor/leaflet.draw.0.2.3/dist/leaflet.draw.css"/>
<link rel="stylesheet" href="/web/assets/vendor/leaflet.label/leaflet.label.css"/>

And

<script src="https://maps.googleapis.com/maps/api/js?key=BlahBlahBlahCrazyWordsMoustache"></script>
<script src="~/assets/vendor/leaflet.1.6.0/leaflet.js"></script>
<script src="~/assets/vendor/leaflet.label/leaflet.label.js"></script>
<script src="~/assets/vendor/leaflet.draw.0.2.3/dist/leaflet.draw.js"></script>
<script src="~/assets/vendor/leaflet.googleMutant/Leaflet.GoogleMutant.js"></script>

This is all fine and dandy for the most part - the Google imagery is loading how I want it.

However, there are buttons present that add polygons and polylines, seemingly to a GeoJSON layer, and measure them, as well as functionality to edit and delete them. Now that functionality is seriously compromised - the label/popup which used to appear on mouseover is gone, only the first shape will edit, and nothing will delete. I do get an error message, which is:

> Uncaught Error: The provided object is not a Layer.
> at i.addLayer (leaflet.js:5)
> at i.showLabel (leaflet.label.js:9)
> at i._showLabel (leaflet.label.js:9)
> at i.fire (leaflet.js:5)
> at i._fireDOMEvent (leaflet.js:5)
> at i._handleDOMEvent (leaflet.js:5)
> at HTMLDivElement.s (leaflet.js:5)
> addLayer @ leaflet.js:5
> showLabel @ leaflet.label.js:9
> _showLabel @ leaflet.label.js:9
> fire @ leaflet.js:5
> _fireDOMEvent @ leaflet.js:5
> _handleDOMEvent @ leaflet.js:5
> s @ leaflet.js:5


texts:
I am using ESRI's spatial framework for Hadoop, in order to add spatial support to Hive.

One scenario that I find very often, is to have layers in different CRS, which I want to relate in the same query. Also, although most things are in WGS84 these days, sometimes it is useful to use a projected CRS, in order to measure things such as areas or distances.

Thus I find that reprojection is a very common need, when working with spatial data. Within PostGIS, I just apply ST_transform to one geometry, to transform it into another CRS. For instance, this code would reproject the centroid of Barcelona, from whatever crs to an UTM grid:

SELECT ST_Transform(ST_Point(2.16991941.387917), 5554);

In the ESRI's function list, I dont see any mention of ST_Transform. Am I missing something or this no way to achieve this using this framework? Is there any workaround around this?

I believe is a shame if one is not able to perform this operation on a scalable way, since doing it using another tool (for instance PostGIS), could be a bottleneck in a scalable workflow...


texts: Given a polygon layer, census tracts, how can I generate a table of centroid coordinates in QGIS, or add columns with the centroid coordinates or generate a new layer of polygon centroid points?


texts: I have a huge raster tif file 40GB, when I use this file in ArcMap it is very slow, I want to use this file only for the first levels, so I don't care about resolution. I don't want a compression I just want to reduce the resolution of the TIF. The tif has 0.3 cell size


texts: QgsMapCanvas() giving QCursor error (QGIS 3.16)

texts:
I am trying to load a standalone QGIS map window from a python script (using Python 3.7.0 64Bit interpreter in VS Code) but when I call QgsMapCanvas() the canvas shows for a fraction of a second and then I get the following error:

QCursor: Cannot create bitmap cursor; invalid bitmap(s)

Here is the code I am trying to run:
import os, sys

from PyQt5 import *
from PyQt5.QtGui import *
from PyQt5.QtCore import *

from qgis.core import *
from qgis.gui import *
from qgis.utils import *

QgsApplication.setPrefixPath(r"C:\OSGeo4W64\apps\qgis-ltr", True)
qgs = QgsApplication([], False)
qgs.initQgis()

canvas = QgsMapCanvas()
canvas.show()

qgs.exitQgis()

I deliberately removed all additions of layers from the script just to get the canvas to show, which it does if I run it from the python console in qgis but not from VS Code. It freezes when I try to add any content to it from the QGIS python console.
I am fairly sure I have added all possible environment variables for PYTHONPATH, PYTHONHOME, QT_PLUGIN_PATH, QT_QPA_PLATFORM_PLUGIN_PATH, and added to the PATH variables. Am I missing something really obvious?


texts: How to find nearest line from point using PostGIS?

texts:
I have a table of Roads in PostGIS which geometry is "geometry(MultiLineString,4326)".

Its attribute are

  gid   osm_id  name    ref type    oneway  bridge  maxspeed    geom

Now i want to know the closest Road from latitude and longitude.

Also I want to set some distance from that Line/Road like 20m both side and if the point is within that distance then it should give the name of that line else not.


texts: Feature not added to KML Vector Layer in Firefox

texts:
UPDATE simplified with JSFiddle example: here is the most simplistic JSFiddle showing what I mean:
http://jsfiddle.net/eflowbeach/75eq1mqe/

You can see in Chrome and IE that the feature appears, in Firefox nothing displays. The KML is valid.

Original Question:
Has anyone else had a KML vector layer display perfectly fine in Chrome (v35) and IE10+, but fails to add the features to the layer in Firefox (v30)? I'm using OpenLayers 2.12 and have also tested in 2.13.1 without success.

var test = new OpenLayers.Layer.Vector("test",
  {
    projection : new OpenLayers.Projection("EPSG:4326"),
    strategies : [new OpenLayers.Strategy.Fixed()],
    protocol : new OpenLayers.Protocol.HTTP(
    {
      url : "myFile.kml",
      format : new OpenLayers.Format.KML(
      {
        extractStyles : true,
        extractAttributes : true
      })
    })
  });
  me.map.addLayer(test);

When you look at the features of this layer in the console:

console.log(test[0].features);

Chrome and IE have 11 features in Firefox the feature array is empty.

Here's the link to the KML file in question:
https://drive.google.com/file/d/0B_iQDBiLH9cfQWIzS3A4anZqcVk/edit?usp=sharing

I've had this happen to other files too, so I was wondering if it was something wrong with my format statement.


texts: Add items to map, collectively, one at a time, per attribute ID

texts:
I have successfully been able to filter, for instance shakemaps, by attribute id - I have successfully achieved this multiple ways - but all ending with the same problem.

With the where filter
Definition expression
Iterating through all attribute ids and bringing them back.

The problem: All exists for only allowing/displaying one per attribute id at a time... My goal is to feed the attribute ids into a checkbox list (which I have done), but allowing for items via attribute id to be added to the map as they are checked, collectively, one at a time - currently I can not seem to get this to work with the aspect of having multiple or more then one appear at a time on the map.

i.e. the below filter (attempted logic 1) & CodePen:
.....
function filterByID(event) {
const selectedID = event.target.getAttribute("data-id");
eqLayerView.filter = {
where: "id = '" + selectedID + "'"
};
}
     view.whenLayerView(fl)
       .then(function(layerView) {

       eqLayerView = layerView;

       eqLayerView.filter = {
         where: "id = ''"
       };

.............

i.e. another attempted logic (adding multiple at a time here, line by line, or via array):

    layer.definitionExpression = "id = 'us70008jr5'",
    layer.definitionExpression = "id = 'cgi988jr52'",

i.e. 3rd attempt with a suggestion here on GIS exchange: Loop through attribute ids of FeatureLayer

    layer
      .load()
      .then(() => {
        // create a query from the layer
        const query = layer.createQuery();
        query.returnDistinctValues = true;
        query.where = "grid_value > 2"; // or 1=1 if you need them all
        // the field you want distinct values for
        query.outFields = ["id"];
        return layer.queryFeatures(query);
      })
      .then(({ features }) => {
        // extract the ids to a list
        const ids = features.map(({ attributes }) => attributes.id);
        return ids;
      })
      .then((ids) => {
        // You can store them how you want, in this case,
        // I put them in a dropdown select menu
        const fragment = document.createDocumentFragment();
        ids.forEach((id) => {
          const option = document.createElement('option');
          option.value = id;
          option.innerText = id;
          fragment.appendChild(option);
        });
        list.appendChild(fragment);
        map.add(layer);
    });

All attempted logic above result in the toggling of a shakemap by attribute id to be displayed only one at a time — by toggling a new on, the previous turns off, I need the ability to have multiple being able to exist on the map at once.


texts: Map Publishing in ArcGIS Pro gives errors 00216 and 00230?

texts:
I'm receiving two errors (00216 & 00079) (below) that need resolution before I will be able to share a map I've created to ArcGIS Online (from ArcGIS Pro's interface). The map document was originally created in ArcMap 10.5(?) and then imported into Pro. I'm using the Canada Lambert Conformal Conic projection with the input layer.

I can't for the life of me seem to get past either. I've added a few different base maps but that just yields a different error (00216 becomes 00230) and does not resolve 00079.


texts: Displaying WFS Features Based on Zoom Level

texts:
I have a WFS feature layer with ~25,000 features and have been experiencing significant performance issues with my OpenLayers map. My first approach to solving this has been to pre-compute feature clusters on the server side, instead of relying on the clustering in OpenLayers.

I wrote a script to cluster my features based on the resolution at each zoom level, and I placed the clusters for all of the zoom levels in my PostGIS database. My hope was to use a filter like the one below to select only those features in the database that are tagged with the zoom level equal to the current zoom level of the map:

function init() {

        var zoomLevelIndicator;
        var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
        renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;

        map = new OpenLayers.Map ("map", {
            controls:[
                new OpenLayers.Control.Navigation(),
                new OpenLayers.Control.PanZoomBar(),
                new OpenLayers.Control.Permalink(),
                new OpenLayers.Control.ScaleLine({geodesic: true}),
                new OpenLayers.Control.Permalink('permalink'),
                new OpenLayers.Control.MousePosition(),
                new OpenLayers.Control.Attribution()],

            maxExtent: new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34),
            maxResolution: 156543.0339,
            numZoomLevels: 19,
            units: 'm',
            projection: new OpenLayers.Projection("EPSG:900913"),
            displayProjection: new OpenLayers.Projection("EPSG:4326"),

            layers: [new OpenLayers.Layer.OSM(),
                     new OpenLayers.Layer.Vector("Filtered by Zoom Level", {
                        strategies: [new OpenLayers.Strategy.BBOX()],
                        protocol: new OpenLayers.Protocol.WFS({
                            version: "1.1.0",
                            url: "http://localhost:8080/geoserver/wfs",
                            srsName: "EPSG:900913",
                            featureNS: "http://localhost:8080/myworkspaceuri",
                            featurePrefix: "MyWorkspace",
                            featureType: "clusteredtestpoints"
                        }),
                        renderers: renderer,

                        filter: new OpenLayers.Filter.Comparison({
                            type: OpenLayers.Filter.Comparison.EQUAL_TO,
                            property: "zoomlevel",
                            value: zoomLevelIndicator

                        })

                    })

                ]

        } );

        map.events.register('zoomend', this, function (event) {
            zoomLevelIndicator = map.getZoom();
        });

            var switcherControl = new OpenLayers.Control.LayerSwitcher();
            map.addControl(switcherControl);
            switcherControl.maximizeControl();

        if( ! map.getCenter() ){
            var lonLat = new OpenLayers.LonLat(lon, lat).transform(new OpenLayers.Projection("EPSG:4326"), map.getProjectionObject());
            map.setCenter (lonLat, zoom);
        }
    }

Specifically, I am trying to update the zoomLevelIndicator variable each time a zoom event occurs. When I test this, the features for the initial zoom level show up, but do not update when I change the zoom level.

Does anybody have any ideas as to why this is not working?

Thanks,

-Rob


texts: Assigning rubberband as polygon using PyQGIS

texts:
I am trying to display a rubberband on my mapCanvas in QGIS but it seems that setToMapCanvas will only accept a QRect object. I am trying to show on the canvas the visible area of a maplayout with a rotation. I tried rotating the QRect with QTransform but I can only get the rotation to rotate from 0,0 (i.e. Top Left Corner) or a wildly random location depending on scale. Or I transform the QRect with a rotation and return it to a PolygonF only to find no way to set it to the map canvas with rubberband.

When I try to return the Polygon I get TypeError: setToMapCanvas(self, QTransform): argument 1 has unexpected type 'QPolygonF'
When I try to transform the QRect I get TypeError: setToMapCanvas(self, QTransform): argument 1 has unexpected type 'transform'

mtp = self.iface.mapCanvas().mapSettings().mapToPixel()
p1 = mtp.transform(QgsPoint(rect.left(), rect.top()))
p2 = mtp.transform(QgsPoint(rect.right(), rect.bottom()))
rect = QTransform().mapToPolygon(QRect(p1.x(), p1.y(), p2.x() - p1.x(), p2.y() - p1.y()))
rect_r = rect.rotate(45)
return QPolygonF(rect_r)

When I use this code to rotate the rubberband

self.__cleanup()
extent = self.mapitem.extent()
center = self.iface.mapCanvas().extent().center()
self.corner = QPointF(center.x() - 0.5 * extent.width(), center.y() - 0.5 * extent.height())
self.rect = QRectF(self.corner.x(), self.corner.y(), extent.width(), extent.height())
self.mapitem.setExtent(QgsRectangle(self.rect))
self.rubberband = QgsRubberBand(self.iface.mapCanvas(), QgsWkbTypes.PolygonGeometry)
self.rubberband.setToCanvasRectangle(self.__canvasRect(self.rect))
#    tveinot 20200122 set the rubberband rotation and rotation point
self.rubberband.setTransformOriginPoint((0.5 * self.rect.width()),(0.5 * self.rect.height()))
self.rubberband.setRotation(self.dialogui.spinbox_Rotation.value())
self.rubberband.setColor(QColor(127, 127, 255, 127))

The rotation point varies depending on my scale and location.
If I change (0.5 * self.rect.width()) to center.x()/y() respectively 1 degree of rotations sends the rectangle off screen.

What am I doing wrong?


texts: Error generating chart: No features contain non-null values of "system:time_start" Create NO2 chart with unit ton/m2 in GEE Sentinel-5p

texts:
I'm using the below code to create NO2(Sentinel-5p) time series charts for my region of interest by mean reduction. However, upon charting, I get an error message
Error generating chart: No features contain non-null values of "system:time_start"
Please, somebody advise me how to fix this error?
Here is my code:
    var geometry =
        /* color: #d63000 */
        /* displayProperties: [
          {
            "type": "rectangle"
          }
        ] */
        ee.Geometry.Polygon(
            [[[112.19326448481443, 37.22101939910812],
              [112.19326448481443, 36.467942373936516],
              [113.19250536005858, 36.467942373936516],
              [113.19250536005858, 37.22101939910812]]], null, false),

    // Create a FeatureCollection from the list and print it.
    var geometry = ee.FeatureCollection(geometry);
    print(geometry);
    Map.centerObject(geometry);
    Map.addLayer(geometry);

//Importing image and geometry:
    var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
      .select('NO2_column_number_density')
      .filterDate('2019-01-01', '2019-01-28')
      .filterBounds(geometry);

    // calculation mol/m^2 to ton/m^2.
    var mol_to_ton = collection.map(function (image) {
      return image.multiply(0.000046);
    });

    //Visualization parameters
    var band_viz = {
      min: 0,
      max: 0.00000002,
      palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
    };

    Map.addLayer(mol_to_ton.mean(), band_viz, 'S5P N02_ton');

    // Create and print the chart.
    print(ui.Chart.image.series(mol_to_ton.mean(), geometry, ee.Reducer.mean(), 1113, 'system:time_start'));


texts: How to convert an Erdas training set to .csv?

texts:
I am working with Erdas Imagine's Signature Editor to perform maximum likelihood classification.  I need to use the same training set created in Erdas for other classification using R.  Is there a way to extract the raster band pixel values under the .aoi polygons and export that information to a .txt or .csv?  I am looking for the same functionality as the ArcGIS Sample tool where the individual raster band pixel values are extracted underneath points and written to a table (e.g. Figure 2).  Alternatively, is it possible to export the signature editor table (i.e. including R,G,B,NIR bands and the class label) to a .csv?  Any workarounds on how to integrate Erdas's signature editor with the analytical functionality of R would be very helpful.

For example, I need the functionality of the signature editor

Figure1

Yet I need to convert these values to a format that I can read into R.  Note that I need all four bands of the input image along with a class ID (i.e. grid code).  This table was generated in ArcGIS.

Figure 2


texts: How do I install PyProj module for Python 2.7 on Windows 7?

texts: I am having trouble installing the pyproj module for Python 2.7 on Windows 7.  I have tried both the 32 and 64 bit installers from this list but the install crashes each time I run it.  I have the Python 2.7 version that installs with ArcMap 10.2.2 loaded on my machine.  I need pyproj to run the Scanaerial plugin in JOSM to locate water sources in my county.  Any ideas on what could be causing the crash, or how to install pyproj properly?


texts: Is there any efficient way to create a raster with pixel value representing a unique id calculated based on pixel lat/lon?

texts:
I have 10m raster for each state in US. Current pixel value represents something. I want to have 'same' raster but pixel size representing a unique id calculated based on pixel lat and lon. I am wondering whether there is any efficient way to do it, ArcGIS existing function or python C# programming are both fine.

I cannot do anything with Raster Calculator in ArcGIS since the pixel value of current raster has nothing to do with unique id or lat/lon.

So far, I tried to create lat/lon pair and calculate such id with programming. But in order keep a connection with the raster, I also read current pixel value out. So I have one record for one pixel, formatted as lon, lat, pixelVal, uniqueId, and I write it into a csv file. Since for 10m raster, the number of pixels is very very large. I am wondering whether I can do it on raster level, instead of using this way. Thanks in advance!


texts: How to redistribute population data using Python in ArcGIS?

texts:
I am trying to redistribute population data from census centroids among to housing density. Would you please help me to do that?
I do have census area centroids (points) and the boundaries belong to each centroid. Point data has population value. Also I have a raster set housing density. I would like to redistribute population data (points) in their own site-boundaries among housing density. Output data needs to be a raster data 100×100 meter pixels.

I am familiar with FORTRAN and VB but I am not a programmer. I have access to ArcGIS 10.1 however I have never used Python.


texts: arcpy clip_managment migration from arcpy 10.1 to 10.2

texts:
i've created script using arpy with arcgis 10.1. One of the goal is to clip raster via shp file, which contains one feature. it worked perfect.
after migration to arcgis 10.2 i've noticed some problem. result of clipping raster in exact shape of feature. it is just clip to extent of clipping feature.

arcpy.Clip_management(inFC, "#", cut_director + "/" + outFC, shp + "/" + shp_clip, "#", "ClippingGeometry")

i've tried to do clipping via arctoolbox (it  does the job) and then copy the process as python snippet but after changing input, output and clipping feature i received error.

thanks in advance for help

edit: i've found answer. when change nodata parametr from "#" (second one in code above) to "256" it works.


texts: Pixelated hillshade from 5m elevation data

texts:
I recently tried to create a hillshade of a valley in the Scottish Cairngorm mountains using some Digimap 5m elevation data (specifically the OS Terrain™ 5 DTM).

The data was downloaded in .asc format and I converted this to raster using the ArcToolbox Ascii -> Raster function. Then I mosaic'd the individual tiles together and used the hillshade tool to produce this:

The view is approximately 10 miles across from left to right, should it be this fuzzy?

I'm very new to ArcMap and GIS in general.


texts:
According to this: http://www.gdal.org/ogr/drv_mitab.html the TAB file writer sets default bounds, with the caveat that 'we should fix this'.

Bounds causes issues with respect to coordinate precision, see here: http://www.dse.vic.gov.au/property-titles-and-maps/maps-imagery-and-data/data/vicmap/datums-and-coordinate-systems/victorian-bounding-box

I would like to know if anyone knows how I can set bounds using OGR, and specifically, with GeoServer's WFS export using the OGR driver.


texts:
I would like to create a marker line in which a font marker is placed within a dashed line space at equal intervals.  Something like line-style-with-alternating-dots-and-long-dashes

c - c - c -

I can get this to somewhat work by creating a marker line with a custom dash pattern:
Dash=10 Space=3 mm
Then setting the font marker with a interval of 13mm and offset of 11.5 (midpoint of space).

The problem I have is that it appears that when the line become more complex there appears to be a misalignment of the dash space and marker offset

Regards


texts: Cannot see OpenLayers under Plugin in QGIS?

texts:
Although I have installed the OpenLayers plugin in QGIS 2.4, I cannot see it in the list when I click on Plugins.

Does anyone else have this issue?


texts: How to get size of SDE geodatabase feature class

texts:
Is there a way in catalog or arcpy to determine the amount of hard drive space a SDE feature class takes up on disk?

In catalog the 'size' column is blank and Esri help states:

  The size, modification date, and read/write columns apply to several file-based objects, such as maps and shapefiles. No values will appear in these columns for data stored in a database or file-based objects that don't support retrieving that information from the file system.

As for arcpy, I have not see a property of any Describe object that can call size or read/write for that matter.  The answer question How to get size of file geodatabase feature class on disk? provides the solution for a file GDB, but how can it be done for SDE feature classes?


texts: Creating new field based on other layers' attributes - "for" cycle error

texts:
I'd like to create a new field in vector layer A, based on attributes in vector layers B and C but PyQGIS returns a series of wrong numbers, particularly just the first feature value of the field C (I suppose the error is in the for cycle's script/sequence). Could you help me?
featuresB=layerB.getFeatures()
featuresA=layerA.getFeatures()
featuresC=layerC.getFeatures()
for fA in featuresA:
    for fC in featuresC:
         for fB in featuresB:
             if fA.geometry().area()*2>=fB.geometry().area():
                print fC[field_x]
             else:
                print fA[field_y]


texts: Should we always calculate length and area in lat,lng to get accurate sizes/lenghts of real world measurements

texts:
Heres is the simple case, that we discussed at the office today.

Given a map in local projection where units are meter. Since we are in denmark and we have some very local projections at our hand I will make an example using two projections here.

First one: http://epsg.io/32632  (UTM 32, a very large local projection)

Second one: http://epsg.io/4097 (DKTM1, a very local projection).

The first Point of interest in both coordinate system:

http://epsg.io/map#srs=32632&x=511397.912332&y=6177295.039367&z=21&layer=gmaps-satellite
http://epsg.io/map#srs=4097&x=211402.25&y=1179643.35&z=21&layer=gmaps-satellite

and the second point of interest in both systems:

http://epsg.io/map#srs=32632&x=508321.748862&y=6177119.058277&z=21&layer=gmaps-satellite
http://epsg.io/map#srs=4097&x=208324.91&y=1179467.30&z=21&layer=gmaps-satellite

Now that we have defined the two points in both system we read out the coordinates for both in both systems and calculate the length between them using Pythagoras.

Distance in DKTM1 (EPSG:4097) is  : 3082.372235
Distance in UTM32 (EPSG:32632) is : 3081.193264

So the next question is, what is the accurate actual distance in the real world if we send out a guy on the runway to measure the distance between these two points?

How would i calculate this using mathematics? Reproject to lat,lng and do something there?

I will buy a beer for the one that provides the best example, and you can pick it up at either of the two locations :)


texts: How can I remove rows that have duplicate cells in 2 dif files - while choosing which file has the dup removed?

texts: I have 2 files and I merged them together into one file-geodatabase feature class. I want to remove entries from File1, Column C, that also appear in File2, ColumnD. I attempted to use the geoprocessing tool DeleteIdentical but it deletes from EITHER File1 or File2 (aka Column C or Column D). I need to specify that only rows from File1 (Column C) should be deleted. I can distinguish in the merged File1_2_Merge which column belongs to which file. I am using ArcGIS 10.2 Desktop. All 3 files are file-geodatabase feature classes.


texts:
Recently I started using DB Manager to query Spatialite DB. I was happy to see that QGIS DB Manager had a "save query" option, but it seems it saves them only for the current session of QGIS. When I reload the program and reconnect to DB all my saved queries are gone.

So, is this the way it's suppose to work? I was hoping that the queries were being saved inside the DB.


texts: Replacing layer data sources using replaceDataSource of ArcPy?

texts:
I'm attempting to replace a datasource in an outline map frame using ArcPy. The rest of the script seems to run OK, but after the script has run the datasource the datasource of that layer shows an error. A section of my code is below:

for df in arcpy.mapping.ListDataFrames(mxd):
    if df.name == 'Overview' or df.name == 'Keyplan':
        df.scale = 500000
    for lyr in arcpy.mapping.ListLayers(mxd,"",df):
        if lyr.datasetName == 'blah':
            arcpy.mapping.RemoveLayer(df,lyr)
            lyr.replaceDataSource(r'\\directory','FILEGDB_WORKSPACE','FCName')
            arcpy.mapping.AddLayer(df,lyr,"TOP")

arcpy.RefreshTOC()
arcpy.RefreshActiveView()
fullpath = os.path.join(fileDir,'ALTERED_' + fileName)
mxd.saveACopy(fullpath)
del mxd

I have had a look at http://forums.arcgis.com/threads/23803-issues-updating-data-source-via-arcpy.mapping?p=78411&viewfull=1 and used some of the tricks there, but I can't seem to make it work. Picture of layer datasource

before:

and after

You'll see it has the correct layer name, but no actual layer. Any ideas what could be causing this. I thought at one point relative paths, but apparently they do not matter when using UNC paths, as I am.

This appears to work when I do it in ArcGIS itself using the "CURRENT" map, so pretty sure it's not a pathname error.

Any ideas?


texts:
I tried to convert one PDF layer to tiff using GDAL but the resulting tif still has all layers.
This is what I tried.
gdal_translate -of GTiff burglesum.pdf burglesumout.tif --config GDAL_PDF_LAYERS "Layer.Image"


texts: The best way to convert points to polygons using ArcObjects

texts:
I have a set of points from GPS I want to convert them to polygons..I have two approaches to do this.1- convert points to polygon..I do this with below code:

        private IPolygon PolygonFromPoint()
    {

        try
        {
            IGeometryBridge2 pGeoBrg = new GeometryEnvironment() as IGeometryBridge2;
            IPointCollection4 pPointColl = new PolygonClass();

            int numPoints = pList.Count;
            WKSPoint[] aWKSPointBuffer = new WKSPoint[numPoints];
            for (int i = 0; i < pList.Count; i++)
            {
                WKSPoint A = new WKSPoint();
                A.X = pList[i].X;
                A.Y = pList[i].Y;
                aWKSPointBuffer[i] = A;
            }
            pGeoBrg.SetWKSPoints(pPointColl, ref aWKSPointBuffer);
            IPolygon pPointColl2 = pPointColl as IPolygon;
            pPointColl2.SpatialReference = spatialrefrence;
            OutputPolygon = pPointColl2;
            status = true;
            return pPointColl2;
        }
        catch (Exception e)
        {
            MessageBox.Show("PolygonFromPoint: " + e.Message);
            status = false;
            return null;
        }
    }

This works some times but some times I get the error that says the polygon is not colosed.I could not manage to fix this problem

Another way is converting to closed polyline and then convert them to a polygon.I managed to convert them to polyline using geoprosseor but can not convert it to polygons.(I can not use feature to polygon either because of licence problam).

So what ways do you suggest me to do this task?I can not find any samples or something like that.Can you please help me do this?

I am using ArcObjects 10.2.1 and C#
thanks a lot


texts: Uploading a map to Qgiscloud seems to work fine. But uploading data presents two problems. Qgiscloud complains that my version of OGR/GDAL is outdated. And if I try to upload data, qgis complains that it is unable to open datasource. I'm using Linux Mint 13 and both Qgis 1.8.0 and 1.9.0 on different computers. OGR/GDAL is 1.7.3 and I see no way to update that. Any suggestions? Thanks.


texts: Multi Value Hot Spot Analyis in ArcGIS Desktop?

texts:
Is there a geoprocessing tool that will find areas that are hot or cold spots in multiple fields in the table?

I have a polygon layer with two fields. One field is the average organic soil depth for bogs and the other field is the average organic soil depth for fens. I want to know if there are any areas where both fields spatially have statistically high values compared to their respective means and vise versa.

I managed to do two separate hot spot analysis on each field and made note of the polygons that both had a hot spot, but is there a one stop tool or analysis specifically for this task?

I'm currently using a geodatabase in ArcGIS 10.3, but I'm open to using other programs and tools.


texts: Getting Intellisense / AutoCompletion for ArcGIS API for JavaScript 4 in Visual Studio 2015?

texts: In Visual Studio 2015, how do I get intellisense / autocompletion for the ESRI ArcGIS API for JavaScript version 4.3. I haven't been able to find any documentation that explains how to add intellisense.


texts: How can I display GeoEye1 imagery in colour in Opticks?

texts: Anyone familiar with the use of the open-source remote sensing software Opticks? Im new to opticks, I load GeoEye1 imagery but it displayed as greyscale. It supposed to be displayed as bands combination true colour or false colour composites.


texts: Append to text file using ArcGIS and Python

texts:
I have wrote a small routine using Python for ArcGIS (9.3) that will get the count of 2 data layers, calc a percentage, and then write the text back to a text file (CSV) to create a sort of "report" that can be opened in excel and formatted.  I want to be able to run this routine on many different dataset pairs and append the data to the end of my CSV report.

The problem is that I am having is that when I attach this Python script to an ArcGIS Script Tool, it automatically erases the file.  I am testing that with

if os.path.exists(file):

right off the bat and this response is always "file does not exist".  I have tried setting

gp.OverwriteOutput = 0

with no success, and am appending to the file using

fwrite = open(file, 'a')

Does anyone have a workaround that will allow me to choose an output file but NOT overwrite the output?


texts: Spatial extent not showing in Geonetwork metadata view

texts:
I have a problem displaying spatial extent in Geonetwork when I open some metadata record (as you can see from the image example).
I get the right extent I previously defined in metadata editor, but when I come to the metadata view page, basemap doesn't show up.

There was similar problem explained here in gis.stackexchange, but that doesn't work for my case.

I am using Geonetwork 3.04, just installed it on Windows XP. It is not connected to any database nor server, cause I'm not very familiar with computing and everything (I just needed Geonetwork as a catalog for spatial data).


texts: I have a tab delimited file and a shapefile which demonstrates the forest areas an area. When I open them in a project the shapefile covers the points of my samples. Is there a way to counteract this? This is done regardless of which file I open first.


texts: Overview, multiple rasters in map canvas

texts:
I work with multiple rasters and merging them into a VRT isn't an option. Is there a way to get info on what raster layer I'm looking at as a "mouse over" event within map canvas?

I'm on version 2.18 and 3.x as well.


texts: Parting lines based on intersection with polygons

texts:
I have two tables:

roads containing one row (line) for each road with column road_code uniquely identifying the name of the road;
districts containing one row (polygon) for each district with column district_code uniquely identifying each district.

I want to create a new table of lines called roads_by_district, where each road from the table roads is parted based on the intersection with the districts.

So table roads_by_district should look like

road_code | district_code | geom(line)
H1          District_A      0105...
H2          District_A      0105...
H1          District_B      0105...
H2          District_B      0105...

I guess I should use a query with ST_Intersects and ST_Contains but I don't want to return the entire road H1 that partially intersects with district District_A, but only the part of H1 that is totally contained by District_A.


texts:
I want to add the functionality of switching the option to pan the map and put it in a checkbox. But I'm not quite sure how I can do that.

The reason I want to do this is: I have a layer which displays features and has a select and hover control. So there's no way to pan the map, because the clickEvent triggers direct action from the Handlers (click,hover). I want to add a checkbox where I can activate/deactivate the possibility to pan the map, so I can pan the map when there are features present, which normally would trigger a reaction to the mouseover the feature.

The way I've tried to do it is this one but it didn't work.
I attached the function to the checkbox but it never worked.

function allowPan() {

    if(panAct == false){
            highlightCtrl.activate();
            selectedControl.activate();
            alert("controls activated!!")
            panAct = true;
        }
    if(panAct == true){
        highlightCtrl.deactivate();
        selectedControl.deactivate();
        alert("controls deactivated!!");
        panAct = false;
        }
}


texts: Creating one line of best fit for multiple overlapping .gpx files

texts: I am new at QGIS and am trying to create a map of ski trails. I mapped out all the trails but in some places I passed by the trail more than once. On the .gpx file it shows as two separate lines and am wondering if there is a way to merge the lines into one? I want to try and create a line of best fit of all the trips I have done.


texts: Transformation from Sphere_Natural_Earth to WGS1984

texts:
I have a set of files that are stored in Tom Patterson's Natural Earth Projection, and I would like to convert them to WGS1984. This will require a transformation, but I can't locate one and don't know where to start with building one, if that's the route this needs to go.

I seek help with methods for reprojecting.

To avoid any confusion, this is NOT the same projection that the Natural Earth data comes in (see the link for disambiguation).

So that you can see the Datum etc., the .prj string is:

PROJCS["Sphere_Natural_Earth",GEOGCS["GCS_Sphere_GRS_1980_Mean_Radius",DATUM["D_Sphere_GRS_1980_Mean_Radius",SPHEROID["Sphere_GRS_1980_Mean_Radius",6371008.7714,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Natural_Earth"],PARAMETER["False_Easting",90.0],PARAMETER["False_Northing",90.0],PARAMETER["Central_Meridian",180.0],UNIT["Meter",1.0]]