texts:
I would like to extract a MODIS land cover dataset, among other things from DAAC. I found this python tool which allows me to feed a URL and automatically download the file, and also this library which allows opening hdf datasets. I aim to later convert hdf datasets extracted from DAAC with gdal or some other python tool.

However, I have a problem. I have no idea what the index of the DAAC database is showing. I see the following abbreviations at the index but have no idea what they mean, and which one of them can lead me to the dataset I want.

ASTT
COMMUNITY
ECOSTRESS
MEASURES
MOLA
MOLT
MOTA
VIIRS

Could anyone help me navigate through the DAAC database to get to what I want? Is there a document somewhere that can explain what these things mean?


texts: Using field as input for multiple ring buffer in ArcGIS Desktop?

texts:
I have a feature class with points and I need to create two concentric circle with different distances for each point. Moreover the two surfaces do not have to overlap.

Can I use "multiple ring buffer"? And how I am able to add as an input features a field that  contains the distance to buffer each feature?


texts: SLD - styling geometry collections in geoserver

texts:
I’m trying to style a geometry collection column (PostGIS)  that stores both polygons and points in each record.

I already read the geoserver manual on how to style mixed geometry types and I don’t think this is possible when you have mixed geometry types in the same record.

I’m using geoserver 2.12.0.

Can anyone please help?


texts: How to change the default resampling method in ArcGIS 10.2.2 to Create Referenced Mosaic Dataset

texts:
When creating a reference mosaic dataset in ArcGIS 10.2.2., the default resampling is Nearest Neighbor but I'm using continuous data so I need to change the resampling method.

When using the Create Referenced Mosaic Dataset tool, I don't see an option for changing the default resampling method.

I've tried a couple things, in Environment Settings > Raster Storage, it seems you can only change the pyramid resampling method. In Advanced ArcMap Settings, changing the default resampling technique to something other than Nearest Neighbor doesn't seem to work either. After changing these both to Bilinear, and re-running the Create Referenced Mosaic Dataset tool, checking in the results window I still see Resampling Method: NEAREST.

For some more background information, these are the steps I have taken: 1) Create a file geodatabase 2) Create a raster catalog 3) Load rasters into the catalog with the load from workspace option 4) Create a referenced mosaic dataset.


texts: Open-source tools for run-off depth calculation

texts:
I am working on modeling surface run-off (in ArcGIS 10.2) due to a rain event (e.g., 5 year rain for a duration of 30 mins.) over an approx. area of 30 km2 for which I have a DEM and necessary soil as well as land-cover data (Manning's n values).

I want to build my own model, but in the meantime want to ask if any open-source tool calculates run-off depth at defined time intervals? The final result would be to generate run-off depths at certain time intervals and generate a movie out of it. Would be great if you could share your experiences with the available open-source hydrological models.


texts: Is 2Gb RAM enough to run ArcGIS for Desktop?

texts:
I am on a budget and I want to buy a laptop solely for the purposes of using Arcmap and nothing else.

The laptop I am going to buy is as follows:

HP Stream 11-d061na 11.6
Intel Celeron N2840 dual core processor.
2.16 GHz/2.58 GHz with Burst
1 MB cache
2GB RAM.
32GB solid state hard drive.
Microsoft Windows 8.1.
11.6 inch screen.
High definition display.
Resolution 1366 x 768 pixels
SD memory card reader

I think that it fits most of the criteria as listed on the website requirements (http://desktop.arcgis.com/en/desktop/latest/get-started/system-requirements/arcgis-desktop-system-requirements.htm#GUID-0D22816B-298C-4ABF-BC4E-5C32D605F055). But before I make the purchase I need a second opinion incase I am wrong!


texts: Looping through records to create ArcGIS reports

texts:
I am using ArcGIS 10.1 Standard with Python 2.7. I am trying to loop through records in a feature class using a Search cursor and for each record, run a custom ArcGIS report. I have created the report layout file (rlf) and can make this work outside a loop. The problem I have run into when I try to run within the loop is that the Search Cursor works with feature classes while the ExportReport function works with layers. I've tried using the MakeFeatureLayer function to make the conversion from feature class to feature layer within the loop, but seem to be missing something in the arguments for MakeFeatureLayer. This is the code I have used successfully to create the one-off report.

import arcpy
mxd = arcpy.mapping.MapDocument(r"G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\StreamCrossingReports\Figures\MXD\AGOL\CreateReports.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
fccrossings = arcpy.mapping.ListLayers(mxd, "", df)[0]
arcpy.mapping.ExportReport(fccrossings,
                       r"G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\StreamCrossingReports\Reports\ArcGIS_Report_Templates\Report_SCRField_v6.rlf",
                       r"G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\StreamCrossingReports\Reports\SCR\Test\Test.pdf", "SELECTED")

del mxd

And this is what I have come to using the loop.

import arcpy

fccrossings = r"G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\StreamCrossingReports\Download_AGOL\NF\SCR_04Nov2015_ReducedPhotos\366c11ade3cb4b95a3772296a1b0d8f1.gdb\StreamCrossingNF"

with arcpy.da.SearchCursor(fccrossings, ["OBJECTID"], '"OBJECTID" = 1337') as cursor:
    for row in cursor:
        report = arcpy.MakeFeatureLayer_management(row, "templayer")
        arcpy.mapping.ExportReport(report,
                    r"G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\StreamCrossingReports\Reports\ArcGIS_Report_Templates\Report_SCRField_v6.rlf",
                    r"G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\StreamCrossingReports\Reports\SCR\Test\Test.pdf")

del cursor

This is the error message I receive...

Traceback (most recent call last):
  File "G:\505573\MGT\40-ENG\4E-ENV\GIS\Misc\Python Scripts\CreateSCR_v2", line 7, in
    report = arcpy.MakeFeatureLayer_management(row, "templayer")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.1\arcpy\arcpy\management.py", line 5748, in MakeFeatureLayer
    raise e
RuntimeError: Object: Error in executing tool


texts: Classifying coordinates on land/on water using QGIS very slow processing issue

texts:
I'm very new to QGIS. I have a problem related to the "select by location" tool of QGIS.

I have loaded 10 different layers into my map containing data of waters (rivers, canals etc.). Now, I want to classify my own data, consisting of ID's with latitude/longitude, as "on land" or "on water". I have chosen the intersect function of QGIS to solve this problem. However, my data contains millions of datapoints which causes the process to be very slow. I have checked the "Use only selected features" checkbox, and the plan is to repeat this step for every single layer. The first layer is "creating a new selection", the next layers should be "add to current selection". Now, the problem is that it takes more than 1 hour for the first layer already. Is there a way to speed up the process? Or I am doing something wrong?


texts: Approximate distance between two points (longitude / latitude) without Haversine

texts:
How can I judge from the long/lat coordinates of a point A alone whether it is possibly/certainly not within a short distance of x meters from another coordinate B, using R (without applying geosphere::distHaversine or any other advanced calculation)?

"Certainly" in this case means that the long/lat of A tells me that it cannot be within the critical distance to B, while "possibly" means the complementary case.

Example:

x <- 100
point_A <- c("long"=8.678180,"lat"=50.114390)
point_B <- data.frame("id"=c(1,2,3,4),"long"=c(8.678459,8.618162,8.610878,7.208705),
"lat"=c(50.114124,50.226831,50.230530,50.799290))
point_B_certainly_not_in_distance <- subset(point_B, ...) # wanted
point_B_possibly_in_distance <- subset(point_B, ...) # wanted

In this example, only coordinate with id 1 should be identified as possibly within distance, whereas id 2, 3 and 4 should be identified as certainly not in distance.

To clarify the purpose: Imagine that point_B has >10^6 rows. In that case, a simplified "pre-subsetting" before applying the Haversine formula makes great sense.


texts: Gdalwarp creates a "perspective drawing" effect rather than an ortho-rectified image

texts:
I have been happily using the "utilities as libraries" functionality of gdal 2.1+ to ortho-rectify an image. First I use gdal.Translate to convert the source file to a geotiff, and then I warp it to (supposedly) match the coordinates. Here is the code:

tmp = gdal.Translate("tmp1.tif", src_img, format="GTiff", outputSRS="EPSG:4326", GCPs=gcpList)
warped = gdal.Warp("warped1.tif", tmp, srcSRS="EPSG:4326", dstSRS="EPSG:4326")

However, when I open warped1.tif at the end, I get something that looks like this:

The strange thing, though, is that when I just open tmp.tif, which has NOT been warped by gdal.Warp, in OpenEV (which seems to warp geotiffs on its own) I get a perfectly orthorectified image (though when I open tmp.tif in windows photo viewer it looks exactly like the source image):

Am I missing something obvious here?

EDIT: I suppose part of the reason for my question is that I don't understand on a deep level what gdal.Translate does. I know that it makes the image file into a GeoTiff based on the ground control points that I provide, but how is the location information stored in the GeoTiff? Is it just stored as a bunch of ground control points, or is it stored as a homography matrix along with a coordinate which could then be used to reproject the image to the correct perspective?

Whichever it is, I am assuming that both gdalwarp and OpenEV use this information to warp the image. What I find strange is that OpenEV does it perfectly while gdalwarp does not.


texts: How can I create maps in ArcGIS with a Baidu basemap?

texts:
I am working on a project mapping select China cities with some automotive related data. The client has requested we show Baidu as the basemap since it is more up-to-date than the OSM basemap we have been using. I have been able to bring in the Baidu basemap into ArcGIS, but I have reason to believe it isn't being displayed properly (at the correct location in relation to a WGS-84 basemap).
Can anyone shed light on this? Perhaps I need to make an adjustment once I bring in this Baidu vector layer to have it render in the right location?
Link to add Baidu basemap:
https://developpaper.com/arcmap-and-arcgis-pro-load-baidu-map/
The first image shows the WGS lat/lon of Tiananmen Square with a green dot and the Baidu lat/lon of Tiananmen Square with a pink square. It is set on the OSM basemap.

The second image shows the same data points, zoomed out. This will provide context for the third image.

The third image shows the same data points, now set on the Baidu basemap that I brought in. You can see how far off it is compared to the OSM basemap, and how the pink square does not align with it.


texts: QGIS shapefiles have the same CRS but don't line up. (Census TIGER Data & State Data)

texts:
I have the Census TIGER data for Oregon from

http://www.census.gov/geo/maps-data/data/tiger-data.html

Then the bottom one of 2010 Census Population & Housing Unit Counts — Blocks.

The shapefile loads fine into QGIS 1.8 and shows it as EPSG 4269/NAD83.

Then I have this city data from the State of Oregon Website.

http://navigator.state.or.us/sdl/data/shapefile/m2/cities.zip

It claims to be in EPSG 4269/NAD83 also.

When I first loaded it up the scale of both were WAY off.  The projection of the Cities you had to zoom out so far to see that the state population map became a tiny dot. (i.imgur.com/QZLQjfO.jpg)  I re-saved the cities shapefile and chose to save it as a EPSG 4269/NAD83 and reloaded it.  Now they both show on the same scale but the cities is far south of the Tiger file. (i.imgur.com/wI1OzSG.jpg)  I tried enabling on the fly projections to see if that would fix it and no go.  I manually loaded up some geocoded addresses under the NAD83 and they fit perfectly onto the TIGER/Census map so I think there is something wacky with the city data but I can't figure it out.

I am pretty new to QGIS and very much a novice at GIS in general.  Any tips or ideas on ways to get these to line up would be helpful.

I have found that if I save other files to 3644 they match up with the Cities data.  But if I save the cities data as any style let it be 3644 or 4269 it no longer matches up it becomes like Image 2 where it is way south of the other files.  This leads me to think that the Cities file has some sort of lat/long correction that when saving to other CRS or when re-projecting as other CRS causes it to show up very far south.  Now investigating how to remove/fix this attribute any and all help is appreciated.


texts: ArcGIS Time Slider - show hourly data only

texts:
I have point data that shows incidents in the years from 2011 to 2015 and it includes the specific time an incident occurred.

I would like to use the time slider functionality in ArcGIS Online to show the incidents that occur in each hour over a period of 24 hours.

For example: I would like it to show all incidents that occurred between years 2011 - 2015 between 13:00 - 14:00, 14:00 - 15:00, 15:00 - 16:00 etc.

How do I achieve this?


texts: How to change "metres" from "degrees" without change of location?

texts: I am working in QGIS (version 3.10). I know that in order to make certain distance analysis (in my case I want to do a heat map) a layer should have meters as units. Initially from the automatic SRC EPSG: 4326- WGS 84 I change my point layer' SRC to one set on metres (3857- WGS84/ Pseudo-Mercator). It does work in the sense the units do change to meters, however the points appear somewhere else. I tried to change the SRC of the project, start from scratch with the new SRC, different SRC combinations (using different SRCs for the project and/ or point layer) but just cannot get the units and the location both right.


texts: Convert shapefile to SVG while keeping attributes

texts:
Is there a way to export shapefile to SVG while keeping attributes?
I tried the mapshaper and QGIS' Mapprinter plugin but with no result, I get the geometry only but I need attributes.
I want when I export to fill the svg id with values from a column so I can join them afterward


texts: Adding hyperlinks to pop-up data in qgis2web?

texts:
Does anyone know if it is possible to add hyperlinks to image fields etc. in the pop-up windows of Qgis2web?

I've had a look in the layers javascript file, but the image source path isn't rendered using html, so I can't simply append <a href" to it.

The image is coming in via the usual method of having a field that contains the path of the image file.

I'd like to achieve linking of the image, like in the example below. Alternatively, I could just populate the 'more info' field with the URL, but if I did that I'd like to avoid just doing it by having a field that contains the entire URL, as some may be quite long (these will only be internal links).


texts: Rinex Header GPS Position Coordinate System

texts:
I am extracting data from rinex observation files and would like to convert the approximate position to latitude and longitude using Proj.4. The rinex standard states that these coordinates are WGS84, which I assume refers to the ellipsoid. What are the Proj4 parameters I need to do the conversion? If this is not possible is there another route I can take? Specifically I'm using the pyproj Python library. The correct coordinates should be close to (45.4343,-123.5516).

I am using RTKLib to process the observations, but want to automate downloading the base station/CORS files.

I did just find this code that performs the conversion. I've ported it to Python and it works well for my needs. Still if there is a way to get Proj4 to do this that would be ideal as I use that library for other projects.

Here is the Python I've tested:

In [29]: prj = pyproj.Proj(proj='geocent',ellps='WGS84',datum='WGS84')

In [30]: prj(-2478053.2006,-3736604.9657)
Out[30]: (1e+30, 1e+30)

In [31]: prj(-2478053.2006,-3736604.9657,inverse=True)
Out[31]: (-22.260730649602785, -33.5664935140521)

Here is an example header.

     2.00           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
                                        20180413 230424 UTC PGM / RUN BY / DATE
28jxmqr6                                                    MARKER NAME
                                                            MARKER NUMBER
                    off                                     OBSERVER / AGENCY
100037              BE-GPS-3300         2.1.47              REC # / TYPE / VERS
                                                            ANT # / TYPE
 -2478053.2006 -3736604.9657  4521684.4805                  APPROX POSITION XYZ
        0.0000        0.0000        0.0000                  ANTENNA: DELTA H/E/N
     1     1                                                WAVELENGTH FACT L1/2
     4    C1    L1    C2    L2                              # / TYPES OF OBSERV
  2018     4    13    19     8   50.5000000     GPS         TIME OF FIRST OBS
  2018     4    13    19    23   28.0000000     GPS         TIME OF LAST OBS
                                                            END OF HEADER


texts: qgis.org/ubuntugis-nightly KEYEXPIRED while updating QGIS3?

texts:
I'm trying to update QGIS 3 (2.99) today using

sudo apt-get update

But I get an error complaining about the keys expired. The Error message is:

Err:8 http://qgis.org/ubuntugis-nightly xenial InRelease
  The following signatures were invalid: KEYEXPIRED 1502999973  KEYEXPIRED 1502999973  KEYEXPIRED 1502999973

...

W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: http://qgis.org/ubuntugis-nightly xenial InRelease: The following signatures were invalid: KEYEXPIRED 1502999973  KEYEXPIRED 1502999973  KEYEXPIRED 1502999973
W: Failed to fetch http://qgis.org/ubuntugis-nightly/dists/xenial/InRelease  The following signatures were invalid: KEYEXPIRED 1502999973  KEYEXPIRED 1502999973  KEYEXPIRED 1502999973
W: Some index files failed to download. They have been ignored, or old ones used instead.

I tried to add-key again from the server as follows:

sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys 089EBE08314DF160
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys 073D307A618E5811

But the error remains.

Any suggestion on how to fix the error?

Also, if really the keys expired as suggested, do I need to remove them?


texts:
I have generated a stream from my DEM. I also have a shapefile of the river which i am working on. I want to overlay the river on the generated stream from the DEM and check if i could see whether the river overlayed will fall on the streams generated and then could generated watershed.

But the challenge is that, since most of the streams generated are not connecting to the path where the river shapefile is overlayed, i am finding it difficult to an entire watershed for the river. Only subwatersheds can i generate. I want to find a way of connecting the gap/space between two streams which are capable of forming subwatershed where the river passes.

So i just want to know whether such procedure is feasible and if so how could it be done.

I have attached a picture of the generated stream and the river shapefile (red line) for attention


texts: How can I add a graticule similar to the DeLorme Atlas to a map in QGIS?

texts: I've got a map in NJ State Plane (Feet) projection and want to add a grid to it with lat and long intervals of 1 minute, and 0.1 minute tick marks. I'm trying to do this in the grid tab in the print composer but getting wonky results (probably some user error here!


texts: Applying equivalent of image.expression function to imageCollection using Google Earth Engine

texts:
I'm just starting to work with Google Earth Engine and I'm trying to calculate a simple index based on an expression. I found this image.expression function, but can't figure out what is the equivalent for an imageCollection.
var s2 = ee.ImageCollection("COPERNICUS/S2_SR")
    .filterDate('2020-01-01', '2021-01-01')
    .filterBounds(geometry)
    .filter(ee.Filter.lt('CLOUD_COVERAGE_ASSESSMENT', 10));

var image = ee.Image(s2.first())

// **HOW DO I DO THIS TO AN IMAGE COLLECTION, INSTEAD OF JUST THE SINGLE IMAGE**?
var index= image.expression(
    '(1-((B6*B7*B8A)/B4)**0.5)*((B12-B8A)/((B12+B8A)**0.5)+1)', {
      'B4': image.select('B4'),
      'B6': image.select('B6'),
      'B7': image.select('B7'),
      'B8A': image.select('B8A'),
      'B12': image.select('B12')
});

Map.addLayer(index, {min: 0, max: 0.99, palette: ['FF0000', '00FF00']});


texts: Calculate longest distance, in KM, between county centroids and their furthest point with PostGIS

texts:
Using this shapefile of the US counties at 1:500,000 resolution, I need to know, per county, the length from the centroid of that county to its furthest point away on its polygon. I would like the result in kilometers.

I have uploaded the shapefile to a PostGIS table with:

shp2pgsql -W LATIN1 gz_2010_us_050_00_500k.shp public.county_shapes | psql -h localhost -d counties -U jerickson

I did not know which SRID to use so I did not specify one.

Then, modifying this question/answer, I tried the following query, limiting to name = 'Missoula (a county in Montana that I am familiar with) for now for testing:

SELECT
       t.gid AS county_id,
       max(ST_Distance(dump.geom, ST_Centroid(t.geom))) AS distance
FROM county_shapes t
JOIN ST_DumpPoints(t.geom) dump ON true
where name = 'Missoula'
GROUP BY 1
ORDER BY 1;

While this seems to be correct (logically), I am not sure how to get the resulting value (0.9038766878461234) into kilometers. Using Google Earth to approximate the value, I would expect it to be at least 75 km.

How can I modify my process to get the distance from the centroid to the farthest point on the county polygon for each county in the shapefile?


texts: Download salinity data from Government of Canada

texts:
Since I would like to use the salinity data with lat/lon in the Gulf of St. Lawrence, Canada to compare with the fishing data. I'm trying to download the 2010-2019 average salinity data. As the table in this database and also the picture showed. (https://www.arcgis.com/home/webmap/viewer.html?url=https%3A%2F%2Fgisp.dfo-mpo.gc.ca%2Farcgis%2Frest%2Fservices%2FFGP%2FTeleost_BottomTemperatureSalinityFond%2FMapServer&source=sd)

And here is the open data portal. https://open.canada.ca/data/en/dataset/10987662-c496-4ba8-a6b9-21cb5a134da2
Since my analysis requirement (I'll use R software), I want to download the data and mix it with my catch data under respective lat/lon. But I can't find how to get the raw data table in excel/csv format.


texts: Filtering Layers based on types in ArcGIS

texts:
I want to loop through all the layers in TOC in ArcGIS and report Raster and Vector Layers. I want to see how many layer are raster and how many are vector. I can access the layers using:

IEnumLayer pEnumLayer = pMap.get_Layers(null, true);

and loop the layers in Enum, but how do I determine which layer is Raster.


texts:
I'm trying to load a layer from a PostGIS database using PyQGIS.
The code is:

import os
from qgis.core import *

project = QgsProject.instance()
project.clear()

uri = QgsDataSourceUri()
uri.setConnection("localhost", "5432", "db_name", "user", "pass")
uri.setDataSource("public", "table_name", "the_geom", "")
vlayer = QgsVectorLayer(uri.uri(False), "table_name", "postgres")
if not vlayer.isValid():
    print("Layer failed to load!")
else:
    project.addMapLayer(vlayer)

The problem is that I don't know what to write in "the_geom" (I tried to insert "point" or something similar).

Is there someone who can tell me?


texts:
I've installed QGIS-OSGeo4W-2.14.11-1-Setup-x86_64.exe on a windows 7 PC
the installation completes without error, we restarted the PC as requested.
We double clicked the below shortcut:
"GRASS GIS 7.2.0" the shortcut points to "C:\Program Files\QGIS 2.14\bin\grass72.bat" -gui but this file "grass72.bat" does not exist.

Where can I obtain it?


texts: How do I use Proj4js or other technique to make custom map projections with google.maps.Projection?

texts:
I'm using Google's Projection class to try to convert my map projection from the default spherical mercator to something like Lambert Conformal Conic. I was previously adapting the code from this example to the v3 API, but am getting some weird drawing errors. Then I remembered Proj4 and decided to try it out.

I'm not sure how to use it with the Projection's two main functions:

fromLatLngToPoint()
fromPointToLatLng()

I understand how Proj4 returns a meter value when you pass it a lat long and a destination projection, but I'm not sure how that could be used with the Projection class's functions (or if this is even the correct approach to doing custom projections in Google maps anyway).

I think the main thing that I'm not sure about is how to use the Point object in these functions. Am I going to need to take into account my map canvas size + my map bounds to do a calculation on my Point's relative position within the frame?

Here's an image of what it looks like, and my code is below:

// my base layer, currently with no tiles
var lambertBase = new google.maps.ImageMapType({
    getTileUrl: function(coord, zoom) {
        return null;
    },
    tileSize: new google.maps.Size(256,256),
    minZoom:  1,
    maxZoom: 8,
    name: 'Lambert'
});

lambertBase.projection = new ConformalConicLambertProjection();

var myLatlng = new google.maps.LatLng(34, -96);

// map options
var mapOptions = {
    zoom: 4,
    maxZoom: 8,
    center: myLatlng,
    panControl: false,
    streetViewControl: false,
    mapTypeControl: false,
    zoomControlOptions: {
        style: google.maps.ZoomControlStyle.SMALL,
        position: google.maps.ControlPosition.TOP_RIGHT
    }
}

// Draw the map into the div with the above options
var vMap = new google.maps.Map(document.getElementById('mapCanvas'), mapOptions);

vMap.mapTypes.set('lambertTiles', lambertBase);
vMap.setMapTypeId('lambertTiles');

// then I draw the county polygons

// projection utils (modified from link above)
var phi0 = 10*Math.PI/180;
var phi1 = 20*Math.PI/180;
var phi2 = 70*Math.PI/180;
var n = (Math.log(Math.cos(phi1)) - Math.log(Math.cos(phi2)))/
    (-Math.log(Math.tan(phi1/2 + Math.PI/4)) + Math.log(Math.tan(phi2/2 + Math.PI/4)));
var rho0 = (Math.cos(phi1)*Math.pow(Math.tan(phi1/2 + Math.PI/4),n))/
    (n*Math.pow(Math.tan(phi0/2 + Math.PI/4),n));

function ConformalConicLambertProjection() {};
ConformalConicLambertProjection.prototype.fromLatLngToPoint = function(latLng) {
    var z = vMap.getZoom();
    var t = Math.PI/180*(100 + latLng.lng())*n;
    var rho = (Math.cos(phi1)*Math.pow(Math.tan(phi1/2. + Math.PI/4.),n))/
        (n*Math.pow(Math.tan((Math.PI/180*latLng.lat())/2. + Math.PI/4.),n));
    var x1 = rho*Math.sin(t);
    var y1 = rho0 - rho*Math.cos(t);
    x = Math.round(Math.pow(2,z)*(128 + 290.9090909090909*x1));
    y = Math.round(Math.pow(2,-1 + z)*(405.68181818181824 - 579.5454545454546*y1));

    return new google.maps.Point(x, y);
};
ConformalConicLambertProjection.prototype.fromPointToLatLng = function(pixel, noWrap) {
    var z = vMap.getZoom();
    var x1 = (-0.0034375*(Math.pow(2,7 + z) - pixel.x))/Math.pow(2,z);
    var y1 = -0.0017254901960784311*Math.pow(2,1 - z)*
        (pixel.y - 405.68181818181824*Math.pow(2,-1 + z));
    var rho = Math.sqrt(2.143305952697287 + Math.pow(x1,2) -
        2.9280067982826044*y1 + Math.pow(y1,2));
    var t = Math.asin(x1/rho);
    var lng = 78.17779218926289*t-100;
    var lat = 28.64788975654*(-3.141592653589793 +
        4*Math.atan(2.004758218840869/Math.pow(rho,1.3644598756425434575)));

    return new google.maps.LatLng(lat, lng, noWrap);
};


texts: Create light overlay from heavy shapefile to use in web app

texts:
A question about soundness of principle as much as method, as I am very new to the subject.

Having to serve over web the polygons in way too large a shapefile (200Mb), I am wondering: could I, with some tool, convert all the data in the shapefile into a single transparent image (png I suppose, covering about 20,000km2) to then use as an easy to load overlay in google maps ?

Edit:
I seem to understand Geoserver and Mapserver would create and serve selected, appropriate tiles of the resulting image on demand.
Not lastly because we don't have a server to run Geoserver at hand, and since the data is static, I was thinking more of an actual image file. Does it make sense ?


texts: Return corner coordinates of images after OpenDroneMap processing

texts:
I have images (with EXIF metadata) made by a drone relative to a particular place; having processed them with OpenDroneMap at this point I have an orthophoto with all the geolocated pixels.
At this point I was wondering if through Python and given the output of OpenDroneMap (GeoTIFF or other), it was possible to trace the coordinates of the angles of the individual photos used for the creation of this output.
Is it possible?


texts: How to zoom out and switch map in Oruxmaps

texts:
I'm using Oruxmaps.
My main maps are at zoom 16 (1:25000). Thing is I can zoom int from different offline maps (zoom 14 for instance) to zoom (my main map 1:25000 ), but when I want to go back to zoom 14, it is not possible, Oruxmaps won't let me.
In my case, I will hike this summer all Pyrenees. So, I will frequently have to move east to west 900km. At zoom 18, I must be swiping 5 minutes with my fingers, and I sometimes lose where I currently am, it is not practical at all.
What I want to achieve is that if I want to move to a distant place is to be able to zoom out to another map zoom 14 or less, move quickly to the other point, and then zoom in again.
To conclude switching offline maps is working when zooming in, but not when zooming out.
Is there a configuration or a trick to achieve this?


texts: Is ArcGIS Pro truncating of some legend titles documented?

texts:
To understand what was causing a legend title to be truncated from "What's on the map?" to "What's on the" in an ArcGIS Pro 2.1.3 layout I developed the test below:

Start ArcGIS Pro
Create a new project using Map.aptx called TestProject in a new folder under C:\Temp
Add one layer to the map - in my test I used a countries shapefile (polygons) from Natural Earth but any shapefile or feature class should show this I think
No matter what your layer is called name it ne_10m_admin_0_countries (or something with more characters in its name)
Insert a Layout - I chose an A4 Portrait but any should be fine
Insert a Map Frame - I chose the only Map and the one with a scale rather than the Default but either should be fine
Make some room on the page for a legend under the map
Insert a Legend and position it by dragging out its rectangle under the map
Save the project and close ArcGIS Pro
Use IDLE (or your preferred Python IDE) to run the script below

import arcpy,os
arcpy.env.overwriteOutput = True
aprx = arcpy.mp.ArcGISProject(r"C:\Temp\TestProject\TestProject.aprx")
lyt = aprx.listLayouts()[0]
legend = lyt.listElements('LEGEND_ELEMENT',"legend")[0]
legend.title = "What's on the map?"
legend.showTitle = True

lyt.exportToPDF(r"C:\Temp\test.pdf")
os.startfile(r"C:\Temp\test.pdf")

On the PDF that opens the legend title of "What's on the map?" appears in full
Close Adobe Acrobat Reader and the Python Shell window of IDLE to release any locks
Open the TestProject on ArcGIS Pro again
Rename the layer to Countries
Save the project and close ArcGIS Pro
Use IDLE (or your preferred Python IDE) to run the script again

This time when the PDF opens the legend title is truncated to "What's on the"

It appears that the length of the longest line of text next to the legend patches can restrict how many characters can appear the legend title.

Is this software behaviour documented anywhere?


texts: Using OpenStreetMap for routing in GeoServer?

texts:
I am new in OSM but now I need create application which will use GIS (Geoserver) and will be provide service to get road from point A to point B.

I found some informations about this but I still have any questions. My is download open street maps and extract to database (PostgreSQL), and for routing use plugin pgRouting. Next connect db to map server and provide api for clients.

Is that correct way?

How osm format could I use for postgre? (osm.pbf? shp.zip?, or something else?)


texts:
Is there any information within a map service regarding the name of the source mxd that was original used to publish/share the data?  Assuming the user shared the data to ArcGIS Online from within ArcMap.

I'd ideally like to use Python and query our services to catalog the various map documents that might have been used to publish/share the service.  This is in support of a data "clean-up" and workflow optimization effort.


texts:
I've been looking through these threads, i found something similar and thought it could be the answer, but it doen't seem to work for me...

I could do this the slow way with Select By Attributes, but It's time consuming and I will have to repeat it several times for other works.

I have a field called "Kultur" (String) with different codes (DG, A, FF, ...). I want to change this information into a numeric information, so i added a new field called "Kultur_Cl" (Short Integer).

In my case if "Kultur" is DG, "Kultur_Cl" should be 1. If "Kultur" is A, "Kultur_Cl" should be 2, ...

Here is what i tried :

And here is what it looks like in ArcMap:

Once it's done i would like to do : If "Kultur_Cl" = 1 then "Kultur_Pkt" = 15, if "Kultur_Cl" = 2 then "Kultur_Pkt" = 18, ... Would this go the same way, or is it different because there would only be numeric data?

Thanks for your help !!


texts: Can I create a selection on Feature Table by attribute with ArcGIS JS API?

texts:
I have a simple app built using the 3.16 API that includes a flat table (so no spatial attributes) using the Feature Table widget. Can I select rows by attribute rather than just by the user clicking on a row? For example, if I wanted the table to display only rows where the "State" field equals 'California' when its initially loaded, but the user could choose to clear the selection to see all rows if desired? Or would I be better off setting a query on the feature layer that the table is built from?

My code is below.

require([
    "esri/map",
    "esri/urlUtils",
    "esri/arcgis/utils",
    "esri/layers/ArcGISDynamicMapServiceLayer",
    "esri/layers/FeatureLayer",
    "esri/layers/LOD",
    "esri/layers/LayerInfo",
    "esri/dijit/FeatureTable",

    "dojo/dom",
    "dojo/_base/array",
    "dojo/ready",
    "dojo/parser",
    "dojo/on",
    "dojo/domReady!",

    "esri/IdentityManager"
],

function (Map, urlUtils, arcgisUtils, ArcGISDynamicMapServiceLayer, FeatureLayer, LOD, LayerInfo,
      FeatureTable, dom, rrayUtils, ready, parser, on) {

    // Wait until DOM is ready *and* all outstanding require() calls have been resolved
    ready(function () {

    // Parse DOM nodes decorated with the data-dojo-type attribute
    parser.parse();

    //create a regular plain old map
    var mapMain = new Map("mapDiv", {
        basemap: "topo",
        center:  [-117.4, 33.3], //[-120, 38],
        zoom: 11,
        showLabels: true
    });

    mapMain.on("load", function() {loadTable(mapMain); });

    function loadTable(myMap){
    console.log("Here comes a table.")

    // Create the feature layer
    myFeatureLayer = new FeatureLayer("https://data.farwestern.com/arcgis/rest/services/scratch/test_SD_RelatedResources/MapServer/2", {
        mode: FeatureLayer.MODE_ONDEMAND,
        visible: true,
        outFields: ["*"]
    });

    mapMain.addLayer(myFeatureLayer)

    // Create the feature table
    myFeatureTable = new FeatureTable({
        featureLayer : myFeatureLayer,
        outFields: ["State", "County", "Name"],
        map : myMap,
        gridOptions: {
        allowSelectAll: true,
        allowTextSelection: true,
        selectionMode: "extended"
        },
        syncSelection : false
    }, 'tableDiv');
    myFeatureTable.startup();
        }

  });
});


texts: How to programatically add KML Icons to Placemark generated from SimpleFeatures using GeoTools?

texts:
I have successfully generated a set of (unstyled) placemarks in a KML file using the GeoTools framework, and the Encoder class. The file renders fine in Google Earth with standard balloon placemarks.

I now want to add some iconstyles to my file. Can anyone give me some guidance here? What is the best way to programatically declare the icons/iconstyles and give each placemark a reference to the corresponding icons? I know the KML syntax for this, but how do I declare/generate this with GeoTools, and make the Encoder class do the final job?


texts:
I used to gdalwarp to reproject raster tiff file, but I did it one by one

How can I manage to do all raster in my folder (about 70 file) ? Is it possible? I am using ubuntu-64bit


texts: Maps, black and white/ gray scale printing

texts:
I'm working on Mapserver to do beautiful maps with AGG rendering

Until now, I've only thought about the best way to get color maps. Now, I have to think how to make color ones compatible with a grey scale printing
http://en.wikipedia.org/wiki/Grayscale

How do you deal with this?

Any advices, tips or resources will be welcome.

Edit : My question is really more relative to visual perception than Mapserver.
It seems I will have to read Semiology of Graphics by Jacques Bertin (reedited this year http://www.esri.com/news/releases/10_4qtr/bertin.html ) or look on Edward Tufte various book/articles http://www.edwardtufte.com/tufte/


texts:
I am currently preparing raster layers in QGIS to be fed into MaxEnt to run a species distribution model.

I downloaded climatic variables in ascii format from I downloaded from Adaptwest (https://adaptwest.databasin.org/pages/adaptwest-climatena) with the following parameters:

Projection    LAMBERT
Datum         WGS84
Spheroid      WGS84
Units         METERS
Zunits        NO
Xshift        0.0
Yshift        0.0

Parameters
49  0  0.0 /* 1st standard parallel
77  0  0.0 /* 2nd standard parallel
-95  0  0.0 /* central meridian
0  0  0.0 /* latitude of projection's origin
0.0 /* false easting (meters)
0.0 /* false northing (meters)

I will need to crop this raster with states boundaries shape files in CSR WSG 84 but I am unable to successfully re-project the ASCII files so that they overlap with the states shape files.

I tried a solution posted in this Q&A site in which I created a CSR using the above parameters:

+proj=lcc +lat_1=49 +lat_2=77 +lat_0= 0 +lon_0=-95 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

but I am still unable to re-project. I get no error messages but the outcome is wrong.

Any idea what am I doing wrong?


texts: Using ArcGIS ModelBuilder in to Select By Attributes?

texts:
I have a file showing for each disaster, which provinces were affected. The number of provinces affected varies and their names or associated code varies as well.

Example, if a flood in 2000 affected 3 provinces called Sichuan, Henan and Yunnan. I need to select those 3 provinces and create a new shapefile. Easy. But I need to do that for over 2000 disasters and the affected provinces always changes (in names and numbers of provinces affected).

I would like to automate this with ModelBuilder but I don't know how to change my SQL request for each disaster.

Do you know if there is a way for ModelBuilder to see in an Excel file, which are the provinces that need to be selected and select them and create the corresponding shapefile?


texts: How to pick adjacent GeoTIFF files out of a set?

texts:
I have a large amount of large GeoTIFF images. Their names are arbitrary and the person that created them does not work here with us anymore.

I need to do some programming experiment. How can I pick a small subset that covers adjacent ground?

Edit:

I have been tasked with creating a prototype of a slippy map with these GeoTIFF files; they are high resolution aerial imagery of the town we're based in. I want to experiment with GeoServer, as advised in another question here, but the files take so much space that I'd like to start with a smaller subset. So I would like to pick a small number of these GeoTIFFs, and to make a sensible slippy map they should be adjacent - covering areas of ground that are next to each other.


texts: Error distance for coordinate transformation (to WGS84)

texts:
I have a dataset with mostly EPSG:31255 and wanted to convert it to WGS84.
I now had to realize that the results are up to 100m off when I plot them on a map or just enter the resulting lat lon in google maps.

Is this within the expected error margin? I use the data for geocoding. So anything > 10m is actually not acceptable

I tried it with pyproj and also with the following tool (which produces different coordinates): cs2cs.mygeodata.eu

This is a sample:

Plainstraße 1, 5020 Salzburg, Austria
input:
EPSG 31255
296772.02
-21921.52
output:
lat: 47.80971875756464
lon: 13.040617960746806


texts: Open source alternative to the "Contours with Barriers" Tool in ArcGIS

texts:
I am trying to create a contour file for the globe for a project and am using open source GIS software packages. I was trained for several years on ArcGIS and there is a tool to create contours with barriers allowing the user to provide a vector file of "barriers" and create contours with these in mind. In my case I am looking for the open source alternative using QGIS, SAGA, GvSIG, GRASS. I have tried creating a mask in grass and applying the contour function to my DEM without any successful results.

here is a link to the tool description:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z000000tv000000.htm

This is my first post, please forgive me if I have left any vital information out.


texts: Custom on-the-fly conversion into shapefiles

texts:
I have a custom binary format, for which I run an external program to output a set of .shp files, which then can be dragged and dropped into QGIS.

What would be the best way of pushing this into some QGIS hook that would do all of this behind the scenes to support something like this:

Custom binary file is dragged into QGIS
QGIS calls the external program, which outputs the .shp files into some temporary folder.
QGIS shows the temporary folder

Browsing through the plugin docs of QGIS 3.4 I didn't find a clear guidance wrt. where I could intercept the drag'n'drop / file-loading process. Also given the excellent Python-extensibility of QGIS, I suspect that there is already some support for my use case.


texts: Reverse geocode through Google Maps API inaccurate?

texts:
I have placed around a 1000 markers in my area through Google My Maps, transformed the KML file to CSV, and now want to transform the coordinates to addresses. I am just testing with the first 5 to test my script (Python) but I am finding that from the 5 test records, only 1 is correct.

For example:
Coordinate: -33.7909042,151.2603029
Google map: http://maps.google.com/maps?force=canvas&t=m&q=loc:-33.7909042+151.2603029

Clearly you can see the address is Rickard St 42

BUT
When you use the API to geocode the coordinates with:
https://maps.googleapis.com/maps/api/geocode/json?latlng=-33.7909042,151.2603029&key=APIKEY

It returns 119 West Street, on the opposite side.

I thought the Google Maps API was really accurate, but the service doesn't reflect the actual map itself. What is going wrong? Is there a better service?


texts: Several algorithms done for topographic correction => black spots - how to solve?

texts:
I have done several algorithms for topographic correction in R including simple cosine, Minnaert, SCS. The input data aspect and slope for the hillshade raster have correctly been done in radians. The hillshade has some spots where the cos (i) is slightly below zero and these spots appear black in the corrected Landsat image. The white surroundings of the black spots are zero or slightly above. That means that cos(i) is around 90° or above in these areas. Such areas receive no direct radiation from the sun but diffuse light.

Can somebody help me with this - I mean do you know an algorithm to solve this problem? Is this phenomenon called overcorrection?

Thanks a lot
S.


texts:
How can I overwrite the url of OpenLayers.Protocol.HTTP? I tried

  searchformPanel.protocol.url

and it works (checked with console.log), but at the end the original url is sent (in the code below: url: '/fs/') (see attached picture).

This is the code:

var searchformPanel = new GeoExt.form.FormPanel({
  border: false,
  width: 250,
  protocol: new OpenLayers.Protocol.HTTP({
    url: '/fs/',
    format: new OpenLayers.Format.GeoJSON()
  }),
  items:[{
    xtype: 'combo',
    id: 'idcombo',
    store: new Ext.data.SimpleStore({
      fields:['fsclass','ollayer'],
      data:[["Boreholes","Boreholes"],["All_layers","All layers"]]
    }),
    displayField: 'ollayer',
    valueField: 'fsclass',
    fieldLabel: 'Layer',
    emptyText: 'select a layer',
    submitValue: false,
    selectOnFocus: true,
    mode: 'local',
    typeAhead: true,
    editable: false,
    forceSelection: true,
    triggerAction: 'all'
  },{
    xtype: 'textfield',
    id: 'idtextfield',
    fieldLabel: 'Find features',
    emptyText: 'enter word',
    name: 'comments__like',
    allowBlank: false
  }],
  listeners:{
    actioncomplete: function(form, action){
      searchShowTip(action.response.features);
    }
  },
  buttons:[{
    text: 'search',
    listeners:{
      click: function(){

        var comboLayer = Ext.getCmp('idcombo').getRawValue();
        var keyword = Ext.getCmp('idtextfield').getRawValue();

        var newUrl = '/fs/' + comboLayer + '?format=GeoJSON&comments__ilike=' + keyword + '&queryable=comments';
        console.log('1:' + newUrl);

        //this gets '/fs/' from the searchformPanel
        console.log('2:' + searchformPanel.protocol.url);

        searchformPanel.protocol.url = newUrl;
        console.log('3:' + searchformPanel.protocol.url);

        searchformPanel.search();
      }
    }
  }]

});

Please any support about this is very welcome, thank you!


texts:
I have found that the map section shown on one viewing is not the same as shown on a second viewing.  The same designated scale produces different dimensions in the second viewing and places that are named on one viewing are not shown on the second viewing.

Why can't Google Maps be relied on to show the same information on multiple viewings?


texts: Is there flood map prediction tutorial available to use with ArcGIS for Desktop?

texts:
I just want to ask if anyone knows a direct tutorial for flood mapping prediction.
The scenario that my supervisor asked me is this:
If we pour water on a specific location, where will it flow and how much is the runoff to cause a flood?


texts: Calculate area above elevation threshold?

texts:
I have a number of 95% isopleths (polygons) that represent animal home ranges.  For each isopleth, I and want to calculate the Area (sq km) above 2000 meters.

My initial thought was to create contours from regional DEMs but because the lines do not close to form polygons, I am unsure how to calculate area.

An ArcGIS solution is preferred for this project.


texts: Using R NetCDF raster stack is missing its CRS, how can I assign one?

texts:
The data I am using was not generated by me, however, I am stuck with using it anyways. My issue is, that in R, the CRS of the raster stack is missing or NA. If I open it with ArcGIS or Unidata (!the file is a NetCDF!) it gives me an odd description I am not aware of. ('char lambert_conformal_conic_2sp')
Overall within the NetCDF are 2 different layers, storing (SWI_001 and SWI_010) pixel values I want to extract for every single day of the year.
I tried to check this webpage for further information about the 'char lambert_conformal_conic_2sp' projection, but it led to nowhere (https://mygeodata.cloud/cs2cs/).
Of course I also tried several different coding approaches, in order to find a way to assign a CRS to the file, similar to:
proj4string(r)=CRS("+proj=aeqd +lat_0=53 +lon_0=24 +x_0=5837287.81977 +y_0=2121415.69617 +ellps=WGS84 +units=m +no_defs ")

f <- projectRaster (r, '+proj=longlat + datum=WGS84')

However, I also tried to assign the same CRS of my file to the poylogon I am trying to extract the final data with, but since the CRS of the NetCDF file is NA to begin with, it doesnt work.
fieldnew <- sf::st_transform(field, crs = "***") #same projection as raster stack

The code I am working with at the moment is:
library(raster)
library(sf)
library(rgdal)
library(exactextractr)
library(ncdf4)
library(ggplot2)
library(reshape2)

files <- list.files('D:\\***\\***\\ASCAT_2018',pattern='*.nc',full.names=TRUE)

r <- stack(files, varname='SWI_001')

field <- read_sf(dsn = 'D:\\***\\***\\Studyarea_Polygon\\Polygon for R', layer = "testfield")

##plot extend of one file
plot(r[[1]], ext=10*extent(fieldnew)) ## check if polygon aligns with raster and if values make sense
plot(fieldnew, add=TRUE)

## extract by weights of polygon Raster,SpatialPoylgon and relevant observations
output <- extract(r [[1:197]], fieldnew, fun=mean, weigths=TRUE, normalizeWeights=FALSE, layer=1, method='bilinear', df=TRUE)

The error messages are:
> f <- raster(files, varname='SWI_001')
Warning:
1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
long_name=CRS definition
spatial_ref=PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6377397.16,299.1527033239178]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",46],PARAMETER["standard_parallel_2",49],PARAMETER["latitude_of_origin",47.5],PARAMETER["central_meridian",13.333333],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]
GeoTransform=-380501 500 0 220501 0 -500
2: In .getCRSfromGridMap4(atts) : cannot create a valid CRS
grid_mapping_name; false_easting; false_northing; scale_factor_at_projection_origin; scale_factor_at_central_meridian; standard_parallel; standard_parallel1; standard_parallel2; longitude_of_central_meridian; longitude_of_projection_origin; latitude_of_projection_origin; straight_vertical_longitude_from_pole; longitude_of_prime_meridian; semi_major_axis; semi_minor_axis; inverse_flattening; earth_radius; +proj; +x_0; +y_0; +k_0; +k_0; +lat_1; +lat_1; +lat_2; +lon_0; +lon_0; +lat_0; +lon_0; +pm; +a; +b; +rf; +a

Is there a way to assign a CRS to my files, or is it a lost cause?
EDIT:
After futher research I found some entries about gdal_warp (gdalUtilis) and CDO (Changing a raster stack projection in R), is it possible to assign/transfrom the missing CRS based on these options?


texts: Proportional symbols for data with wide range

texts:
I have point shapefile with 15 points whose attribute for symbol ranges from 2000 to over a million. There are 3 points with over 1mil while the other 12 points are all in thousands. I want the size of circles to represent the values of this attribute. But when I use proportional symbols, the 3 points with the values greater than 1 million dominate and other 12 points are just tiny points.

So I am thinking I would set an upper limit that would still show these 3 points as big circles but I would still be able to see proportional sizes for other 12 points. How do I go about doing this? And is this the best way to represent these points proportionally?


texts: From where can I get Landsat 7 and Landsat 8 metadata files?

texts: I downloaded satellite images but mtl files were not in the zip. I searched a lot but i dint get from where to download mtl files. I need to convert radiance to temperature in both L7 and L8 images.


texts: Can the QGIS spatial index be used to evaluate points contained in multiple polygons of shapefile layers using PyQGis?

texts:
I want to use the spatial index of QGIS to improve the performance in the spatial query of multiple points in multiple polygons, both are shapefiles, but has not given me good results in execution time.
The execution time is slower than a sequential algorithm with a for loop using the QgsGeometry methods (even separating the spatial index creation time).
On the other hand, the result when querying, the amount of polygons contained in polygons, of two Shapefile layers, the result is as expected and the execution time is in tenths of a second.
Is this a limitation of the QGIS spatial index (QgsSpatialIndex) or am I implementing it wrong for the case of points? here is the code
import time

cliente=pry.mapLayersByName('cliente')[0] # points layer
parcelas=pry.mapLayersByName('parcela')[0] # polygon layer
campoP=parcelas.fields().field('id') #field to identify the polygon

#time to create spatial index
tiempo_inicialCIE = time.time()
index = QgsSpatialIndex() # Spatial index
index = QgsSpatialIndex(cliente.getFeatures())
duracionCIE=time.time()-tiempo_inicialCIE

dic3=[]    # list storing the result, id and number of points contained
lisIds=[]  # list of ids returned by the spatial index
tiempo_inicialIE = time.time()
for i in parcelas.getSelectedFeatures():
    bbx=i.geometry().boundingBox()
    for j in cliente.getFeatures():
        listIds=index.intersects(bbx)
    cv=0
    for c in listIds:
        if i.geometry().contains(cliente.getFeature(c).geometry()):# Evaluate intersecting features
            cv=cv+1
    if cv > 0:
        dic3.append([i[campoP.name()],cv])

duracionIE=time.time()-tiempo_inicialIE

Note: I have already optimized the process using other libraries and reduced the execution time for all entities from  +20 minutes to 4.5 seconds. But I would prefer to do it with PyQGis since it is for a PlugIn.


texts: Seeking GIS tools for simulating landscape fragmentation?

texts:
I am currently looking for any GIS packages (or add-ons) that allow the simulation of landscape fragmentation on raster data.

Maybe also other software than GIS tools?

Does anyone of you have had experience with that issue already and can suggest a good tool for that task?

Info on tutorials etc. are also welcome!


texts:
I have to calculate a NDVI raster with 2 Landsat maps.
The formula is as following :

The problem is that it doesn't show up in my layer list and when I try to export it from GRASS via r.out.gdal.gtiff, I get this:

My colleagues also get this message but they don't have the "Unable to create dataset using driver" part, and it works well on their computer.

I used both Qgis 2.12 Lyon and 2.8.3 Wien and it does the same thing.
Can you guys help me to fix this please ?


texts: Using ArcGIS ModelBuilder iterator to for Raster Calculator?

texts:
I have a raster that I need to divide by X in a model. I have hundreds of values I need to divide it by. What I am trying to do is create a list of values, and use the iterator to loop through them and divide my raster by that value with the Raster Calculator

Ex) If I want 10, 15, and 20, I will make a list of those three, start the model, and it will run three times.

How do I do this?

I've been messing around with the multivalue iterator with no success.

Here's what I have right now:


texts: How to calculate a new bounding box that is 80% of the actual bounding box returned by google map?

texts:
I need to show a marker for each user that falls within the visible portion of google map. To do so, I get the bounding box from google map, and send it down to the database and I receive a list of users that fall within the box. Then loop through the list of users and place a marker on the map for each user.

All works well, however, if the lat/lng falls just bellow the top border, then only the bottom of the marker is visible. If the lat/lng falls next to the side borders, then only 1/2 of the marker is visible.

One way to tackle this problem is to send a reduced version of the bounding box to the database. This way, you'd have a bounding box with margins, but you don't show anything within the margins.

Now, what is the best way of shrinking down the bounding box before sending it down to the database?

Is there any better way of doing this?

Check the image, 1/2 of the marker is hidden. (top)


texts:
I run ogr2ogr to transfer geometrical data from a pgsql (PostGIS) database to a MS SQL one. It works quite well. I'm trying to achieve exactly the same result in Python, as I'm planning to use programming to automate the entire process and add some error reporting to it, but I can't seem to wrap my head around it. I run Python 3 in an Anaconda environment in which I installed GDAL.

How do I run the following code in native Python/GDAL?

ogr2ogr -lco PRECISION=NO -lco SRID=27700 -lco SPATIAL_INDEX=NO -lco OVERWRITE=YES -lco UPLOAD_GEOM_FORMAT=wkt
-f "MSSQLSpatial" "MSSQL:server=server;database=db;trusted_connection=yes"
PG:"dbname=db host=host port=1234 user=user password=pass"
"schema.table" -t_srs "EPSG:27700" -s_srs "EPSG:27700" -a_srs "EPSG:27700" -gt 100000 -progress -nln schema.table

I can load the GDAL library in Python and there is an OGR component, but I can't find any ogr2ogr component.


texts:
I am trying to create a geotiff and the software (raster design) is asking for a code that I thought was EPSG code. however while trying to find out what epsg = us48 I noticed that the dialog actually says geotiff code. Now I am really in a loop.
autocad map contains a coordinate system that is named us48 albers equal area, nad 27 meter orig lat 23, orig long -96, northern standard par 45.5, southern standerd par 29.5
ESRI difines the same system as...

  'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",23.0],UNIT["Meter",1.0]]'>

and it says USGS after it.
This is the same as EPSG 6703,
I cannot find anything that matches in the raster design code.
I found USA US National Atlas Equal Area but it is lambert azimuthal 45, -100 the code for it in raster design is 2163
Apparently the autodesk version is supposed to be the same as epsg as desribed in this document...
http://svn.osgeo.org/fdo/trunk/Providers/GenericRdbms/com/ExtendedCoordSys.txt/

The values in this document are what I am using in the raster design dialog that I thought were epsg.
Can someone straighten me out?
I want an epsg code for US48


texts: Road network per grid cell with weighting

texts:
Im trying to calculate road network per grid cell.

I did it with the help of Measure length within vector grid? & Calculate line distances inside of polygons using ArcGIS for Desktop? but I am also trying to weight these roads by type.

For example, if road is main its impact should be high, but if road is secondary it will be lower. Is this possible?


texts: What causes ValueError: I/O operation on closed file?

texts:
I obtain the following error message:

ValueError: I/O operation on closed file

Here is my code (It selects each polygon contained in the shapefile "selected_polygons.shp" and builds the least-cost path between the previously selected polygon and each polygon contained in the same shapefile)

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:\Users\LCP"
arcpy.env.extent = "costsn.tif"

with arcpy.da.SearchCursor("selected_polygons.shp",("*")) as rowsInPatches_start:

for rowStart in rowsInPatches_start:

    ID_patch_start = int(rowStart[3])

    ## Define SQL expression for the fonction Select Layer By Attribute
    expressionForSelectInPatches_start = "GRIDCODE=%s" % (ID_patch_start)

    ## Process: Select Layer By Attribute
    arcpy.MakeFeatureLayer_management("selected_polygons.shp", "Selected_patch_start", expressionForSelectInPatches_start)

    ## Process: Cost Distance
    outCostDist=CostDistance("Selected_patch_start", "costsn.tif", "", "outCostLink.tif")

    ## Save the output
    outCostDist.save("outCostDist.tif")

    with arcpy.da.SearchCursor("selected_polygons.shp",("*")) as rowsInSelectedPatches_end:

        for rowEnd in rowsInSelectedPatches_end:

            ID_patch_end = int(rowEnd[3])

            ## Define SQL expression for the fonction Select Layer By Attribute
            expressionForSelectInPatches_end = "GRIDCODE=%s" % (ID_patch_end)

            ## Process: Select Layer By Attribute in Patches_end
            arcpy.MakeFeatureLayer_management("selected_polygons.shp", "Selected_patch_end", expressionForSelectInPatches_end)

            ## Process: Cost Path
            outCostPath = CostPath("Selected_patch_end", "outCostDist.tif", "outCostLink.tif", "EACH_ZONE","FID")

            ## Save the output
            outCostPath.save('P_' +  str(int(ID_patch_start)) + '_' + str(int(ID_patch_end)) + ".tif")

            ## Write in a file .txt
            outfile=open('P_' +  str(int(ID_patch_start)) + '_' + str(int(ID_patch_end)) + ".txt", "w")
            rowsTxt = arcpy.SearchCursor('P_' +  str(int(ID_patch_start)) + '_' + str(int(ID_patch_end)) + ".tif")
            for rowTxt in rowsTxt:
                value = rowTxt.getValue("Value")
                count = rowTxt.getValue("Count")
                pathcost = rowTxt.getValue("PATHCOST")
                startrow = rowTxt.getValue("STARTROW")
                startcol = rowTxt.getValue("STARTCOL")
                print value, count, pathcost, startrow, startcol
                outfile.write(str(value) + " " + str(count) + " " + str(pathcost) + " " + str(startrow) + " " + str(startcol) + "\n")
                outfile.close()


texts: Accessing differing geodatabase schemas from arcpy

texts:
I am attempting to write a script that processes a lot of different feature classes in a geodatabase.  Most of the feature classes have the same schema as the user connecting to the database

geodatabase.myschema.featureclass

however some feature classes have a schema belonging to a different user

geodatabase.differentschema.featureclass2

When using arcpy with a .sde connection and user that owns the default schema, I do not need to specify the schema name when referencing feature classes in arcpy tools.

import arcpy, os
db = r"D:\GIS\DBConns\geodatabase.sde"
fc = "featureclass"
fcpath = os.path.join(db, fc)

result = int(arcpy.GetCount_management(fcpath).getOutput(0))
print result

> 147

However if I reference the feature class with a different schema, the above doesn't work

Runtime error Traceback (most recent call last): File "<string">, line 6, in <module>
File "C:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 15370, in GetCount
raise e Execute Error:
Failed to execute. Parameters are not valid.
Error 000732: Input Rows:
D:\GIS\DBConns\geodatabase.sde\featureclass2 does not exist or is not supported
Failed to execute (GetCount).

If I change the fc variable to include the schema name e.g. fc = "differentschema.featureclass2" then it works no problem.

Is there a way to quickly determine the schema name for a feature class from arcpy so that the tools work correctly?

I tried doing an arcpy.ListFeatureClasses("*featureclass2") however this adds a lot of time to the process (the gdb has a lot of feature classes, many of which are in Feature Datasets also), so would like to avoid that if possible.


texts: Mac os x 10.9.4 run geoserver binary in background

texts:
I have had problems running Geoserver 2.5.2 as a WAR package in Tomcat 7 on OS X 10.9. Tomcat runs but when I add Geoserver, Tomcat crashes. This happened after the OS X 10.9.4 update. Well, I didn't bother to use more time to fix this since the Geoserver OS-independant binary seems to work fine.

Only problem I have now is how to have Geoserver run in the background and upon startup. I am used to run programs in the background like this:

/usr/local/geoserver/bin/startup.sh &

This starts geoserver but when I close the terminal window, geoserver stops. This may be a stupid question but I would be glad for any help. Also for a startup script.


texts:
I have UTM coordinates data XY and Z which shows the presence of mountain. For Z, I have just binary data.

1 = presence of mountains, and
0 = absence of mountains (presence of sediments in other words).

I am wondering that how can I plot it in ArcMap?


texts: Constructing CON statement in Raster Calculator of ArcGIS Spatial Analyst

texts:
I have this expression that works in the QGIS Raster Calculator but not ArcGIS Raster Calculator.
How can I reconstruct it?
It is a binary raster that I am working with.
("Outputs\Model Fin" == 1 & "Outputs\Obs Fin" == 1) * 1 + ("Outputs\Model Fin" == 1 & "Outputs\Obs Fin" == 0) * 2 + ("Outputs\Model Fin" == 0 & "Outputs\Obs Fin" == 1) * 3 + ("Outputs\Model Fin" == 0 & "Outputs\Obs Fin" == 0) * 4


texts: Representing dataset with long field name using graduated color?

texts:
I have some questions on symbology of dataset.

I have 47 recordset ranging from -1 to 10. And some of value have zero.

And I want to represent my dataset with graduated color. It usually works well.

But something weird this time. Different graduated color representation is showing whenever i regenerate the dataset.

So map is changing everytime I have clicked the button on the map and It is closing suddenly.

Why does it happen?


texts: Dealing with duplicate locations kriging - a method for replacing with a single location taking the mean of the duplicates? (R)

texts:
I understand that the existence of duplicate points (locations) when Kriging will result in a singular covariance matrix. The only way I've encountered people dealing with this is to keep the first of any duplicate point and removing the rest.

But it seems rather wasteful to me to arbitrarily throw away data like that. So I am trying to write an R script which will take any duplicate locations in my data and replace with a single location and the mean of the variables.

The data I'm working with can be downloaded here (https://www.dropbox.com/sh/xnwp3zz5abnilyo/AABRVJZ0kTmWk0T9Fcp4-bVSa?dl=0/).

And here's how I attempted this:

library(rgdal)
library(sp)
library(maptools)

 #load data

hs1<- readOGR (".", "Hollicombe_S1_L1-5_A1.2")

#remove columns we're not interested in

hs1<- subset(hs1, select = -c(1:16, 23:24))

So I start with hs1 - a SPDF with 552 obs and 6 variables...

#check for duplicate location (present if lengths differ)
length(hs1@coords)
[1] 1104
length(unique(hs1@coords))
[1] 730

This confirms there are duplicates. So now my attempt:

hs1.d <- hs1[duplicated(hs1@coords),] # creates new SPDF with only duplicated locations
hs1.u <- hs1[!duplicated(hs1@coords),] # creates new SPDF with only unique locations

# coerce duplicated locations SPDF to an ordinary data frame

hs1.md<- as.data.frame(hs1.d)

# combine the X&Y into a single "location"
hs1.md <- within(hs1.md,
  Location <- paste(coords.x1, coords.x2, sep = ","))

# aggregate duplicate locations and calculate a mean value for each

means_by_location<-  aggregate (cbind(BioArea,BioVolume,MeanBioHei,MaxBioheig,PerArIn, PerVolIn)~Location,  hs1.md, mean)

#split location back to X&Y

lat_long <- strsplit(means_by_location$Location, ",")
means_by_location$coords.x1 <- sapply(lat_long, function(x) x[1]) #adds X data back
means_by_location$coords.x2 <- sapply(lat_long, function(x) x[2])#adds Y data back
means_by_location$coords.x1 <- as.numeric (means_by_location$coords.x1) #converts to numeric
means_by_location$coords.x2 <- as.numeric (means_by_location$coords.x2)#converts to numeric

# add spatial information back in to create SPDF

coordinates(means_by_location) = ~coords.x1+coords.x2 # adds the locations
proj4string(means_by_location) = CRS(proj4string(hs1)) # sets the CRS

# hs1.md as SPDF containing single rows for previously duplicated locations
# with mean values for each variable

hs1.md <- subset(means_by_location, select = -(1))

#merge hs1.md and hs1.u to create new SPDF without duplicates

hs1 <- spRbind (hs1.u, hs1.md)

So hs1 is now a SPDF with 543 obs (i.e. 9 observations have been removed).

But there still remain duplicate locations and the number of unique locations remains the same :

length(hs1@coords) # total number of locations
[1] 1086

length(unique(hs1@coords)) #number of unique locations
[1] 730

Can anyone spot where I have gone wrong? Or does anyone else have a similar method for dealing with duplicates?


texts: How to remove selected elements from a QGIS legend?

texts:
I'm working with a QGIS project with a number of layers. When assembling a printout all visible layers are automatically added to the legend:

I would like to exercise some control over the legend elements are remove some of the entries:

However, as the entries correspond to the upper levels in the hierarchy, deleting them will remove the keys:

  1 - 20000
   20000 - 40000
   ...

I'm working in Qgis 2.8. I'm not interested in modifying my layer structure but only in changing the legend elements.


texts:
I am a novice to the field of GIS so please forgive my ignorance.

I am required to plot a shapefile (.shp) on a basemap which is an openlayers map. The shapefile is published in GeoServer as a wms layer.

How do I extract the legends information from the shapefile so that I can display it on the webpage?


texts: Mapping multiple historical routes from lists of points

texts:
I'm working with a spreadsheet of data culled from historical archives, in which I've got a list of travelers, each with a start location and an end location (lat/long). Some also have known intermediate stopping points, but not all. What I want to do in qgis or arcgis is map each of these people's travels as a separate line/route made up of two to five connected points. Is it possible to do this? Am I missing an obvious solution?

I'm pushing up against the limits of my autodidact humanities-person GIS skills. I can figure out how to create a layer (from a CSV file) marking one point per person, but beyond that I am stuck. There are a couple of hundred people involved, meaning drawing the routes manually is really not optimal.

I can kludge together more or less the right result in Palladio, the Stanford digital humanities app (http://hdlab.stanford.edu/palladio/), where it ends up looking something like this:
https://imgur.com/CloRIVw

But Palladio is rather awkward/limited for this purpose, and there's no way to join together line segments belonging to the same person's journey, or to export or save the results for use in other GIS software.


texts: Mosaicing non-overlapping rasters in QGIS

texts:
I have a set of georeferenced rasters that I want to combine in a seamless layer in QGIS. I believe what I have to do is clip them by mask layer and then build a virtual raster. There is an issue however. The original maps that I have scanned have zero overlap, so each one starts exactly where the other one ends. What do you think would be the best course of action here:

include the neatlines in the clipping mask - then my eventual map will have black lines along the seams
exclude the neatlines - then there will be gaps up to 10 pixels wide - is there any good way to fill them?
any other suggestions? maybe I am approaching it from entirely wrong angle?

By the way, is there by any chance a tool to trace the neatlines automatically instead of painstakingly drawing the clipping polygons manually?


texts: QGIS Map tools and Map canvas giving QString errors

texts:
I'm trying to learn more about using QGIS Map Canvas widgets and creating custom map tools, so the information from the PyQGIS cookbook is very helpful. After trying to experiment with mywnd.py from this page: https://www.qgis.org/en/docs/pyqgis_developer_cookbook/canvas.html I had a few difficulties.

On first try running the script, I got ImportError: cannot import name QString. To overcome this error I imported like so:

try:
    from PyQt4.QtCore import SIGNAL, Qt, QString
except ImportError:
    QString = str

That allowed me to get through the import. However, after inputting the following lines into the python console as directed by the cookbook, I got NameError: global name 'QString' is not defined

w = mywnd.MyWnd(qgis.utils.iface.activeLayer())


texts: Dissolve: Statistics type "most frequent" value of dissolved Polygons

texts: I want to aggregate my features which share a common boundary in ArcGIS. This is easily done using the Dissolve Tool in ArcGIS. However, I want to have the attributes of the dissolved polygon based on the attribute occuring most often in this case (see below just a descriptive example). I know that I can dissolve my polygons and then do a Spatial Join and there set the Merge Rule to "Mode". But why is this not directly possible in the Dissolve-functionality in the Statistic Type, do I miss something? Doing this two step approach is quite time consuming, especially when you have very big datasets (Spatial Join can compute for long times as we all know).


texts: Map tm_shape using digitized historical map as background

texts:
I am trying to create a map using tm_shape on r

While the process is relatively easy with regular .tif rasters the process seems more complicated when I try to use a raster from a digitized historic map to be the background (historic.tif).

When I load the historic.tif file on QGIS I can then plot my other shapefiles (roads, cities) and make the map. However on r I can't load the file. I get the following message :

This historic.tif is an old map of Liberia that I found in the archives and digitized using QGIS and then saved as a .tif file.

  Error in .local(.Object, ...) :

  Error in .rasterObjectFromFile(x, band = band, objecttype =
  "RasterLayer",  :    Cannot create a RasterLayer object from this
  file. (file does not exist)

library(sf)
library(raster)
library(dplyr)
library(spData)
library(spDataLarge)
library(rgdal)
library(rgeos)

library(tmap)

historic_map <- raster("~/historic.tif")

Does anyone has a solution to propose?


texts: I'm using Geoserver 2.1.1, Is there way to remove the Layer preview?, actually I'm only using the gmaps web services, I want to remove the KML preview or at least add some sort of authentication.


texts: When I create a strip map of a line the polygons overlap slightly even if I choose a 0% overlap. How Can I make the strip map completely flush without manually going through each frame and clipping them?


texts: I'm new to Python/VB Script and am needing some help.  I have tried several different If/Then statements but can't seem to get it straight.  I am hoping for some help.  I have 3 fields.  ZONING, SQFeet_1, and Buildout. Depending on the Value of the ZONING field, will depend on what formula I need to use.  For instance, if ZONING = RS6 then I need to use the formula  (SQFeet_1/6000*.80) or if ZONING = SP then I need to use the formula (SQFeet_1/2000*.80), etc...  How would I write the IF/Then statement to get the results?  Or would I need to use a different statement?


texts: How to get current scale in OpenLayers 3

texts:
I am new with OpenLayers 3. I'm working with the projection EPSG:4326.

I can get current resolution from view.getResolution(), but I am having trouble to get current scale of map. Please, give me a solution as soon as possible.


texts:
I am supposed to copy a shapefile into an existing database (*.gdb) using QGIS.
There is no problem opening the database, neither when pasting the elements into it.

However, once I try to save the edits, I get the following error :

Could not commit changes to layer CDN_GDB_Nozon_v28 CDN_CE_CCO

Errors: ERROR: 155 feature(s) not added.

Provider errors:

OGR error creating feature -312: Failed at writing Row to Table in CreateFeature. (The spatial index grid size is invalid.)

How can I solve it?


texts: Leaflet GeoJSON and JQuery Autocomplete issue

texts:
I am using jquery autocomplete to search polygons. the JQuery UI uses an array as a data source.  I created 4 arrays, the first 3 for each of the three GeoJSON layers. They load okay. Then I concat them into the 4th array, which I hoped to be my data source for the Jquery UI.

However it fails to work. The GeoJSON is async and appears to load their array after the JQuery UI loads. So the autocomplete doesn't get the array data source.

What I need is to have all three layers arrays complete and combined, then to fire off the JQuery UI using this as it's source. I looked at chaining promises but that's beyond me.

var arr = [];
var arr1 = [];
var arr2 = [];
var arr3 = [];

var senate = L.geoJSON(null, {onEachFeature: forEachFeature, style: style2}).addTo(map);
$.getJSON(url, function(data) {
senate.addData(data);

for (i = 0; i < data.features.length; i++) {
    arr1.push({name:data.features[i].properties.NAME, Dist:data.features[i].properties.DISTRICT, Data:"senate"});
}
    console.log(arr1.length +" records in senate");
});

var assembly = L.geoJSON(null, {onEachFeature: forEachFeature, style: style2});
$.getJSON(url2, function(data) {
    assembly.addData(data);

    for (i = 0; i < data.features.length; i++) {
        arr2.push({name:data.features[i].properties.NAME, Dist:data.features[i].properties.DISTRICT, Data:"assembly"});
    }
    console.log(arr2.length +" records in assembly");
});

var congress = L.geoJSON(null, {onEachFeature: forEachFeature, style: style2});
$.getJSON(url3, function(data) {
    congress.addData(data);

    for (i = 0; i < data.features.length; i++) {
        arr3.push({name:data.features[i].properties.NAME, Dist:data.features[i].properties.DISTRICT, Data:"congress"});
    }
    console.log(arr3.length +" records in congress");
});

console.log(" Records in assembly 240");
var arr = arr1.concat(arr2, arr3);
console.log(arr.length +" records combined should be 240");

// Jquery fails as no data in the source.    //$( function() {
// $(document).ready(function(){     //$(window).bind("load", function() {
document.addEventListener("DOMContentLoaded", function(event){
console.log(arr.length +' Data in array - jquery');
var availableTags = [];

for (i = 0; i < arr.length; i++) {

My console.logs shows the order:

Records in assembly 240
0 records combined should be 240  << before JQuery
0 Data in array - jquery  <<shows Autocomplete has no data
27 records in congress  << No GeoJSON data loading
63 records in senate
150 records in assembly


texts:
I want to use EPSG:25832 as the default projection for all my layers. The projection is not working correctly because the Position (center point) should be in Germany but I end up in France. When I want to pan to Germany I get an "InvalidStateError". I did the following:

loaded the latest proj4j Version with

<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.3.12/proj4-src.js"></script>

and

<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.3.12/proj4.js"></script>

defined the projection (parameter from http://spatialreference.org/ref/epsg/etrs89-utm-zone-32n/):

 var myProjectionName = 'EPSG:25832'; proj4.defs(myProjectionName,
  "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs");

created the projection:

var myProjection = new ol.proj.Projection({
    code: myProjectionName,
    extent: [265948.8191, 6421521.2254, 677786.3629, 7288831.7014],
    units: 'm'
});
ol.proj.addProjection(myProjection);

I got the extent parameter also from http://spatialreference.org/ref/epsg/etrs89-utm-zone-32n/

created the map:

var map = new ol.Map({
  layers: layers,
  target: 'map',
  view: new ol.View({
  projection: myProjection, // Is this necessary?
  center: ol.proj.transform([7.20,51.67], 'EPSG:4326', myProjectionName),
  zoom: 10
  })
});

I saw different approaches on how the map was created (with and without the projection parameter). Is the projection needed or is it sufficient to add the "ol.proj.transform" at the Center parameter?

Before I used

center: ol.proj.transform([7.20,51.67], 'EPSG:4326', 'EPSG:3857'),

which worked fine for displaying the data. Now I am working with a WFS-T and it is necessary to use EPSG:25832.
What is wrong with my projection?


texts: Seeking QGIS equivalent to Focal Statistics of ArcGIS Spatial Analyst?

texts:
I'm attempting to use QGIS v2.0.1 to accomplish the same tasks that I can do using ArcMap v10.1 and I'm running into some trouble finding equivalent tools. In ArcMap I am using the Focal Statistics (mean) Spatial Analyst tool and I do not know if GRASS or SAGA has an equivalent in QGIS.

Does anyone know what the tool would be called if it exists?


texts: Changing style of a LineString in OpenLayers 6

texts:
I am using OpenLayers 6.5 and I have a linestring added as a map layer. I am having trouble accessing the feature (linestring) in order to change the style .
Here is how I create it
var lineString = new ol.layer.Vector({
    id: 'MarkersTrail' + myid,
    source: new ol.source.Vector({
        features: [
            new ol.Feature({
                id: myid,
                geometry: new ol.geom.LineString(points),
                name: myid
            })
       ]
   }),
   style: new ol.style.Style({
       stroke : new ol.style.Stroke({
           color: '#0000ff',
           width: 5
       })
   })

  });
   map.addLayer(lineString);
   lineString.setVisible(false);

When I then try to access it my feature is null?
layer is determined by map.getLayers().
var feature = layer.getSource().getFeatureById(myid);

What am I doing wrong?


texts: mapbox gl- Toggle highlighting of specific feature of a layer in map style

texts:
I tried it with "Setpaintproperty" and "data driven functions" but I dont know the exact syntax..

    map.on('mousemove', function(e) {
    var states = map.queryRenderedFeatures(e.point, {
     layers: ['n_Region6 Municipals']
    });
  map.setPaintProperty('n_Region6 Municipals', 'fill-color','
  {"property":"NAME_2","type":"category","stops": +
  [[states[0].properties.NAME_2 +,"blue"]]}');// im not sureif this is
                                                 allowed

I see some Instructions on the web using filter but it makes you add 2 or more layers and I dont want that.I only what to use the existing layer in the map style.is that posssible in Mapbox GL?how to code it simply?


texts: ArcGIS field calculator Python generate random point names from attribute table

texts:
I am new to arcpy and am looking to generate 15 random point values, which are names from an existing attribute column into a new column. I must use Python in the field calculator. My existing column name is: CityNames which includes 350 cites, and my new column name is: randomCities.

 def RandomPts(randomCities):
 import random
 'return 15 random points'
 return random.choice(CityNames,15)


texts:
I am using MapGuide + Fusion for the basic Sheboygan-google example. I need to project the Sheboygan.MapDefinition to SheboyganMercator.MapDefinition so I can correctly display the example over google maps. However the WKT in the example to project it is not validated by maestro. Also maestro does not support projecting to Pseudo-Mercator at this time.

How can I create the projection manually?


texts: Avoid intersection is missing QGIS 2.10, 2.8.2

texts:
I need to digitize a lot of polygons and I want to ease my work with the 'avoid intersection' button so that the new polygons do not overlap. In the manual I found (for version 2.6) that the button is hidden in 'snap options', but I can't find it.

How do I activate 'avoid intersection'?


texts: Creating a Landsat Mosaic from frames in different UTM Zones

texts:
As part of my current project I have calculated NDVI change over time for four connected Landsat frames. Now that I have this done I was planning on mosaicking them together so I could conduct some analysis (mainly to relate NDVI change to landcover type, topography and other landscape variables).

The issue I am having is that the southern most frame is in UTM Zone 12N while the three frames north of it are in UTM Zone 11N. This causes the mosaic to shift pixels around a bit. I created a mosaic in ENVI of just the three northern frames and it worked perfectly, with pixels staying where they should be and values remaining the same.

I am looking for advice on how to proceed.

What is the best way to mosaic all four frames together that has minimal impact on my subsequent analysis (eg. not changing or interpolating new values)?

Or should I run my analysis twice (one for the UTM Zone 12N frames and another for the UTM Zone 11N frame) in order to insure the quality of my data and avoid any possible warping that reprojecting could introduce?

I would like to mosaic all four frames together if I could, but not at the cost of lowering the quality of my data.


texts: How do I go from the WKT results of Calculate UTM Zone (Cartography) to a spatial reference object?

texts:
Calculate UTM Zone (Cartography) assigns a UTM zone to each feature within a feature class as a field value. The value returned is a 'UTM spatial reference string', which is a long string that looks something like this:

PROJCS["GCS North American 1983 UTM Zone 10S (Calculated)",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-123.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

I'd like to use this information to create a spatial reference object in the projection indicated in the UTM spatial reference string. I'm not sure how to go about this, however. It looks like parsing out a value in between the brackets after PROJCS would get me close. However, given the different format of the string (spaces vs underscores, etc.) I'm not sure how reliably I'd be able to get it to work anyways. Is there a better way to get a spatial reference object that each feature falls within?

Maybe there's a better way to get UTM zone per feature, but I would like the code to stand alone and not require a UTM Zone feature class to function.


texts:
I downloaded the Leaflet Routing Machine from https://github.com/perliedman/leaflet-routing-machine. Unzipped and opened the index.html in my browser to see how it looks, and got an error as shown in the image below:

ERR_FILE_NOT_FOUND file://router.project-osrm.org/viaroute?instructions=true&loc=57.74,11.94&loc=57.6792,11.949&jsonp=_l_routing_machine_0

index.html

<!DOCTYPE html>

    <meta charset="utf-8" />
    <title>Leaflet OSRM Example</title>
    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.2/leaflet.css" />
    <link rel="stylesheet" href="../css/leaflet-routing-machine.css" />
    <link rel="stylesheet" href="index.css" />


texts: Merging two point data sets with artificial grid cells in QGIS

texts:
I have two datasets. The SCAD has conflict data with the exact latitude and longitude, the other one has the exact coordinates for foreign aid. My research question is, whether foreign aid decreases the likelihood of a conflict. Therefore, I would like to have both data in one dataset and the unit of analysis should be 1km x 1km artificial grid cells.
In the end I would like to have artificial grid cells as unit of analysis. With the exact coordinates I cannot do any data analysis.


texts: Properly Converting from UTM to Lat/Long with gdalwarp

texts:
I have a file in GeoTiff format in UTM coordinates.  I convert the file to LatLong, then back to UTM, and the coordinates aren't what I started with.  I have the printout of the process below.

gdalsrsinfo 6b6.tif produces the following:

PROJ.4 : '+proj=utm +zone=16 +datum=NAD83 +units=m +no_defs '

OGC WKT :
PROJCS["NAD_1983_UTM_Zone_16N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26916"]]

gdalwarp -t_srs '+proj=longlat +datum=WGS84 +no_defs' 6b6.tif output.tif

gdalinfo output.tif
Driver: GTiff/GeoTIFF
Files: output.tif
Size is 1384, 1379
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-86.750459052801645,40.250430422810176)
Pixel Size = (0.000091263992370,-0.000091263992370)
Metadata:
  AREA_OR_POINT=Area
  BandDefinitionKeyword=NONE
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -86.7504591,  40.2504304) ( 86d45' 1.65"W, 40d15' 1.55"N)
Lower Left  ( -86.7504591,  40.1245774) ( 86d45' 1.65"W, 40d 7'28.48"N)
Upper Right ( -86.6241497,  40.2504304) ( 86d37'26.94"W, 40d15' 1.55"N)
Lower Right ( -86.6241497,  40.1245774) ( 86d37'26.94"W, 40d 7'28.48"N)
Center      ( -86.6873044,  40.1875039) ( 86d41'14.30"W, 40d11'15.01"N)
Band 1 Block=1384x1 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38
  Metadata:
    RepresentationType=ATHEMATIC

and then converting back, gdalwarp -t_srs EPSG:26916 output.tif back.tif

gdalinfo back.tif
Driver: GTiff/GeoTIFF
Files: back.tif
Size is 1197, 1552
Coordinate System is:
PROJCS["NAD83 / UTM zone 16N",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-87],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","26916"]]
Origin = (521222.733228680794127,4455620.891041336581111)
Pixel Size = (9.022681340919917,-9.022681340919917)
Metadata:
  AREA_OR_POINT=Area
  BandDefinitionKeyword=NONE
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  521222.733, 4455620.891) ( 86d45' 1.65"W, 40d15' 2.78"N)
Lower Left  (  521222.733, 4441617.690) ( 86d45' 3.31"W, 40d 7'28.60"N)
Upper Right (  532022.883, 4455620.891) ( 86d37'24.49"W, 40d15' 1.54"N)
Lower Right (  532022.883, 4441617.690) ( 86d37'27.00"W, 40d 7'27.37"N)
Center      (  526622.808, 4448619.290) ( 86d41'14.11"W, 40d11'15.14"N)
Band 1 Block=1197x1 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38
  Metadata:
    RepresentationType=ATHEMATIC


texts: Openlayers-3 - Extracting data for styling from my GeoJSON

texts:
I am trying to understand how to get data from my GeoJSON while rendering the map, so I can style every polygon/point differently.

What i'd like is to create style-groups with the property named "Functie" from my GeoJSON.
It just doesn't seem to work -> there is no styling applied to the element and the element does not appears on the map. If I don't add any styling, the points do appear with the basic styling.

My GeoJSON:

var puntobject = {
    "type": "Feature",
    "properties": {
        "property_a": 4,
        "property_d": 123,
        "property_c": 1242345,
        "property_d": 55151,
    },
    "geometry": {
        "type": "Point",
        "coordinates": [5.533485496901196, 51.609474175170639, 0.0]
    }
}

My code so far, based on an example. This code does not get the specific property, but the type. But this doesn't seem to work either -> no styling or visible elements?

var styleFunction = function(feature,geom) {
console.log(feature.getType());
return styles[feature.getGeometry().getType()];
};

var styles = {
    'Feature': [new ol.style.Style({
        image: circle
    })],
    'LineString': [new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'green',
            width: 1
        })
    })]
};

var puntObjectenSource = new ol.source.Vector ({
  features: (new ol.format.GeoJSON()).readFeatures(puntobject)
});
//Puntobjecten style
var puntObjectenStyle = new ol.style.Style ({
   image:circle
});
//Puntobjecten layer
var puntObjecten = new ol.layer.Vector ({
  source: puntObjectenSource,
  style: styleFunction
});

And the final javascript in my html:

   var map = new ol.Map({
      target: 'map',
      layers: [
        new ol.layer.Tile({
          //source: new ol.source.MapQuest({layer: 'sat'})
          source: new ol.source.OSM({layer:'hot'})
        }),
        puntObjecten
      ],
      view: new ol.View({
        center: ol.proj.transform(lonlat, 'EPSG:4326', 'EPSG:3857'),
        zoom: 15
      })
    });

I am trying to get the properties out of my JSON -> ON LOAD. This way I can add the styling on rendering per element and add some colors.
How can I get this done?

P.S.
I am new to Openlayers and having HUGE difficulties trying to understand the documentation on the website. Any tips on how to understand the documentation are welcome.


texts:
layer_gesamt = new OpenLayers.Layer.WMS(
            "Geoserver layers - Tiled", "http://mygeoserver/xyz/wms",
                {
                    LAYERS: 'myLayer',
                    STYLES: '',
                    format:  'image/png',
                    transparent: true,
                    tiled: true,
                    tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom,
                },
                {
                    buffer: 0,
                    displayOutsideMaxExtent: false,
                    isBaseLayer: true,
                    yx : {epsg : true}
                }

          );
  layer_gesamt.setVisibility(false);

  map.addLayers([
     layer_gesamt,
     layer_selectable,
                 ]);
  map.setCenter(config.init_center.getCenterLonLat(),config.init_zoom);

I want the layer_gesamt in the example NOT to be displayed,  but to be loaded. That's why I set  layer_gesamt.setVisibility(false);. However, the layer is still displayed. What can I do?


texts: Can I edit shapefile while using a definition query in QGIS?

texts:
Is it possible to edit shapes while using a "definition query" or "query builder" in QGIS? I'm trying to do this in the same way I used to in ArcGIS.

Thanks!


texts: Showing traffic signs at different scales using ArcGIS Desktop?

texts:
I would like to have a single point feature class containing traffic signs but the challenge is a cartographic one; displaying these signs correctly at varying scales.

If I create several signs in their actual geographic location, say around a road intersection, they will only show correctly at that scale. As soon as I zoom out, the points bunch-up and show as a cluster on top of that intersection.  My current workaround is to place the signs well away of from the road but in the general location so that as i zoom in and out the sighs are shown more less in the approximate location.  I also limit the visibility to a small range of scales.  Not ideal.

I would like to be able to zoom out and have the points be "pushed out" to a scale-depended but conceptually correct location. The complexity of this is compounded by the fact that the signs need to be pushed out into certain directions based on which side of the road they are on, while ensuring they are not moved across an intersecting nearby road, etc.

Can this be done without having to create a separate layer for each scale?

The target audience for this product are internal web GIS users. The matter is further complicated by the fact that I want the users to be able to edit, create and delete these signs so this limits the use of custom programming and tools; the data and its symbology needs to be dynamic and be able to be exported as a feature service for use on web GIS.

The actual geographic GPS location of the sign is not important, as long as it is conceptually correct at varying scales.

I am open to ideas even if it means a completely different strategy (within the parameter specified below)

ArcGIS Standard, Spatial Analyst


texts:
I need to update a table of points with the FIP code of the nearest county polygon from a separate table.  I've managed to update the column for those points that Intersect a county polygon with this:

update points p set fips = c.fips
from counties c where ST_Intersects(c.the_geom, p.the_geom)

but I've got about 60 coastal points that, because of generalization in the county polygons, don't overlap a polygon feature.  I need to update the fips of these with the NEAREST county polygon's fips code, and I've been futzing around with the <-> Indexed NN, but I just can't get it to work for me.

So, where p.fips IS NULL I need to set p.fips = c.fips of the nearest county feature.


texts: Does anyone know of any good sources for US data (counties, congressional districts, etc) in GeoJSON format? Has anyone converted public-domain shapefiles/sources to GeoJSON that they would be willing to share?


texts:
Newbie here confused by my first experience in QGIS

I have a handful of lines defining a bike route.

Originally they were saved in a .gml file and later translated into .kml so they could be viewed in GE.

I can open both .gml and .kml as layers in QGIS, but neither allows me to edit the lines

If I save either the .gml or .kml layers as a Shapefile, then load that shapefile as a QGIS layer, then I can edit the lines.

I have also tried loading the .kml file, making a new blank shapefile layer, selecting and copying the lines on the .kml layer and tried to paste onto the blank shapefile layer, which also was not possible.

I don't understand if lines are understood by QGIS on kml layer well enough to be saved as Shapefile, why they are not editable.   I feel I must be missing something basic that explains this and perhaps it has other implications that I should understand.


texts: Downloading high resolution Google Earth Images

texts:
How can I download historical google earth images for years like 2015 and 2020?
I cannot find date filter in universal map downloader.